innerproduct.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 "shrink.h"
   9 #include "extend.h"
  10 #include "exponentialsum.h"
  11 #include "innerproduct.h"
  12
  13
  14 /****************************************************************************
  15  * Note: Arguments are copied and not modified!
  16  ****************************************************************************/
  17
  18 double complex innerproduct(int n1, int k1, int *h1, int **G1, int **GBar1, int Q1, int *D1, int **J1, int n2, int k2, int *h2, int **G2, int **GBar2, int Q2, int *D2, int **J2) {
  19   // n1 should equal n2 (we assume this)
  20
  21   int a, b, c;
  22
  23   int n = n1;
  24   int k = k1;
  25   int *h = calloc(n, sizeof(int));
  26   int **G = calloc(n, sizeof(int*));
  27   int **GBar = calloc(n, sizeof(int*));
  28   for(a=0; a<n; a++) {
  29     h[a] = h1[a];
  30     G[a] = calloc(n, sizeof(int));
  31     GBar[a] = calloc(n, sizeof(int));
  32     for(b=0; b<n; b++) {
  33       G[a][b] = G1[a][b];
  34       GBar[a][b] = GBar1[a][b];
  35     }
  36   }
  37
  38   int Q1copy = Q1;
  39   int *D1copy = calloc(k1, sizeof(int));
  40   int **J1copy = calloc(k1, sizeof(int*));
  41   for(a=0; a<k1; a++) {
  42     D1copy[a] = D1[a];
  43     J1copy[a] = calloc(k1, sizeof(int));
  44     for(b=0; b<k1; b++) {
  45       J1copy[a][b] = J1[a][b];
  46     }
  47   }
  48
  49   int Q2copy = Q2;
  50   int *D2copy = calloc(k2, sizeof(int));
  51   int **J2copy = calloc(k2, sizeof(int*));
  52   for(a=0; a<k2; a++) {
  53     D2copy[a] = D2[a];
  54     J2copy[a] = calloc(k2, sizeof(int));
  55     for(b=0; b<k2; b++) {
  56       J2copy[a][b] = J2[a][b];
  57     }
  58   }
  59
  60   int alpha;
  61   int **R = calloc(k, sizeof(int*));
  62   for(a=0; a<k; a++)
  63     R[a] = calloc(k2, sizeof(int));
  64   int **RT = calloc(k2, sizeof(int*));
  65   for(a=0; a<k2; a++)
  66     RT[a] = calloc(k, sizeof(int));
  67   char shrinkResult;
  68
  69   for(a=k2; a<n2; a++) {
  70     alpha = dotProductMod(h2, GBar2[a], n2, 2);
  71
  72     shrinkResult = shrink(n, &k, h, G, GBar, &Q1copy, &D1copy, &J1copy, GBar2[a], alpha);
  73     if(shrinkResult == 'E') {// if EMPTY
  74       free(h);
  75       deallocate_mem(&G, n); deallocate_mem(&GBar, n);
  76       free(D1copy);
  77       deallocate_mem(&J1copy, k);
  78       free(D2copy);
  79       deallocate_mem(&J2copy, k2);
  80       deallocate_mem(&R, k);
  81       deallocate_mem(&RT, k2);
  82       return(0.0+0.0*I);
  83     }
  84   }
  85
  86   // Now K = K1 intersect K2 = (n1, k1, h1, G1, GBar1)
  87   int *y = calloc(k2, sizeof(int));
  88   int *tmpVec = calloc(n2, sizeof(int));
  89   addVectorMod(h,h2,tmpVec,n2,2);
  90   for(a=0; a<k2; a++) {
  91     y[a] = dotProductMod(tmpVec,GBar2[a],n2,2);
  92     for(b=0; b<k; b++) {
  93       R[b][a] = dotProductMod(G[b],GBar2[a],n2,2);
  94     }
  95   }
  96
  97   // IF YOU UNCOMMENT THIS THEN YOU WILL BE CHANGING h2!!!
  98   /* for(a=0; a<k2; a++) { */
  99   /*   for(b=0; b<n2; b++) { */
 100   /*     h2[b] = (h2[b] + y[a]*G2[a][b])%2; */
 101   /*   } */
 102   /* } */
 103   /* printf("h2 (should now match h):\n"); */
 104   /* printVector(h2, n2); */
 105   /* printf("h:\n"); */
 106   /* printVector(h, n2); */
 107
 108   free(tmpVec);
 109
 110   // update Q2 and D2 using Eqs. 52 and 53 with y
 111   for(a=0; a<k2; a++) {
 112     Q2copy = Q2copy + D2copy[a]*y[a];
 113     for(b=a+1; b<k2; b++) {
 114       Q2copy = Q2copy + J2copy[a][b]*y[a]*y[b];
 115     }
 116   }
 117   Q2copy = Q2copy%8;
 118
 119   for(a=0; a<k2; a++) {
 120     for(b=0; b<k2; b++)
 121       D2copy[a] = D2copy[a] + J2copy[a][b]*y[b];
 122     D2copy[a] = D2copy[a]%8;
 123   }
 124
 125   free(y);
 126
 127   // update D2copy and J2copy using Eqs. 49 and 50 with R
 128   // (note that R is a kxk2 matrix so this may reduce the dimension of D2copy and J2copy)
 129   y = calloc(k, sizeof(int));
 130   for(c=0;c<k;c++) {
 131     y[c] = 0;
 132     for(a=0;a<k2;a++) {
 133       y[c] = y[c] + R[c][a]*D2copy[a];
 134       for(b=a+1;b<k2;b++)
 135         y[c] = y[c] + J2copy[a][b]*R[c][a]*R[c][b];
 136     }
 137     y[c] = y[c]%8;
 138   }
 139   free(D2copy);
 140   D2copy = calloc(k, sizeof(int));
 141   memcpy(D2copy, y, sizeof(int)*k);
 142   free(y);
 143
 144   transp(R,RT,k,k2);
 145   int **tmpMatrix = calloc(k2, sizeof(int*));
 146   for(a=0; a<k2; a++)
 147     tmpMatrix[a] = calloc(k, sizeof(int));
 148   multMatrixMod(J2copy,RT,tmpMatrix,k2,k2,k2,k,8);
 149   deallocate_mem(&J2copy, k2);
 150   J2copy = calloc(k, sizeof(int*));
 151   for(a=0;a<k;a++)
 152     J2copy[a] = calloc(k, sizeof(int));
 153   multMatrixMod(R,tmpMatrix,J2copy,k,k2,k2,k,8);
 154
 155   deallocate_mem(&tmpMatrix, k2);
 156
 157   deallocate_mem(&R, k);
 158   deallocate_mem(&RT, k2);
 159
 160   // Now (Q1, D1, J1) and (Q2, D2, J2) are defined in the same basis
 161
 162   int Q = (Q1copy - Q2copy)%8 < 0 ? (Q1copy - Q2copy)%8 + 8 : (Q1copy - Q2copy)%8;
 163   int *D = calloc(k, sizeof(int));
 164   int **J = calloc(k, sizeof(int*));
 165
 166   for(a=0; a<k; a++) { // again, if k=0 then no need to fill these
 167     J[a] = calloc(k, sizeof(int));
 168     D[a] = (D1copy[a] - D2copy[a])%8 < 0 ? (D1copy[a] - D2copy[a])%8 + 8 : (D1copy[a] - D2copy[a])%8;
 169     for(b=0; b<k; b++)
 170       J[a][b] = (J1copy[a][b] - J2copy[a][b])%8 < 0 ? (J1copy[a][b] - J2copy[a][b])%8 + 8 : (J1copy[a][b] - J2copy[a][b])%8;
 171
 172   }
 173
 174   free(h);
 175   deallocate_mem(&G, n); deallocate_mem(&GBar, n);
 176
 177   double complex amplitude = pow(2,-0.5*(k1+k2))*exponentialsum(&k, &Q, D, J);
 178
 179   free(D);
 180   deallocate_mem(&J, k);
 181
 182   free(D1copy);
 183   deallocate_mem(&J1copy, k);
 184
 185   free(D2copy);
 186   deallocate_mem(&J2copy, k);
 187
 188   return(amplitude);
 189
 190 }