7 #include "exponentialsum.h"
10 #include "measurepauli.h"
11 #include "innerproduct.h"
13 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
15 // order of matrix elements is [row][column]!!!
20 int N; // number of qubits
24 printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n");
28 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)
32 int alpha[N], beta[N], gamma[N], delta[N];
37 double complex coeffb60 = (-16.0+12.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
38 double complex coeffb66 = (96.0-68.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
39 double complex coeffe6 = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)*cpow(2.0,2.5);
40 double complex coeffo6 = (-14.0+10.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-14.0*PI*I/8.0)*cpow(2.0,2.5);
41 double complex coeffk6 = (7.0-5.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-8.0*PI*I/8.0)*4.0*csqrt(2.0)*cpow(2.0,0.5);
42 double complex coeffphiprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
43 double complex coeffphidprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
46 int n1 = 6; int k1 = 6; int (*(G1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h1[] = {0, 0, 0, 0, 0, 0}; int Q1 = 0; int D1[] = {0, 0, 0, 0, 0, 0}; int (*(J1[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
48 int n2 = 6; int k2 = 6; int (*(G2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h2[] = {0, 0, 0, 0, 0, 0}; int Q2 = 4; int D2[] = {4, 4, 4, 4, 4, 4}; int (*(J2[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
50 int n3 = 6; int k3 = 5; int (*(G3[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar3[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h3[] = {1, 0, 0, 0, 0, 0}; int Q3 = 4; int D3[] = {0, 0, 0, 0, 0}; int (*(J3[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
52 int n4 = 6; int k4 = 5; int (*(G4[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar4[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h4[] = {0, 0, 0, 0, 0, 0}; int Q4 = 4; int D4[] = {4, 4, 4, 4, 4}; int (*(J4[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
54 int n5 = 6; int k5 = 1; int (*(G5[])) = { (int[]) {1, 1, 1, 1, 1, 1}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1} }; int (*(GBar5[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1} }; int h5[] = {1, 1, 1, 1, 1, 1}; int Q5 = 6; int D5[] = {2}; int (*(J5[])) = { (int[]) {4} };
56 int n6 = 6; int k6 = 5; int (*(G6[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar6[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h6[] = {1, 0, 0, 0, 0, 0}; int Q6 = 0; int D6[] = {0, 0, 0, 0, 0}; int (*(J6[])) = { (int[]) {0, 4, 0, 0, 4}, (int[]) {4, 0, 4, 0, 0}, (int[]) {0, 4, 0, 4, 0}, (int[]) {0, 0, 4, 0, 4}, (int[]) {4, 0, 0, 4, 0} };
58 int n7 = 6; int k7 = 5; int (*(G7[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar7[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h7[] = {1, 0, 0, 0, 0, 0}; int Q7 = 0; int D7[] = {0, 0, 0, 0, 0}; int (*(J7[])) = { (int[]) {0, 0, 4, 4, 0}, (int[]) {0, 0, 0, 4, 4}, (int[]) {4, 0, 0, 0, 4}, (int[]) {4, 4, 0, 0, 0}, (int[]) {0, 4, 4, 0, 0} };
60 int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
61 double complex Gamma[(int)pow(7,N/6)]; // prefactor in front of resultant state
62 G = calloc(pow(7,N/6),sizeof(int*)); GBar = calloc(pow(7,N/6),sizeof(int*));
63 h = calloc(pow(7,N/6),sizeof(int*));
65 J = calloc(pow(7,N/6),sizeof(int*)); D = calloc(pow(7,N/6),sizeof(int*)); Q = calloc(pow(7,N/6),sizeof(int));
67 K = calloc(pow(7,N/6), sizeof(int));
69 double complex origGamma[(int)pow(7,N/6)];
70 int *origK, *origQ, **origD, ***origJ;
71 int ***origG, ***origGBar, **origh;
73 origG = calloc(pow(7,N/6),sizeof(int*)); origGBar = calloc(pow(7,N/6),sizeof(int*));
74 origh = calloc(pow(7,N/6),sizeof(int*));
76 origJ = calloc(pow(7,N/6),sizeof(int*)); origD = calloc(pow(7,N/6),sizeof(int*)); origQ = calloc(pow(7,N/6),sizeof(int));
78 origK = calloc(pow(7,N/6), sizeof(int));
80 int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
83 for(j=0; j<pow(7,N/6); j++) { // there will be 7^(N/6) combinations when using k=3 tensor factors
89 for(k=0; k<N/6; k++) {
90 K[j] += (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
98 G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
99 h[j] = calloc(N, sizeof(int));
102 J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
103 for(k=0; k<K[j]; k++)
104 J[j][k] = calloc(K[j], sizeof(int));
107 origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
108 origh[j] = calloc(N, sizeof(int));
111 origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
112 for(k=0; k<K[j]; k++)
113 origJ[j][k] = calloc(K[j], sizeof(int));
117 G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
118 origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
121 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'
122 int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
123 for(k=0; k<N/6; k++) {
125 Q[j] += ((combination%7)==6)*Q7 + ((combination%7)==5)*Q6 + ((combination%7)==4)*Q5 + ((combination%7)==3)*Q4 + ((combination%7)==2)*Q3 + ((combination%7)==1)*Q2 + ((combination%7)==0)*Q1;
127 Gamma[j] *= (((combination%7)==6)*coeffphidprime + ((combination%7)==5)*coeffphiprime + ((combination%7)==4)*coeffk6 + ((combination%7)==3)*coeffo6 + ((combination%7)==2)*coeffe6 + ((combination%7)==1)*coeffb66 + ((combination%7)==0)*coeffb60);
129 Kcombo = (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
130 for(l=0; l<Kcombo; l++) {
131 // D1 has a different number of rows 'l' than D2 and D3 so you need to use something like 'switch' to check combination%3 without going out of bound of J1
132 switch(combination%7) {
134 D[j][Kcounter+l] = D1[l];
137 D[j][Kcounter+l] = D2[l];
140 D[j][Kcounter+l] = D3[l];
143 D[j][Kcounter+l] = D4[l];
146 D[j][Kcounter+l] = D5[l];
149 D[j][Kcounter+l] = D6[l];
152 D[j][Kcounter+l] = D7[l];
158 for(m=0; m<Kcombo; m++) {
159 // J1 has a different number of rows 'l' than J2 and J3 so you need to use something like 'switch' to check combination%3 without going out of bound of J1
160 switch(combination%7) {
162 J[j][Kcounter+l][Kcounter+m] = J1[l][m];
165 J[j][Kcounter+l][Kcounter+m] = J2[l][m];
168 J[j][Kcounter+l][Kcounter+m] = J3[l][m];
171 J[j][Kcounter+l][Kcounter+m] = J4[l][m];
174 J[j][Kcounter+l][Kcounter+m] = J5[l][m];
177 J[j][Kcounter+l][Kcounter+m] = J6[l][m];
180 J[j][Kcounter+l][Kcounter+m] = J7[l][m];
189 for(l=0; l<n1; l++) { // assuming n1=n2=n3
190 h[j][k*n1+l] = ((combination%7)==6)*h7[l] + ((combination%7)==5)*h6[l] + ((combination%7)==4)*h5[l] + ((combination%7)==3)*h4[l] + ((combination%7)==2)*h3[l] + ((combination%7)==1)*h2[l] + ((combination%7)==0)*h1[l];
192 // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
193 for(l=0; l<Kcombo; l++) {
194 for(m=0; m<n1; m++) { // assuming n1=n2=n3
195 G[j][Kcounter+l][k*n1+m] = ((combination%7)==6)*G7[l][m] + ((combination%7)==5)*G6[l][m] + ((combination%7)==4)*G5[l][m] + ((combination%7)==3)*G4[l][m] + ((combination%7)==2)*G3[l][m] + ((combination%7)==1)*G2[l][m] + ((combination%7)==0)*G1[l][m];
196 GBar[j][Kcounter+l][k*n1+m] = ((combination%7)==6)*GBar7[l][m] + ((combination%7)==5)*GBar6[l][m] + ((combination%7)==4)*GBar5[l][m] + ((combination%7)==3)*GBar4[l][m] + ((combination%7)==2)*GBar3[l][m] + ((combination%7)==1)*GBar2[l][m] + ((combination%7)==0)*GBar1[l][m];
199 Kcounter = Kcounter + Kcombo;
201 /* printf("intermediate G[%d]:\n", j); */
202 /* printMatrix(G[j], N, N); */
203 /* printf("intermediate GBar[%d]:\n", j); */
204 /* printMatrix(GBar[j], N, N); */
205 //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
207 //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
209 combination /= 7; // shift to the right by one (in base-7 arithmetic)
213 // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
215 for(k=0; k<(N/6); k++) {
216 Kcombo = (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
217 //printf("Kcounter=%d\n", Kcounter);
218 // G and GBar rows that are outside the first 'k' spanning basis states
219 for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
220 //printf("l=%d\n", l);
221 for(m=0; m<n1; m++) { // assuming n1=n2=n3
222 /* printf("m=%d\n", m); */
223 /* printf("Kcounter+l=%d\n", Kcounter+l); */
224 /* printf("k*n1+m=%d\n", k*n1+m); */
225 G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%7)==6)*G7[l][m] + ((combination%7)==5)*G6[l][m] + ((combination%7)==4)*G5[l][m] + ((combination%7)==3)*G4[l][m] + ((combination%7)==2)*G3[l][m] + ((combination%7)==1)*G2[l][m] + ((combination%7)==0)*G1[l][m];
226 GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%7)==6)*GBar7[l][m] + ((combination%7)==5)*GBar6[l][m] + ((combination%7)==4)*GBar5[l][m] + ((combination%7)==3)*GBar4[l][m] + ((combination%7)==2)*GBar3[l][m] + ((combination%7)==1)*GBar2[l][m] + ((combination%7)==0)*GBar1[l][m];
229 Kcounter = Kcounter + (n1-Kcombo);
231 /* printf("intermediate G[%d]:\n", j); */
232 /* printMatrix(G[j], N, N); */
233 /* printf("intermediate GBar[%d]:\n", j); */
234 /* printMatrix(GBar[j], N, N); */
239 memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
241 for(k=0; k<K[j]; k++) {
242 memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
245 memcpy(origh[j], h[j], N*sizeof(int));
246 memcpy(origD[j], D[j], K[j]*sizeof(int));
248 printf("G[%d]:\n", j);
249 printMatrix(G[j], N, N);
250 printf("GBar[%d]:\n", j);
251 printMatrix(GBar[j], N, N);
253 printf("h[%d]:\n", j);
254 printVector(h[j], N);
256 printf("J[%d]:\n", j);
257 printMatrix(J[j], K[j], K[j]);
259 printf("D[%d]:\n", j);
260 printVector(D[j], K[j]);
262 printf("Q[%d]=%d\n", j, Q[j]);
266 memcpy(origGamma, Gamma, pow(7,N/6)*sizeof(double complex));
268 memcpy(origQ, Q, pow(7,N/6)*sizeof(int));
270 while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
273 if(Paulicounter > N) {
274 printf("Error: Number of Paulis is greater than N!\n");
278 // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
282 omega += 3; // -I = I^3
290 printf("omega=%d\n", omega);
292 printVector(gamma, N);
294 printVector(beta, N);
298 for(j=0; j<pow(7,N/6); j++) { // the kets
300 printf("========\n");
302 printf("K=%d\n", K[j]);
304 printVector(h[j], N);
305 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
307 printMatrix(G[j], N, N);
309 printMatrix(GBar[j], N, N);
310 printf("Q=%d\n", Q[j]);
312 printVector(D[j], K[j]);
314 printMatrix(J[j], K[j], K[j]);
315 Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
316 printf("\nafter:\n");
317 printf("K=%d\n", K[j]);
319 printVector(h[j], N);
320 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
322 printMatrix(G[j], N, N);
324 printMatrix(GBar[j], N, N);
325 printf("Q=%d\n", Q[j]);
327 printVector(D[j], K[j]);
329 printMatrix(J[j], K[j], K[j]);
335 double complex amplitude = 0.0 + 0.0*I;
336 for(i=0; i<pow(7,N/6); i++) { // the bras
337 for(j=0; j<pow(7,N/6); j++) {
338 double complex 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]);
339 amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
343 printf("amplitude:\n");
344 if(creal(amplitude+0.00000001)>0)
345 printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
347 printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
355 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
358 int newomega, newalpha, newbeta, newgamma, newdelta;
361 if(scanf("%d", &newomega) != EOF) {
363 for(i=0; i<numqubits; i++) {
364 if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
365 printf("Error: Too few input coeffs!\n");
368 if(newalpha+newbeta+newgamma+newdelta > 1) {
369 printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
372 alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;