strongsim6.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <complex.h>
   6 #include "matrix.h"
   7 #include "exponentialsum.h"
   8 #include "shrink.h"
   9 #include "extend.h"
  10 #include "measurepauli.h"
  11 #include "innerproduct.h"
  12
  13 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
  14
  15 // order of matrix elements is [row][column]!!!
  16
  17 int main()
  18 {
  19
  20   int N;              // number of qubits
  21   scanf("%d", &N);
  22
  23   if(N%6 != 0) {
  24     printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n");
  25     return 1;
  26   }
  27
  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)
  29   scanf("%d", &T);
  30
  31   int omega;
  32   int alpha[N], beta[N], gamma[N], delta[N];
  33   int Paulicounter = 0;
  34
  35   int i, j, k, l, m;
  36
  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);
  44
  45   // b60
  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} };
  47   // b66
  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} };
  49   // e6
  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}  };
  51   // o6
  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}  };
  53   // k6
  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} };
  55   // phiprime
  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}  };
  57   // phidoubleprime
  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}  };
  59
  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*));
  64
  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));
  66
  67   K = calloc(pow(7,N/6), sizeof(int));
  68
  69   double complex origGamma[(int)pow(7,N/6)];
  70   int *origK, *origQ, **origD, ***origJ;
  71   int ***origG, ***origGBar, **origh;
  72
  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*));
  75
  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));
  77
  78   origK = calloc(pow(7,N/6), sizeof(int));
  79
  80   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
  81
  82
  83   for(j=0; j<pow(7,N/6); j++) { // there will be 7^(N/6) combinations when using k=3 tensor factors
  84
  85     combination = j;
  86
  87     K[j] = 0.0;
  88
  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);
  91       combination /= 7;
  92     }
  93     combination = j;
  94     origK[j] = K[j];
  95
  96     Gamma[j] = 1.0;
  97
  98     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
  99     h[j] = calloc(N, sizeof(int));
 100
 101     if(K[j] > 0) {
 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));
 105     }
 106
 107     origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
 108     origh[j] = calloc(N, sizeof(int));
 109
 110     if(K[j] > 0) {
 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));
 114     }
 115
 116     for(k=0; k<N; k++) {
 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));
 119     }
 120
 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++) {
 124
 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;
 126
 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);
 128
 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) {
 133         case 0:
 134           D[j][Kcounter+l] = D1[l];
 135           break;
 136         case 1:
 137           D[j][Kcounter+l] = D2[l];
 138           break;
 139         case 2:
 140           D[j][Kcounter+l] = D3[l];
 141           break;
 142         case 3:
 143           D[j][Kcounter+l] = D4[l];
 144             break;
 145         case 4:
 146           D[j][Kcounter+l] = D5[l];
 147           break;
 148         case 5:
 149           D[j][Kcounter+l] = D6[l];
 150           break;
 151         case 6:
 152           D[j][Kcounter+l] = D7[l];
 153           break;
 154         default:
 155           printf("error");
 156           return 1;
 157           }
 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) {
 161           case 0:
 162             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
 163             break;
 164           case 1:
 165               J[j][Kcounter+l][Kcounter+m] = J2[l][m];
 166               break;
 167           case 2:
 168             J[j][Kcounter+l][Kcounter+m] = J3[l][m];
 169             break;
 170           case 3:
 171             J[j][Kcounter+l][Kcounter+m] = J4[l][m];
 172             break;
 173             case 4:
 174               J[j][Kcounter+l][Kcounter+m] = J5[l][m];
 175               break;
 176           case 5:
 177             J[j][Kcounter+l][Kcounter+m] = J6[l][m];
 178             break;
 179           case 6:
 180             J[j][Kcounter+l][Kcounter+m] = J7[l][m];
 181             break;
 182           default:
 183             printf("error");
 184             return 1;
 185           }
 186         }
 187       }
 188
 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];
 191       }
 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];
 197         }
 198       }
 199       Kcounter = Kcounter + Kcombo;
 200
 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));
 206
 207       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 208
 209       combination /= 7; // shift to the right by one (in base-7 arithmetic)
 210     }
 211     //printf("!\n");
 212
 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
 214     combination = 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];
 227         }
 228       }
 229       Kcounter = Kcounter + (n1-Kcombo);
 230
 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); */
 235
 236       combination /= 7;
 237     }
 238     for(k=0; k<N; k++) {
 239       memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 240     }
 241     for(k=0; k<K[j]; k++) {
 242       memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 243     }
 244
 245     memcpy(origh[j], h[j], N*sizeof(int));
 246     memcpy(origD[j], D[j], K[j]*sizeof(int));
 247
 248     printf("G[%d]:\n", j);
 249     printMatrix(G[j], N, N);
 250     printf("GBar[%d]:\n", j);
 251     printMatrix(GBar[j], N, N);
 252
 253     printf("h[%d]:\n", j);
 254     printVector(h[j], N);
 255
 256     printf("J[%d]:\n", j);
 257     printMatrix(J[j], K[j], K[j]);
 258
 259     printf("D[%d]:\n", j);
 260     printVector(D[j], K[j]);
 261
 262     printf("Q[%d]=%d\n", j, Q[j]);
 263
 264   }
 265   //exit(0);
 266   memcpy(origGamma, Gamma, pow(7,N/6)*sizeof(double complex));
 267
 268   memcpy(origQ, Q, pow(7,N/6)*sizeof(int));
 269
 270   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
 271
 272     Paulicounter++;
 273     if(Paulicounter > N) {
 274       printf("Error: Number of Paulis is greater than N!\n");
 275       return 1;
 276     }
 277
 278     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 279     // Y_i = -I*Z*X
 280     for(i=0; i<N; i++) {
 281       if(delta[i]){
 282         omega += 3; // -I = I^3
 283         beta[i] = delta[i];
 284         gamma[i] = delta[i];
 285       }
 286     }
 287
 288     printf("*******\n");
 289     printf("*******\n");
 290     printf("omega=%d\n", omega);
 291     printf("X:\n");
 292     printVector(gamma, N);
 293     printf("Z:\n");
 294     printVector(beta, N);
 295     printf("*******\n");
 296     printf("*******\n");
 297
 298     for(j=0; j<pow(7,N/6); j++) { // the kets
 299
 300       printf("========\n");
 301       printf("before:\n");
 302       printf("K=%d\n", K[j]);
 303       printf("h:\n");
 304       printVector(h[j], N);
 305       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 306       printf("G:\n");
 307       printMatrix(G[j], N, N);
 308       printf("GBar:\n");
 309       printMatrix(GBar[j], N, N);
 310       printf("Q=%d\n", Q[j]);
 311       printf("D:\n");
 312       printVector(D[j], K[j]);
 313       printf("J:\n");
 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]);
 318       printf("h:\n");
 319       printVector(h[j], N);
 320       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 321       printf("G:\n");
 322       printMatrix(G[j], N, N);
 323       printf("GBar:\n");
 324       printMatrix(GBar[j], N, N);
 325       printf("Q=%d\n", Q[j]);
 326       printf("D:\n");
 327       printVector(D[j], K[j]);
 328       printf("J:\n");
 329       printMatrix(J[j], K[j], K[j]);
 330
 331     }
 332
 333   }
 334
 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;
 340     }
 341   }
 342
 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)));
 346   else
 347     printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
 348
 349   return 0;
 350
 351 }
 352
 353
 354
 355 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 356 {
 357
 358   int newomega, newalpha, newbeta, newgamma, newdelta;
 359   int i;
 360
 361   if(scanf("%d", &newomega) != EOF) {
 362     *omega = newomega;
 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");
 366         exit(0);
 367       }
 368       if(newalpha+newbeta+newgamma+newdelta > 1) {
 369         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 370         exit(0);
 371       }
 372       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 373     }
 374     return 1;
 375   } else
 376     return 0;
 377
 378 }