final v1.0 files
[weak_simulation_stab_extent.git] / weaksim.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include <complex.h>
6 #include <time.h>
7 #include "matrix.h"
8 #include "exponentialsum.h"
9 #include "shrink.h"
10 #include "extend.h"
11 #include "measurepauli.h"
12 #include "innerproduct.h"
13 #include "sparsify.h"
14 #include "omp.h"
15
16 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
17
18 // order of matrix elements is [row][column]!!!
19
20 int main(int argc, char *argv[])
21 {
22
23   if(argc != 5 && argc != 6) {
24     printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\" \"seed\" \"(optional) sample number prefactor\"\n");
25     exit(0);
26   }
27
28   double additiveErrorDelta = atof(argv[1]);        // additive error delta
29   double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
30   int coherentSampling = atoi(argv[3]);           // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
31   int seed = atoi(argv[4]); // random number seed
32   double samplePrefactor = 1.0;
33   if(argc == 6)
34     samplePrefactor = atof(argv[5]); // factor that multiplies k, the number of samples used
35
36   int N;              // number of qubits
37   scanf("%d", &N);
38
39   int T;              // number of T gate magic states (set to the first 'K' of the 'N' qubits -- the rest are set to the '0' computational basis state)
40   scanf("%d", &T);  
41
42   printf("phi = %lf\n", phi);
43
44   int omega;
45   int alpha[N], beta[N], gamma[N], delta[N];
46   int Paulicounter = 0;
47
48   int i, j, k, l, m;
49
50   //|T> = e^(pi*i/8)/2*sqrt(4-2*sqrt(2))* (sqrt(-i)*(|0>+i|1>)/sqrt(2) + (|0>+|1>)/sqrt(2))
51
52   //double additiveErrorDelta = 0.1;
53
54   double complex coeffa = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*cexp(-PI*I*0.25)*I/sqrt(2.0)*(-I+cexp(-0.25*PI*I))*(-I+cexp(I*phi)))); // factor of cexp(PI*I/8.0)*cexp(-PI*I*0.25) comes from converting (|0>+|1>)/sqrt(2) under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
55   double complex coeffb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*phi)))); // factor of cexp(PI*I/8.0) comes from converting |0> under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
56   // alternative coefficient to use instead of coeffb to get overall entangled state
57   double complex coeffbb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*0.25*PI))));
58
59   int n1 = 1; int k1 = 1; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
60   int n2 = 1; int k2 = 1; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {0}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
61
62   long* stabStateIndices;
63   int numStabStates;
64
65   //srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
66   srand((unsigned)seed); // seeding the random number generator for sparsify()
67
68   if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling, samplePrefactor))
69     return 1;
70
71   /* printf("checking: numStabStateIndices:\n"); */
72   /* for(i=0; i<numStabStates; i++) */
73   /*   printf("%ld ", stabStateIndices[i]); */
74   /* printf("\n"); */
75   /* fflush(stdout); */
76
77   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
78   double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
79   G = calloc(numStabStates,sizeof(int*)); GBar = calloc(numStabStates,sizeof(int*));
80   h = calloc(numStabStates,sizeof(int*));
81   
82   J = calloc(numStabStates,sizeof(int*)); D = calloc(numStabStates,sizeof(int*)); Q = calloc(numStabStates,sizeof(int));
83
84   K = calloc(numStabStates, sizeof(int));
85
86   double complex origGamma[(int)numStabStates];
87   int *origK, *origQ, **origD, ***origJ;
88   int ***origG, ***origGBar, **origh;
89
90   origG = calloc(numStabStates,sizeof(int*)); origGBar = calloc(numStabStates,sizeof(int*));
91   origh = calloc(numStabStates,sizeof(int*));
92   
93   origJ = calloc(numStabStates,sizeof(int*)); origD = calloc(numStabStates,sizeof(int*)); origQ = calloc(numStabStates,sizeof(int));
94
95   origK = calloc(numStabStates, sizeof(int));
96
97   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
98   
99   double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
100
101   for(j=0; j<numStabStates; j++) {
102
103     combination = stabStateIndices[j];
104
105     K[j] = 0.0;
106     
107     for(k=0; k<N; k++) {
108       K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
109       combination /= 2;
110     }
111     combination = stabStateIndices[j];
112     origK[j] = K[j];
113
114     Gamma[j] = 1.0;
115     Gamma[j] *= L1Norm/((double)numStabStates);
116
117     // the coefficients which are a product of 'coeffa', 'coeffb', 'coeffbb' (that are subsequently multiplied into Gamma[j]) is multiplied by 'norm'
118     //Gamma[j] *= norm;
119
120     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
121     h[j] = calloc(N, sizeof(int));
122
123     if(K[j] > 0) {
124       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
125       for(k=0; k<K[j]; k++)
126         J[j][k] = calloc(K[j], sizeof(int));
127     }
128
129     origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
130     origh[j] = calloc(N, sizeof(int));
131
132     if(K[j] > 0) {
133       origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
134       for(k=0; k<K[j]; k++)
135         origJ[j][k] = calloc(K[j], sizeof(int));
136     }
137
138     for(k=0; k<N; k++) {
139       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
140       origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
141     }
142
143     int Kcounter = 0; // Kcounter keeps track of the K<=N that we have added already to the G rows etc. for each combination that is indexed by the digits (base 3) of 'j' in that we go through with 'k'
144     int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
145
146     // if combination contains at least one instance of the second state, i.e. contains the 0 digit in binary, then we want to have it have one instance of coeffb instead of coeffbb
147     for(k=0; k<N; k++) {
148       if(combination%2==1) {
149         Gamma[j] *= coeffb/coeffbb;
150         break; // break out of loop
151       }
152       combination /= 2; // shift to the right by one (in base-2 arithmetic)
153     }
154     combination = stabStateIndices[j];
155
156     for(k=0; k<N; k++) {
157
158       Q[j] += ((combination%2)==1)*Q2 + ((combination%2)==0)*Q1;
159       
160       Gamma[j] *= (((combination%2)==1)*coeffbb + ((combination%2)==0)*coeffa); // only assign coeffbb instead of coeffb; coeffb replaces one instance of coeffbb before this loop
161
162       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
163       for(l=0; l<Kcombo; l++) {
164           // D1 may have a different number of rows 'l' than D2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
165           switch(combination%2) {
166           case 0:
167             D[j][Kcounter+l] = D1[l];
168             break;
169           case 1:
170             D[j][Kcounter+l] = D2[l];
171             break;
172           default:
173             printf("error");
174             return 1;
175           }
176         for(m=0; m<Kcombo; m++) {
177           // J1 may have a different number of rows 'l' than J2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
178           switch(combination%2) {
179           case 0:
180             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
181             break;
182           case 1:
183             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
184             break;
185           default:
186             printf("error");
187             return 1;
188           }
189         }
190       }
191
192       for(l=0; l<n1; l++) { // assuming n1=n2
193         h[j][k*n1+l] = ((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l];
194       }
195       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
196       for(l=0; l<Kcombo; l++) {
197         for(m=0; m<n1; m++) { // assuming n1=n2
198           G[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
199           GBar[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
200         }
201       }
202       Kcounter = Kcounter + Kcombo;
203       
204       /* printf("intermediate G[%d]:\n", j); */
205       /* printMatrix(G[j], N, N); */
206       /* printf("intermediate GBar[%d]:\n", j); */
207       /* printMatrix(GBar[j], N, N); */
208       //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
209
210       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
211       
212       combination /= 2; // shift to the right by one (in base-2 arithmetic)
213     }
214     //printf("!\n");
215
216     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
217     combination = stabStateIndices[j];
218     for(k=0; k<(N); k++) {
219       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
220       //printf("Kcounter=%d\n", Kcounter);
221       // G and GBar rows that are outside the first 'k' spanning basis states
222       for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
223         //printf("l=%d\n", l);
224         for(m=0; m<n1; m++) { // assuming n1=n2=n3
225           /* printf("m=%d\n", m); */
226           /* printf("Kcounter+l=%d\n", Kcounter+l); */
227           /* printf("k*n1+m=%d\n", k*n1+m); */
228           G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
229           GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
230         }
231       }
232       Kcounter = Kcounter + (n1-Kcombo);
233
234       /* printf("intermediate G[%d]:\n", j); */
235       /* printMatrix(G[j], N, N); */
236       /* printf("intermediate GBar[%d]:\n", j); */
237       /* printMatrix(GBar[j], N, N); */
238       
239       combination /= 2;
240     }
241     for(k=0; k<N; k++) {
242       memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
243     }
244     for(k=0; k<K[j]; k++) {
245       memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));      
246     }
247
248     memcpy(origh[j], h[j], N*sizeof(int));
249     memcpy(origD[j], D[j], K[j]*sizeof(int));
250
251   }
252   //exit(0);
253   memcpy(origGamma, Gamma, numStabStates*sizeof(double complex));
254
255   memcpy(origQ, Q, numStabStates*sizeof(int));
256
257   double complex amplitude = 0.0 + 0.0*I;
258   double complex lastamplitude = 0.0 + 0.0*I;
259
260   double complex probability = 1.0 + 0.0*I;
261
262   // the first measurement should be the identity
263   omega = 0;
264   for(i=0; i<N; i++) {
265     alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
266   }
267
268   for(j=0; j<numStabStates; j++) { // the kets
269     Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
270   }
271
272   double complex newamplitude;
273   
274   #pragma omp parallel
275   {
276   #pragma omp for schedule(static) collapse(2)
277   for(i=0; i<numStabStates; i++) { // the bras
278     for(j=0; j<numStabStates; j++) {
279       newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
280       amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
281     }
282   }
283   }
284
285   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
286   
287     Paulicounter++;
288     if(Paulicounter > N) {
289       printf("Error: Number of Paulis is greater than N!\n");
290       return 1;
291     }
292     
293     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
294     // Y_i = -I*Z*X
295     for(i=0; i<N; i++) {
296       if(delta[i]){
297         omega += 3; // -I = I^3
298         beta[i] = delta[i];
299         gamma[i] = delta[i];
300       }
301     }
302
303
304     for(j=0; j<numStabStates; j++) { // the kets
305
306       Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
307
308     }
309
310     lastamplitude = amplitude;
311     amplitude = 0.0 + 0.0*I;
312     #pragma omp parallel
313     {
314     #pragma omp for schedule(static) collapse(2)
315     for(i=0; i<numStabStates; i++) { // the bras
316       for(j=0; j<numStabStates; j++) {
317         newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
318         amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
319       }
320     }
321     }
322     printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
323     printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
324
325     probability *= amplitude/lastamplitude;
326     //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
327
328     printf("probability = %lf\n", cabs(probability));
329     printf("!\n");
330   }
331
332   //printf("numStabStates=%d\n", numStabStates);
333   printf("L1Norm=%lf\n", L1Norm);
334
335   printf("\namplitude:\n");
336   if(creal(amplitude+0.00000001)>0)
337     printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
338   else
339     printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
340   //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
341   printf("\nabs(amplitude):\n");
342   printf("%lf\n", cabs(amplitude));
343
344   printf("probability = %lf\n", cabs(probability));
345
346   return 0;
347
348 }
349
350
351
352 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
353 {
354     
355   int newomega, newalpha, newbeta, newgamma, newdelta;
356   int i;
357
358   if(scanf("%d", &newomega) != EOF) {
359     *omega = newomega;
360     for(i=0; i<numqubits; i++) {
361       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
362         printf("Error: Too few input coeffs!\n");
363         exit(0);
364       }
365       if(newalpha+newbeta+newgamma+newdelta > 1) {
366         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
367         exit(0);
368       }
369       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
370     }
371     return 1;
372   } else
373     return 0;
374     
375 }