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
23 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)
27 int alpha[N], beta[N], gamma[N], delta[N];
33 int n1 = 1; int k1 = 0; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {0}; int (*(J1[])) = { (int[]) {0} };
34 int n2 = 1; int k2 = 0; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {1}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
36 int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
37 double complex Gamma[(int)pow(2,N)]; // prefactor in front of resultant state
38 G = calloc(pow(2,N),sizeof(int*)); GBar = calloc(pow(2,N),sizeof(int*));
39 h = calloc(pow(2,N),sizeof(int*));
41 J = calloc(pow(2,N),sizeof(int*)); D = calloc(pow(2,N),sizeof(int*)); Q = calloc(pow(2,N),sizeof(int));
43 K = calloc(pow(2,N), sizeof(int));
45 double complex origGamma[(int)pow(2,N)];
46 int *origK, *origQ, **origD, ***origJ;
47 int ***origG, ***origGBar, **origh;
49 origG = calloc(pow(2,N),sizeof(int*)); origGBar = calloc(pow(2,N),sizeof(int*));
50 origh = calloc(pow(2,N),sizeof(int*));
52 origJ = calloc(pow(2,N),sizeof(int*)); origD = calloc(pow(2,N),sizeof(int*)); origQ = calloc(pow(2,N),sizeof(int));
54 origK = calloc(pow(2,N), sizeof(int));
56 int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
59 for(j=0; j<pow(2,N); j++) { // there will be 2^N combinations when using k=1 tensor factors
64 G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
65 h[j] = calloc(N, sizeof(int));
69 J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
72 origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
73 origh[j] = calloc(N, sizeof(int));
76 origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
81 Gamma[j] *= (1.0/sqrt(2.0))*((combination & 0x1)*cexp(0.25*PI*I) + (~combination & 0x1)*cexp(0.0*PI*I));
83 G[j][k] = calloc(N, sizeof(int*)); GBar[j][k] = calloc(N, sizeof(int*));
86 J[j] = calloc(K[j], sizeof(int*));
88 origG[j][k] = calloc(N, sizeof(int*)); origGBar[j][k] = calloc(N, sizeof(int*));
91 origJ[j] = calloc(K[j], sizeof(int*));
93 h[j][k] = (combination & 0x1)*h2[0] + (~combination & 0x1)*h1[0];
94 for(l=0; l<n1; l++) { // assuming n1=n2
95 G[j][k][k*n1+l] = (combination & 0x1)*G2[0][l] + (~combination & 0x1)*G1[0][l];
96 GBar[j][k][k*n1+l] = (combination & 0x1)*GBar2[0][l] + (~combination & 0x1)*GBar1[0][l];
98 memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
100 combination /= 2; // shift to the right by one (in base-2 arithmetic)
103 memcpy(origh[j], h[j], N*sizeof(int));
107 memcpy(origGamma, Gamma, pow(2,N)*sizeof(double complex));
109 while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
112 if(Paulicounter > N) {
113 printf("Error: Number of Paulis is greater than N!\n");
117 // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
121 omega += 3; // -I = I^3
127 for(j=0; j<pow(2,N); j++) { // the kets
129 Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
135 double complex amplitude = 0.0 + 0.0*I;
136 for(i=0; i<pow(2,N); i++) { // the bras
137 for(j=0; j<pow(2,N); j++) {
138 // check to see if second arguments are modified!!! They shouldn't be!
139 double complex newamplitude = conj(origGamma[i])*Gamma[j]*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]);
140 amplitude = amplitude + newamplitude;
144 printf("amplitude:\n");
145 if(creal(amplitude+0.00000001)>0)
146 printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
148 printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
156 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
159 int newomega, newalpha, newbeta, newgamma, newdelta;
162 if(scanf("%d", &newomega) != EOF) {
164 for(i=0; i<numqubits; i++) {
165 if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
166 printf("Error: Too few input coeffs!\n");
169 if(newalpha+newbeta+newgamma+newdelta > 1) {
170 printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
173 alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;