8 #include "exponentialsum.h"
16 void printLinkedList(struct Node *node);
17 void append(struct Node** head_ref, int new_data);
18 void append2(struct Node** head_ref, struct Node** last_ref, int new_data);
19 void freeList(struct Node** head_ref);
21 complex double exponentialsum(int *k, int *Q, int *D, int **J) {
27 struct Node* setS = NULL;
28 struct Node* setWalker = NULL;
29 struct Node* setWalker2 = NULL;
31 struct Node* setDa = NULL; // first element in \mathcal D
32 struct Node* setDb = NULL; // second element in \mathcal D (sets should have the same size)
34 int setSleftover = -1; // -1 means nothing left over
36 int** R; // basis change matrix R and R^T
38 R = calloc(*k, sizeof(int*));
39 RT = calloc(*k, sizeof(int*));
40 // initially set it to the k-dimensional identity matrix
42 R[a] = calloc(*k, sizeof(int));
44 RT[a] = calloc(*k, sizeof(int));
48 int* tmpVec; // temporary vector used in Eq. 49 to update D
49 int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
50 tmpMatrix = calloc(*k, sizeof(int*));
52 tmpMatrix[a] = calloc(*k, sizeof(int));
55 if(D[a] == 2 || D[a] == 6) {
60 // Choose 'a' to be the first element in the linked list setS; a=setS->data
64 while(setWalker->next != NULL) {
65 setWalker = setWalker->next;
66 R[setWalker->data][setS->data] += 1;
69 tmpVec = calloc(*k, sizeof(int));
72 /* //tmpVec[c] = 0; */
73 /* for(a=0;a<*k;a++) // don't confuse this for-loop a with 'a'=setS->data */
74 /* tmpVec[c] = tmpVec[c] + R[c][a]*D[a]; */
75 /* for(a=0;a<*k;a++) */
76 /* for(b=a+1;b<*k;b++) */
77 /* tmpVec[c] = tmpVec[c] + J[a][b]*R[c][a]*R[c][b]; */
78 /* tmpVec[c] = tmpVec[c]%8; */
81 tmp += R[c][setS->data]*D[setS->data];
83 tmp += J[c][setS->data]*R[c][c]*R[c][setS->data];
85 tmp += J[setS->data][c]*R[c][setS->data]*R[c][c];
88 memcpy(D, tmpVec, sizeof(int)*(*k));
97 tmp += R[c][c]*J[c][d]*R[d][d];
98 tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
99 tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
101 tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
105 tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
108 tmp += R[c][setS->data]*J[setS->data][setS->data]*R[d][setS->data];
110 tmpMatrix[c][d] = tmp;
114 memcpy(J[c], tmpMatrix[c], sizeof(int)*(*k));
115 /* transp(R,RT,*k,*k); */
116 /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
117 /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
119 // remove everything but first element from setS
122 //append(&setS, tmp);
124 setSleftover = setS->data;
127 struct Node* setE = NULL;
130 //if((setS != NULL && a != setS->data) || setS == NULL)
131 if((setSleftover != -1 && a != setSleftover) || setSleftover == -1)
135 struct Node* setM = NULL;
136 int setMCount = 0; // number of elements in M
140 while(setE != NULL) {
141 // let 'a' be setE->data (the first element in setE)
142 struct Node* setK = NULL;
144 while(setWalker->next != NULL) {
145 setWalker = setWalker->next;
146 if(J[setE->data][setWalker->data] == 4)
147 append(&setK, setWalker->data);
150 if(setK == NULL) { // Found new monomer {'a'}
151 append(&setM, setE->data);
154 if(setE->next != NULL) {
155 setE->data = setE->next->data;
156 setE->next = setE->next->next;
160 // let 'b' be setK->data (the first element in setK)
161 setWalker = setE->next;
162 // reset 'R' to the k-dimensional identity matrix
163 for(a=0; a<*k; a++) {
164 memset(R[a], 0, sizeof(int)*(*k));
165 //for(b=0; b<*k; b++)
170 while(setWalker != NULL) {
171 if(setWalker->data == setK->data) {
172 setWalker = setWalker->next;
173 continue; // 'c' \in E minus {'a', 'b'} where 'a'=setE->data, 'b'=setK->data and 'c'=setWalker->data
175 R[setWalker->data][setK->data] = J[setE->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
176 R[setWalker->data][setE->data] = J[setK->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
177 setWalker = setWalker->next;
180 tmpVec = calloc(*k, sizeof(int));
184 /* for(a=0;a<*k;a++) {// don't confuse this for-loop a with 'a'=setS->data */
185 /* tmp+= R[c][a]*D[a]; */
186 /* //for(a=0;a<*k;a++) */
187 /* for(b=a+1;b<*k;b++) */
188 /* tmp += J[a][b]*R[c][a]*R[c][b]; */
190 if(c!=setK->data && c!=setE->data)
192 tmp += R[c][setK->data]*D[setK->data];
193 tmp += R[c][setE->data]*D[setE->data];
194 if(setK->data>setE->data)
195 tmp += J[setE->data][setK->data]*R[c][setE->data]*R[c][setK->data];
196 else if(setE->data>setK->data)
197 tmp += J[setK->data][setE->data]*R[c][setK->data]*R[c][setE->data];
199 tmp += J[c][setK->data]*R[c][c]*R[c][setK->data];
200 else if(setK->data<c)
201 tmp += J[setK->data][c]*R[c][setK->data]*R[c][c];
203 tmp += J[c][setE->data]*R[c][c]*R[c][setE->data];
204 else if(setE->data<c)
205 tmp += J[setE->data][c]*R[c][setE->data]*R[c][c];
208 memcpy(D, tmpVec, sizeof(int)*(*k));
211 // Eq. 50 update of J
215 if(c!=setK->data && c!=setE->data) {
216 if(d!=setK->data && d!=setE->data) {
217 tmp += R[c][c]*J[c][d]*R[d][d];
218 tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
219 tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
220 tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
221 tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
223 tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
224 tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
227 if(d!=setK->data && d!=setE->data) {
228 tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
229 tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
232 tmp += R[c][setK->data]*J[setK->data][setK->data]*R[d][setK->data];
233 tmp += R[c][setE->data]*J[setE->data][setK->data]*R[d][setK->data];
234 tmp += R[c][setK->data]*J[setK->data][setE->data]*R[d][setE->data];
235 tmp += R[c][setE->data]*J[setE->data][setE->data]*R[d][setE->data];
237 tmpMatrix[c][d] = tmp;
241 /* transp(R,RT,*k,*k); */
242 /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
243 /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
245 // Now {'a','b'} form a dimer
247 append(&setDa,setE->data);
248 append(&setDb,setK->data);
250 // E = E minus {'a', 'b'}
251 setWalker = setE; // 'a'=setE->data
252 while(setWalker != NULL) {
253 if(setWalker->next->data == setK->data) {
254 setWalker2 = setWalker->next;
255 setWalker->next = setWalker->next->next; // delete 'b'
259 setWalker = setWalker->next;
262 setE = setE->next; // delete 'a'
269 complex double W = 0.0 + 0.0*I;
272 if(setSleftover == -1) {
273 W = cexp(I*PI*0.25*(*Q));
275 for(a=0; a<setMCount; a++) {
276 W *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
277 setWalker = setWalker->next;
282 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])));
283 setWalker = setWalker->next;
284 setWalker2 = setWalker2->next;
287 complex double W0 = 0.0 + 0.0*I;
288 complex double W1 = 0.0 + 0.0*I;
290 W0 = cexp(I*PI*0.25*(*Q));
292 for(a=0; a<setMCount; a++) {
293 W0 *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
294 setWalker = setWalker->next;
299 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])));
300 setWalker = setWalker->next;
301 setWalker2 = setWalker2->next;
303 // W1 = cexp(I*PI*0.25*(*Q+D[setS->data]));
304 W1 = cexp(I*PI*0.25*(*Q+D[setSleftover]));
306 for(a=0; a<setMCount; a++) {
307 // W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setS->data])));
308 W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setSleftover])));
309 setWalker = setWalker->next;
314 // 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])));
315 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])));
316 setWalker = setWalker->next;
317 setWalker2 = setWalker2->next;
322 deallocate_mem(&tmpMatrix, *k);
323 deallocate_mem(&R, *k);
324 deallocate_mem(&RT, *k);
326 freeList(&setDa); freeList(&setDb); freeList(&setM); freeList(&setE); freeList(&setS);
333 void append(struct Node** head_ref, int new_data)
335 struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
336 struct Node *last = *head_ref;
338 new_node->data = new_data;
339 new_node->next = NULL;
341 // if the linked list is empty, then make the new node as head
342 if (*head_ref == NULL) {
343 *head_ref = new_node;
347 // else traverse till the last node
348 while (last->next != NULL) {
351 last->next = new_node;
356 /*void freeList(struct Node** head_ref)
359 struct Node *second_last_walker = *head_ref;
360 struct Node *walker = *head_ref;
362 // else traverse till the last node
364 while(walker->next != NULL) {
365 while(walker->next != NULL) {
366 //printf("!%d\n", walker->data);
367 second_last_walker = walker;
368 walker = walker->next;
369 //printf("%d\n", walker->data);
370 //printf("!%d\n", second_last_walker->data);
373 second_last_walker->next = NULL;
375 //printf("!!%d\n", second_last_walker->data);
383 void freeList(struct Node** head_ref)
386 struct Node *walker = *head_ref;
387 struct Node *walker2;
389 // else traverse till the last node
391 walker2 = walker->next;
392 while(walker2 != NULL) {
395 walker2 = walker->next;
403 void printLinkedList(struct Node *node)
410 printf(" %d ", node->data);