X-Git-Url: https://s3miclassical.com/gitweb/?p=strong_simulation_stabilizer_rank.git;a=blobdiff_plain;f=exponentialsum.c;fp=exponentialsum.c;h=0e0bd6e8cd890ec0ade0f6496abef63872405725;hp=0000000000000000000000000000000000000000;hb=f36d60f8c85e0879a7d2b52e816638d2b4fdf86a;hpb=a85d4d567fc589503c5e263e1520295c8d433a45

diff --git a/exponentialsum.c b/exponentialsum.c
new file mode 100644
index 0000000..0e0bd6e
--- /dev/null
+++ b/exponentialsum.c
@@ -0,0 +1,415 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "exponentialsum.h"
+
+// linked list
+struct Node {
+  int data;
+  struct Node* next;
+};
+
+void printLinkedList(struct Node *node);
+void append(struct Node** head_ref, int new_data);
+void append2(struct Node** head_ref, struct Node** last_ref, int new_data);
+void freeList(struct Node** head_ref);
+
+complex double exponentialsum(int *k, int *Q, int *D, int **J) {
+
+  int a, b, c, d;
+
+  int tmp;
+
+  struct Node* setS = NULL;
+  struct Node* setWalker = NULL;
+  struct Node* setWalker2 = NULL;
+
+  struct Node* setDa = NULL; // first element in \mathcal D
+  struct Node* setDb = NULL; // second element in \mathcal D (sets should have the same size)
+
+  int setSleftover = -1;  // -1 means nothing left over
+
+  int** R;  // basis change matrix R and R^T
+  int** RT;
+  R = calloc(*k, sizeof(int*));
+  RT = calloc(*k, sizeof(int*));
+  // initially set it to the k-dimensional identity matrix
+  for(a=0; a<*k; a++) {
+    R[a] = calloc(*k, sizeof(int));
+    R[a][a] = 1;
+    RT[a] = calloc(*k, sizeof(int));
+    RT[a][a] = 1;
+  }
+
+  int* tmpVec; // temporary vector used in Eq. 49 to update D
+  int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
+  tmpMatrix = calloc(*k, sizeof(int*));
+  for(a=0; a<*k; a++)
+    tmpMatrix[a] = calloc(*k, sizeof(int));
+
+  for(a=0; a<*k; a++) {
+    if(D[a] == 2 || D[a] == 6) {
+      append(&setS, a);
+    }
+  }
+
+  // Choose 'a' to be the first element in the linked list setS; a=setS->data
+
+  if(setS != NULL) {
+    setWalker = setS;
+    while(setWalker->next != NULL) {
+      setWalker = setWalker->next;
+      R[setWalker->data][setS->data] += 1;
+    }
+
+    tmpVec = calloc(*k, sizeof(int));
+    for(c=0;c<*k;c++) {
+      tmp = 0;
+      /* //tmpVec[c] = 0; */
+      /* for(a=0;a<*k;a++) // don't confuse this for-loop a with 'a'=setS->data */
+      /* 	tmpVec[c] = tmpVec[c] + R[c][a]*D[a]; */
+      /* for(a=0;a<*k;a++) */
+      /* 	for(b=a+1;b<*k;b++) */
+      /* 	  tmpVec[c] = tmpVec[c] + J[a][b]*R[c][a]*R[c][b]; */
+      /* tmpVec[c] = tmpVec[c]%8; */
+      if(c!=setS->data)
+	tmp += R[c][c]*D[c];
+      tmp += R[c][setS->data]*D[setS->data];
+      if(setS->data>c)
+	tmp += J[c][setS->data]*R[c][c]*R[c][setS->data];
+      else if(setS->data<c)
+	tmp += J[setS->data][c]*R[c][setS->data]*R[c][c];
+      tmpVec[c] = tmp%8;
+    }
+    memcpy(D, tmpVec, sizeof(int)*(*k));
+    free(tmpVec);
+
+    // Eq. 50 update of J
+    for(c=0;c<*k;c++) {
+      for(d=0;d<*k; d++) {
+	tmp = 0;
+	if(c!=setS->data) {
+	  if(d!=setS->data) {
+	    tmp += R[c][c]*J[c][d]*R[d][d];
+	    tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
+	    tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
+	  }else {
+	    tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
+	  }
+	} else {
+	  if(d!=setS->data) {
+	    tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
+	  }
+	}
+	tmp += R[c][setS->data]*J[setS->data][setS->data]*R[d][setS->data];
+	tmp %= 8;
+	tmpMatrix[c][d] = tmp;
+      }
+    }
+    for(c=0; c<*k; c++)
+      memcpy(J[c], tmpMatrix[c], sizeof(int)*(*k));
+    /* transp(R,RT,*k,*k); */
+    /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+    /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
+
+    // remove everything but first element from setS
+    //tmp = setS->data;
+    //freeList(&setS);
+    //append(&setS, tmp);
+    //setS->next = NULL;
+    setSleftover = setS->data;
+  }
+
+  struct Node* setE = NULL;
+
+  for(a=0;a<*k;a++) {
+    //if((setS != NULL && a != setS->data) || setS == NULL)
+    if((setSleftover != -1 && a != setSleftover) || setSleftover == -1)
+      append(&setE, a);
+  }
+
+  struct Node* setM = NULL;
+  int setMCount = 0; // number of elements in M
+
+  int r = 0;
+
+  while(setE != NULL) {
+    // let 'a' be setE->data (the first element in setE)
+    struct Node* setK = NULL;
+    setWalker = setE;
+    while(setWalker->next != NULL) {
+      setWalker = setWalker->next;
+      if(J[setE->data][setWalker->data] == 4)
+	append(&setK, setWalker->data);
+    }
+
+    if(setK == NULL) { // Found new monomer {'a'}
+      append(&setM, setE->data);
+      setMCount++;
+
+      if(setE->next != NULL) {
+	setE->data = setE->next->data;
+	setE->next = setE->next->next;
+      }
+      else setE = NULL;
+    } else {
+      // let 'b' be setK->data (the first element in setK)
+      setWalker = setE->next;
+      // reset 'R' to the k-dimensional identity matrix
+      for(a=0; a<*k; a++) {
+	memset(R[a], 0, sizeof(int)*(*k)); 
+	//for(b=0; b<*k; b++)
+	//  R[a][b] = 0;
+	R[a][a] = 1;
+      }
+	  
+      while(setWalker != NULL) {
+	if(setWalker->data == setK->data) {
+	  setWalker = setWalker->next;
+	  continue;  // 'c' \in E minus {'a', 'b'} where 'a'=setE->data, 'b'=setK->data and 'c'=setWalker->data
+	}
+	R[setWalker->data][setK->data] = J[setE->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
+	R[setWalker->data][setE->data] = J[setK->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
+	setWalker = setWalker->next;
+      }
+      // update D and J
+      tmpVec = calloc(*k, sizeof(int));
+      for(c=0;c<*k;c++) {
+	tmp = 0;
+	//tmpVec[c] = 0;
+	/* for(a=0;a<*k;a++) {// don't confuse this for-loop a with 'a'=setS->data */
+	/*   tmp+= R[c][a]*D[a]; */
+	/* //for(a=0;a<*k;a++) */
+	/*   for(b=a+1;b<*k;b++) */
+	/*     tmp += J[a][b]*R[c][a]*R[c][b]; */
+	/* } */
+	if(c!=setK->data && c!=setE->data)
+	  tmp += R[c][c]*D[c];
+	tmp += R[c][setK->data]*D[setK->data];
+	tmp += R[c][setE->data]*D[setE->data];
+	if(setK->data>setE->data)
+	  tmp += J[setE->data][setK->data]*R[c][setE->data]*R[c][setK->data];
+	else if(setE->data>setK->data)
+	  tmp += J[setK->data][setE->data]*R[c][setK->data]*R[c][setE->data];
+	if(setK->data>c)
+	  tmp += J[c][setK->data]*R[c][c]*R[c][setK->data];
+	else if(setK->data<c)
+	  tmp += J[setK->data][c]*R[c][setK->data]*R[c][c];
+	if(setE->data>c)
+	  tmp += J[c][setE->data]*R[c][c]*R[c][setE->data];
+	else if(setE->data<c)
+	  tmp += J[setE->data][c]*R[c][setE->data]*R[c][c];
+	tmpVec[c] = tmp%8;
+      }
+      memcpy(D, tmpVec, sizeof(int)*(*k));
+      free(tmpVec);
+
+      // Eq. 50 update of J
+      for(c=0;c<*k;c++) {
+	for(d=0;d<*k; d++) {
+	  tmp = 0;
+	  if(c!=setK->data && c!=setE->data) {
+	    if(d!=setK->data && d!=setE->data) {
+	      tmp += R[c][c]*J[c][d]*R[d][d];
+	      tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
+	      tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
+	      tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
+	      tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
+	    }else {
+	      tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
+	      tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
+	    }
+	  } else {
+	    if(d!=setK->data && d!=setE->data) {
+	      tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
+	      tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
+	    }
+	  }
+	  tmp += R[c][setK->data]*J[setK->data][setK->data]*R[d][setK->data];
+	  tmp += R[c][setE->data]*J[setE->data][setK->data]*R[d][setK->data];
+	  tmp += R[c][setK->data]*J[setK->data][setE->data]*R[d][setE->data];
+	  tmp += R[c][setE->data]*J[setE->data][setE->data]*R[d][setE->data];
+	  tmp %= 8;
+	  tmpMatrix[c][d] = tmp;
+	}
+      }
+      J = tmpMatrix;
+      /* transp(R,RT,*k,*k); */
+      /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+      /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
+
+      // Now {'a','b'} form a dimer
+      r++;
+      append(&setDa,setE->data);
+      append(&setDb,setK->data);
+
+      // E = E minus {'a', 'b'}
+      setWalker = setE; // 'a'=setE->data
+      while(setWalker != NULL) {
+	if(setWalker->next->data == setK->data) {
+	  free(setWalker->next);
+	  setWalker->next = setWalker->next->next; // delete 'b'
+	  break;
+	}
+	setWalker = setWalker->next;
+      }
+      setWalker = setE;
+      setE = setE->next; // delete 'a'
+      free(setWalker);
+    }
+
+    freeList(&setK);
+  }
+
+  complex double W = 0.0 + 0.0*I;
+  
+//if(setS == NULL) {
+  if(setSleftover == -1) {
+    W = cexp(I*PI*0.25*(*Q));
+    setWalker = setM;
+    for(a=0; a<setMCount; a++) {
+      W *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
+      setWalker = setWalker->next;
+    }
+    setWalker = setDa;
+    setWalker2 = setDb;
+    for(a=0; a<r; a++) {
+      W *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data])) + cexp(I*PI*0.25*(D[setWalker2->data])) - cexp(I*PI*0.25*(D[setWalker->data] + D[setWalker2->data])));
+      setWalker = setWalker->next;
+      setWalker2 = setWalker2->next;
+    }
+  } else {
+    complex double W0 = 0.0 + 0.0*I;
+    complex double W1 = 0.0 + 0.0*I;
+    
+    W0 = cexp(I*PI*0.25*(*Q));
+    setWalker = setM;
+    for(a=0; a<setMCount; a++) {
+      W0 *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
+      setWalker = setWalker->next;
+    }
+    setWalker = setDa;
+    setWalker2 = setDb;
+    for(a=0; a<r; a++) {
+      W0 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data])) + cexp(I*PI*0.25*(D[setWalker2->data])) - cexp(I*PI*0.25*(D[setWalker->data] + D[setWalker2->data])));
+      setWalker = setWalker->next;
+      setWalker2 = setWalker2->next;
+    }
+    //    W1 = cexp(I*PI*0.25*(*Q+D[setS->data]));
+    W1 = cexp(I*PI*0.25*(*Q+D[setSleftover]));    
+    setWalker = setM;
+    for(a=0; a<setMCount; a++) {
+      //      W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setS->data])));
+      W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setSleftover])));
+      setWalker = setWalker->next;
+    }
+    setWalker = setDa;
+    setWalker2 = setDb;
+    for(a=0; a<r; a++) {
+      //      W1 *= (1.0 + cexp(I*PI*0.25*(J[setWalker->data][setS->data] + D[setWalker->data])) + cexp(I*PI*0.25*(J[setWalker2->data][setS->data] + D[setWalker2->data])) - cexp(I*PI*0.25*(J[setWalker->data][setS->data] + J[setWalker2->data][setS->data] + D[setWalker->data] + D[setWalker2->data])));
+      W1 *= (1.0 + cexp(I*PI*0.25*(J[setWalker->data][setSleftover] + D[setWalker->data])) + cexp(I*PI*0.25*(J[setWalker2->data][setSleftover] + D[setWalker2->data])) - cexp(I*PI*0.25*(J[setWalker->data][setSleftover] + J[setWalker2->data][setSleftover] + D[setWalker->data] + D[setWalker2->data])));      
+      setWalker = setWalker->next;
+      setWalker2 = setWalker2->next;
+    }
+    W = W0 + W1;
+  }
+
+  deallocate_mem(&tmpMatrix, *k);
+  deallocate_mem(&R, *k);
+  deallocate_mem(&RT, *k);
+
+  freeList(&setDa); freeList(&setDb); freeList(&setM); freeList(&setE); freeList(&setS);
+  
+  return(W);
+  
+}
+
+
+void append(struct Node** head_ref, int new_data) 
+{ 
+    struct Node* new_node = (struct Node*) malloc(sizeof(struct Node)); 
+    struct Node *last = *head_ref;   
+
+    new_node->data = new_data; 
+    new_node->next = NULL;
+
+    // if the linked list is empty, then make the new node as head
+    if (*head_ref == NULL) {
+      *head_ref = new_node; 
+      return; 
+    }   
+       
+    // else traverse till the last node
+    while (last->next != NULL) {
+        last = last->next; 
+    }
+    last->next = new_node;
+    
+    return;     
+} 
+
+/*void freeList(struct Node** head_ref)
+{ 
+       
+  struct Node *second_last_walker = *head_ref;
+  struct Node *walker = *head_ref;   
+
+    // else traverse till the last node
+    if(walker != NULL) {
+      while(walker->next != NULL) {
+	while(walker->next != NULL) {
+	  //printf("!%d\n", walker->data);
+	  second_last_walker = walker;
+	  walker = walker->next;
+	  //printf("%d\n", walker->data);
+	  //printf("!%d\n", second_last_walker->data);
+	}
+	free(walker);
+	second_last_walker->next = NULL;
+	walker = *head_ref;
+	//printf("!!%d\n", second_last_walker->data);
+      }
+    }
+    free(walker);
+    
+    return;     
+    }*/
+
+void freeList(struct Node** head_ref)
+{ 
+       
+  struct Node *walker = *head_ref;
+  struct Node *walker2;   
+
+    // else traverse till the last node
+    if(walker != NULL) {
+      walker2 = walker->next;
+      while(walker2 != NULL) {
+	free(walker);
+	walker = walker2;
+	walker2 = walker->next;
+      }
+    }
+    free(walker);
+    
+    return;     
+}
+
+void printLinkedList(struct Node *node) 
+{
+  if(node == NULL)
+    printf("NULL\n");
+  else {
+    while (node != NULL) 
+      {
+	printf(" %d ", node->data);
+	node = node->next; 
+      }
+  }
+  printf("\n");
+}
+