exponentialsum.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <complex.h>
   6
   7 #include "matrix.h"
   8 #include "exponentialsum.h"
   9
  10 // linked list
  11 struct Node {
  12   int data;
  13   struct Node* next;
  14 };
  15
  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);
  20
  21 complex double exponentialsum(int *k, int *Q, int *D, int **J) {
  22
  23   int a, b, c, d;
  24
  25   int tmp;
  26
  27   struct Node* setS = NULL;
  28   struct Node* setWalker = NULL;
  29   struct Node* setWalker2 = NULL;
  30
  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)
  33
  34   int setSleftover = -1;  // -1 means nothing left over
  35
  36   int** R;  // basis change matrix R and R^T
  37   int** RT;
  38   R = calloc(*k, sizeof(int*));
  39   RT = calloc(*k, sizeof(int*));
  40   // initially set it to the k-dimensional identity matrix
  41   for(a=0; a<*k; a++) {
  42     R[a] = calloc(*k, sizeof(int));
  43     R[a][a] = 1;
  44     RT[a] = calloc(*k, sizeof(int));
  45     RT[a][a] = 1;
  46   }
  47
  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*));
  51   for(a=0; a<*k; a++)
  52     tmpMatrix[a] = calloc(*k, sizeof(int));
  53
  54   for(a=0; a<*k; a++) {
  55     if(D[a] == 2 || D[a] == 6) {
  56       append(&setS, a);
  57     }
  58   }
  59
  60   // Choose 'a' to be the first element in the linked list setS; a=setS->data
  61
  62   if(setS != NULL) {
  63     setWalker = setS;
  64     while(setWalker->next != NULL) {
  65       setWalker = setWalker->next;
  66       R[setWalker->data][setS->data] += 1;
  67     }
  68
  69     tmpVec = calloc(*k, sizeof(int));
  70     for(c=0;c<*k;c++) {
  71       tmp = 0;
  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; */
  79       if(c!=setS->data)
  80         tmp += R[c][c]*D[c];
  81       tmp += R[c][setS->data]*D[setS->data];
  82       if(setS->data>c)
  83         tmp += J[c][setS->data]*R[c][c]*R[c][setS->data];
  84       else if(setS->data<c)
  85         tmp += J[setS->data][c]*R[c][setS->data]*R[c][c];
  86       tmpVec[c] = tmp%8;
  87     }
  88     memcpy(D, tmpVec, sizeof(int)*(*k));
  89     free(tmpVec);
  90
  91     // Eq. 50 update of J
  92     for(c=0;c<*k;c++) {
  93       for(d=0;d<*k; d++) {
  94         tmp = 0;
  95         if(c!=setS->data) {
  96           if(d!=setS->data) {
  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];
 100           }else {
 101             tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
 102           }
 103         } else {
 104           if(d!=setS->data) {
 105             tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
 106           }
 107         }
 108         tmp += R[c][setS->data]*J[setS->data][setS->data]*R[d][setS->data];
 109         tmp %= 8;
 110         tmpMatrix[c][d] = tmp;
 111       }
 112     }
 113     for(c=0; c<*k; c++)
 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); */
 118
 119     // remove everything but first element from setS
 120     //tmp = setS->data;
 121     //freeList(&setS);
 122     //append(&setS, tmp);
 123     //setS->next = NULL;
 124     setSleftover = setS->data;
 125   }
 126
 127   struct Node* setE = NULL;
 128
 129   for(a=0;a<*k;a++) {
 130     //if((setS != NULL && a != setS->data) || setS == NULL)
 131     if((setSleftover != -1 && a != setSleftover) || setSleftover == -1)
 132       append(&setE, a);
 133   }
 134
 135   struct Node* setM = NULL;
 136   int setMCount = 0; // number of elements in M
 137
 138   int r = 0;
 139
 140   while(setE != NULL) {
 141     // let 'a' be setE->data (the first element in setE)
 142     struct Node* setK = NULL;
 143     setWalker = setE;
 144     while(setWalker->next != NULL) {
 145       setWalker = setWalker->next;
 146       if(J[setE->data][setWalker->data] == 4)
 147         append(&setK, setWalker->data);
 148     }
 149
 150     if(setK == NULL) { // Found new monomer {'a'}
 151       append(&setM, setE->data);
 152       setMCount++;
 153
 154       if(setE->next != NULL) {
 155         setE->data = setE->next->data;
 156         setE->next = setE->next->next;
 157       }
 158       else setE = NULL;
 159     } else {
 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++)
 166         //  R[a][b] = 0;
 167         R[a][a] = 1;
 168       }
 169
 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
 174         }
 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;
 178       }
 179       // update D and J
 180       tmpVec = calloc(*k, sizeof(int));
 181       for(c=0;c<*k;c++) {
 182         tmp = 0;
 183         //tmpVec[c] = 0;
 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]; */
 189         /* } */
 190         if(c!=setK->data && c!=setE->data)
 191           tmp += R[c][c]*D[c];
 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];
 198         if(setK->data>c)
 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];
 202         if(setE->data>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];
 206         tmpVec[c] = tmp%8;
 207       }
 208       memcpy(D, tmpVec, sizeof(int)*(*k));
 209       free(tmpVec);
 210
 211       // Eq. 50 update of J
 212       for(c=0;c<*k;c++) {
 213         for(d=0;d<*k; d++) {
 214           tmp = 0;
 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];
 222             }else {
 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];
 225             }
 226           } else {
 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];
 230             }
 231           }
 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];
 236           tmp %= 8;
 237           tmpMatrix[c][d] = tmp;
 238         }
 239       }
 240       J = tmpMatrix;
 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); */
 244
 245       // Now {'a','b'} form a dimer
 246       r++;
 247       append(&setDa,setE->data);
 248       append(&setDb,setK->data);
 249
 250       // E = E minus {'a', 'b'}
 251       setWalker = setE; // 'a'=setE->data
 252       while(setWalker != NULL) {
 253         if(setWalker->next->data == setK->data) {
 254           free(setWalker->next);
 255           setWalker->next = setWalker->next->next; // delete 'b'
 256           break;
 257         }
 258         setWalker = setWalker->next;
 259       }
 260       setWalker = setE;
 261       setE = setE->next; // delete 'a'
 262       free(setWalker);
 263     }
 264
 265     freeList(&setK);
 266   }
 267
 268   complex double W = 0.0 + 0.0*I;
 269
 270 //if(setS == NULL) {
 271   if(setSleftover == -1) {
 272     W = cexp(I*PI*0.25*(*Q));
 273     setWalker = setM;
 274     for(a=0; a<setMCount; a++) {
 275       W *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
 276       setWalker = setWalker->next;
 277     }
 278     setWalker = setDa;
 279     setWalker2 = setDb;
 280     for(a=0; a<r; a++) {
 281       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])));
 282       setWalker = setWalker->next;
 283       setWalker2 = setWalker2->next;
 284     }
 285   } else {
 286     complex double W0 = 0.0 + 0.0*I;
 287     complex double W1 = 0.0 + 0.0*I;
 288
 289     W0 = cexp(I*PI*0.25*(*Q));
 290     setWalker = setM;
 291     for(a=0; a<setMCount; a++) {
 292       W0 *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
 293       setWalker = setWalker->next;
 294     }
 295     setWalker = setDa;
 296     setWalker2 = setDb;
 297     for(a=0; a<r; a++) {
 298       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])));
 299       setWalker = setWalker->next;
 300       setWalker2 = setWalker2->next;
 301     }
 302     //    W1 = cexp(I*PI*0.25*(*Q+D[setS->data]));
 303     W1 = cexp(I*PI*0.25*(*Q+D[setSleftover]));
 304     setWalker = setM;
 305     for(a=0; a<setMCount; a++) {
 306       //      W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setS->data])));
 307       W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setSleftover])));
 308       setWalker = setWalker->next;
 309     }
 310     setWalker = setDa;
 311     setWalker2 = setDb;
 312     for(a=0; a<r; a++) {
 313       //      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])));
 314       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])));
 315       setWalker = setWalker->next;
 316       setWalker2 = setWalker2->next;
 317     }
 318     W = W0 + W1;
 319   }
 320
 321   deallocate_mem(&tmpMatrix, *k);
 322   deallocate_mem(&R, *k);
 323   deallocate_mem(&RT, *k);
 324
 325   freeList(&setDa); freeList(&setDb); freeList(&setM); freeList(&setE); freeList(&setS);
 326
 327   return(W);
 328
 329 }
 330
 331
 332 void append(struct Node** head_ref, int new_data)
 333 {
 334     struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
 335     struct Node *last = *head_ref;
 336
 337     new_node->data = new_data;
 338     new_node->next = NULL;
 339
 340     // if the linked list is empty, then make the new node as head
 341     if (*head_ref == NULL) {
 342       *head_ref = new_node;
 343       return;
 344     }
 345
 346     // else traverse till the last node
 347     while (last->next != NULL) {
 348         last = last->next;
 349     }
 350     last->next = new_node;
 351
 352     return;
 353 }
 354
 355 /*void freeList(struct Node** head_ref)
 356 {
 357
 358   struct Node *second_last_walker = *head_ref;
 359   struct Node *walker = *head_ref;
 360
 361     // else traverse till the last node
 362     if(walker != NULL) {
 363       while(walker->next != NULL) {
 364         while(walker->next != NULL) {
 365           //printf("!%d\n", walker->data);
 366           second_last_walker = walker;
 367           walker = walker->next;
 368           //printf("%d\n", walker->data);
 369           //printf("!%d\n", second_last_walker->data);
 370         }
 371         free(walker);
 372         second_last_walker->next = NULL;
 373         walker = *head_ref;
 374         //printf("!!%d\n", second_last_walker->data);
 375       }
 376     }
 377     free(walker);
 378
 379     return;
 380     }*/
 381
 382 void freeList(struct Node** head_ref)
 383 {
 384
 385   struct Node *walker = *head_ref;
 386   struct Node *walker2;
 387
 388     // else traverse till the last node
 389     if(walker != NULL) {
 390       walker2 = walker->next;
 391       while(walker2 != NULL) {
 392         free(walker);
 393         walker = walker2;
 394         walker2 = walker->next;
 395       }
 396     }
 397     free(walker);
 398
 399     return;
 400 }
 401
 402 void printLinkedList(struct Node *node)
 403 {
 404   if(node == NULL)
 405     printf("NULL\n");
 406   else {
 407     while (node != NULL)
 408       {
 409         printf(" %d ", node->data);
 410         node = node->next;
 411       }
 412   }
 413   printf("\n");
 414 }
 415