deleted incorrect var numPaulis and replaced with numrandomsteps in test2 bash script
[strong_simulation_stabilizer_rank.git] / 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 }