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           setWalker2 = setWalker->next;
 255           setWalker->next = setWalker->next->next; // delete 'b'
 256           free(setWalker2);
 257           break;
 258         }
 259         setWalker = setWalker->next;
 260       }
 261       setWalker = setE;
 262       setE = setE->next; // delete 'a'
 263       free(setWalker);
 264     }
 265
 266     freeList(&setK);
 267   }
 268
 269   complex double W = 0.0 + 0.0*I;
 270
 271 //if(setS == NULL) {
 272   if(setSleftover == -1) {
 273     W = cexp(I*PI*0.25*(*Q));
 274     setWalker = setM;
 275     for(a=0; a<setMCount; a++) {
 276       W *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
 277       setWalker = setWalker->next;
 278     }
 279     setWalker = setDa;
 280     setWalker2 = setDb;
 281     for(a=0; a<r; a++) {
 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;
 285     }
 286   } else {
 287     complex double W0 = 0.0 + 0.0*I;
 288     complex double W1 = 0.0 + 0.0*I;
 289
 290     W0 = cexp(I*PI*0.25*(*Q));
 291     setWalker = setM;
 292     for(a=0; a<setMCount; a++) {
 293       W0 *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
 294       setWalker = setWalker->next;
 295     }
 296     setWalker = setDa;
 297     setWalker2 = setDb;
 298     for(a=0; a<r; a++) {
 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;
 302     }
 303     //    W1 = cexp(I*PI*0.25*(*Q+D[setS->data]));
 304     W1 = cexp(I*PI*0.25*(*Q+D[setSleftover]));
 305     setWalker = setM;
 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;
 310     }
 311     setWalker = setDa;
 312     setWalker2 = setDb;
 313     for(a=0; a<r; a++) {
 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;
 318     }
 319     W = W0 + W1;
 320   }
 321
 322   deallocate_mem(&tmpMatrix, *k);
 323   deallocate_mem(&R, *k);
 324   deallocate_mem(&RT, *k);
 325
 326   freeList(&setDa); freeList(&setDb); freeList(&setM); freeList(&setE); freeList(&setS);
 327
 328   return(W);
 329
 330 }
 331
 332
 333 void append(struct Node** head_ref, int new_data)
 334 {
 335     struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
 336     struct Node *last = *head_ref;
 337
 338     new_node->data = new_data;
 339     new_node->next = NULL;
 340
 341     // if the linked list is empty, then make the new node as head
 342     if (*head_ref == NULL) {
 343       *head_ref = new_node;
 344       return;
 345     }
 346
 347     // else traverse till the last node
 348     while (last->next != NULL) {
 349         last = last->next;
 350     }
 351     last->next = new_node;
 352
 353     return;
 354 }
 355
 356 /*void freeList(struct Node** head_ref)
 357 {
 358
 359   struct Node *second_last_walker = *head_ref;
 360   struct Node *walker = *head_ref;
 361
 362     // else traverse till the last node
 363     if(walker != NULL) {
 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);
 371         }
 372         free(walker);
 373         second_last_walker->next = NULL;
 374         walker = *head_ref;
 375         //printf("!!%d\n", second_last_walker->data);
 376       }
 377     }
 378     free(walker);
 379
 380     return;
 381     }*/
 382
 383 void freeList(struct Node** head_ref)
 384 {
 385
 386   struct Node *walker = *head_ref;
 387   struct Node *walker2;
 388
 389     // else traverse till the last node
 390     if(walker != NULL) {
 391       walker2 = walker->next;
 392       while(walker2 != NULL) {
 393         free(walker);
 394         walker = walker2;
 395         walker2 = walker->next;
 396       }
 397     }
 398     free(walker);
 399
 400     return;
 401 }
 402
 403 void printLinkedList(struct Node *node)
 404 {
 405   if(node == NULL)
 406     printf("NULL\n");
 407   else {
 408     while (node != NULL)
 409       {
 410         printf(" %d ", node->data);
 411         node = node->next;
 412       }
 413   }
 414   printf("\n");
 415 }
 416