--- /dev/null
+#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");
+}
+
projects / strong_simulation_stabilizer_rank.git / blobdiff