strongsim6_relerr.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 "shrinkstar.h"
  10 #include "extend.h"
  11 #include "measurepauli.h"
  12 #include "innerproduct.h"
  13 #include "randomstabilizerstate.h"
  14
  15 #define ZEROTHRESHOLD (0.00000001)
  16
  17 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
  18
  19 // order of matrix elements is [row][column]!!!
  20
  21 int main(int argc, char *argv[])
  22 {
  23
  24   if(argc != 2) {
  25     printf("strongsim6_rellerr argument: \"number of stabilizer state samples\"\n");
  26     exit(0);
  27   }
  28
  29   int NUMSTABSTATESAMPLES = atoi(argv[1]);        // number of stabilizer state samples
  30
  31   int N;                  // number of qubits
  32   scanf("%d", &N);
  33
  34   if(N%6 != 0) {
  35     printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n");
  36     return 1;
  37   }
  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   int omega[N]; // max of N measurements
  43   int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
  44   int Paulicounter = 0;
  45
  46   int i, j, k, l, m;
  47
  48   FILE *fp;
  49   char buff[255];
  50
  51   srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
  52
  53   fp = fopen("Pd.txt", "r");
  54
  55   if(fscanf(fp, "%s", buff) == EOF) {
  56     printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
  57     return 1;
  58   }
  59
  60   double** Pd;
  61
  62   int PdN = atoi(buff);
  63
  64   Pd = calloc(PdN, sizeof(double*));
  65   for(i=0; i<PdN; i++)
  66     Pd[i] = calloc(PdN+1, sizeof(double));
  67
  68   double tmp;
  69
  70   for(i=1; i<PdN; i++) {
  71     tmp = 0.0;
  72     for(j=0; j<=i; j++) {
  73       if(fscanf(fp, "%s", buff) == EOF) {
  74         printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
  75         return 1;
  76       }
  77       Pd[i][j] = atof(buff);
  78       //printf("%e ", Pd[i][j]);
  79       tmp += Pd[i][j];
  80     }
  81     //printf("\n");
  82     //printf("total=%f\n", tmp);
  83   }
  84
  85
  86   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);
  87   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);
  88   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);
  89   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);
  90   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);
  91   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);
  92   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);
  93
  94   // b60
  95   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} };
  96   // b66
  97   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} };
  98   // e6
  99   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}  };
 100   // o6
 101   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}  };
 102   // k6
 103   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} };
 104   // phiprime
 105   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}  };
 106   // phidoubleprime
 107   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}  };
 108
 109
 110   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
 111   double complex Gamma[(int)pow(7,N/6)]; // prefactor in front of resultant state
 112   G = calloc(pow(7,N/6),sizeof(int*)); GBar = calloc(pow(7,N/6),sizeof(int*));
 113   h = calloc(pow(7,N/6),sizeof(int*));
 114
 115   J = calloc(pow(7,N/6),sizeof(int*)); D = calloc(pow(7,N/6),sizeof(int*)); Q = calloc(pow(7,N/6),sizeof(int));
 116
 117   K = calloc(pow(7,N/6), sizeof(int));
 118
 119   int origK, origQ, *origD;
 120   int **origJ;
 121   int **origG, **origGBar;
 122   int *origh;
 123   double complex origGamma;
 124
 125   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
 126
 127
 128   for(j=0; j<pow(7,N/6); j++) { // there will be 7^(N/6) combinations when using k=12 tensor factors
 129
 130     combination = j;
 131
 132     K[j] = 0.0;
 133
 134     for(k=0; k<N/6; k++) {
 135       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);
 136       combination /= 7;
 137     }
 138     combination = j;
 139
 140     Gamma[j] = 1.0;
 141
 142     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
 143     h[j] = calloc(N, sizeof(int));
 144
 145     if(K[j] > 0) {
 146       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
 147       for(k=0; k<K[j]; k++)
 148         J[j][k] = calloc(K[j], sizeof(int));
 149     }
 150
 151     for(k=0; k<N; k++) {
 152       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
 153     }
 154
 155     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'
 156     int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
 157     for(k=0; k<N/6; k++) {
 158
 159       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);
 160
 161
 162       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);
 163
 164       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);
 165       for(l=0; l<Kcombo; l++) {
 166         // 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
 167         switch(combination%7) {
 168         case 0:
 169           D[j][Kcounter+l] = D1[l];
 170           break;
 171         case 1:
 172           D[j][Kcounter+l] = D2[l];
 173           break;
 174         case 2:
 175           D[j][Kcounter+l] = D3[l];
 176           break;
 177         case 3:
 178           D[j][Kcounter+l] = D4[l];
 179             break;
 180         case 4:
 181           D[j][Kcounter+l] = D5[l];
 182           break;
 183         case 5:
 184           D[j][Kcounter+l] = D6[l];
 185           break;
 186         case 6:
 187           D[j][Kcounter+l] = D7[l];
 188           break;
 189         default:
 190           printf("error");
 191           return 1;
 192           }
 193         for(m=0; m<Kcombo; m++) {
 194           // 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
 195           switch(combination%7) {
 196           case 0:
 197             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
 198             break;
 199           case 1:
 200             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
 201             break;
 202           case 2:
 203             J[j][Kcounter+l][Kcounter+m] = J3[l][m];
 204             break;
 205           case 3:
 206             J[j][Kcounter+l][Kcounter+m] = J4[l][m];
 207             break;
 208           case 4:
 209             J[j][Kcounter+l][Kcounter+m] = J5[l][m];
 210             break;
 211           case 5:
 212             J[j][Kcounter+l][Kcounter+m] = J6[l][m];
 213             break;
 214           case 6:
 215             J[j][Kcounter+l][Kcounter+m] = J7[l][m];
 216             break;
 217           default:
 218             printf("error");
 219             return 1;
 220           }
 221         }
 222       }
 223
 224       for(l=0; l<n1; l++) { // assuming n1=n2=n3
 225         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]);
 226       }
 227       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
 228       for(l=0; l<Kcombo; l++) {
 229         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 230           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]);
 231           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]);
 232         }
 233       }
 234       Kcounter = Kcounter + Kcombo;
 235
 236       /* printf("intermediate G[%d]:\n", j); */
 237       /* printMatrix(G[j], N, N); */
 238       /* printf("intermediate GBar[%d]:\n", j); */
 239       /* printMatrix(GBar[j], N, N); */
 240       //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 241
 242       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 243
 244       combination /= 7; // shift to the right by one (in base-7 arithmetic)
 245     }
 246     //printf("!\n");
 247
 248     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
 249     combination = j;
 250     for(k=0; k<(N/6); k++) {
 251       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);
 252       //printf("Kcounter=%d\n", Kcounter);
 253       // G and GBar rows that are outside the first 'k' spanning basis states
 254       for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
 255         //printf("l=%d\n", l);
 256         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 257           /* printf("m=%d\n", m); */
 258           /* printf("Kcounter+l=%d\n", Kcounter+l); */
 259           /* printf("k*n1+m=%d\n", k*n1+m); */
 260           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]);
 261           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]);
 262         }
 263       }
 264       Kcounter = Kcounter + (n1-Kcombo);
 265
 266       /* printf("intermediate G[%d]:\n", j); */
 267       /* printMatrix(G[j], N, N); */
 268       /* printf("intermediate GBar[%d]:\n", j); */
 269       /* printMatrix(GBar[j], N, N); */
 270
 271       combination /= 7;
 272     }
 273
 274     /*printf("G[%d]:\n", j);
 275     printMatrix(G[j], N, N);
 276     printf("GBar[%d]:\n", j);
 277     printMatrix(GBar[j], N, N);
 278
 279     printf("h[%d]:\n", j);
 280     printVector(h[j], N);
 281
 282     printf("J[%d]:\n", j);
 283     printMatrix(J[j], K[j], K[j]);
 284
 285     printf("D[%d]:\n", j);
 286     printVector(D[j], K[j]);
 287
 288     printf("Q[%d]=%d\n", j, Q[j]);*/
 289
 290   }
 291   //exit(0);
 292
 293   while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
 294
 295     if((Paulicounter+1) > N) {
 296       printf("Error: Number of Paulis is greater than N!\n");
 297       return 1;
 298     }
 299
 300     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 301     // Y_i = -I*Z*X
 302     for(i=0; i<N; i++) {
 303       if(delta[Paulicounter][i]){
 304         omega[Paulicounter] += 3; // -I = I^3
 305         beta[Paulicounter][i] = delta[Paulicounter][i];
 306         gamma[Paulicounter][i] = delta[Paulicounter][i];
 307       }
 308     }
 309
 310     /*printf("*******\n");
 311     printf("*******\n");
 312     printf("omega=%d\n", omega);
 313     printf("X:\n");
 314     printVector(gamma, N);
 315     printf("Z:\n");
 316     printVector(beta, N);
 317     printf("*******\n");
 318     printf("*******\n");*/
 319
 320     //for(j=0; j<pow(7,N/6); j++) { // the kets
 321
 322       /*printf("========\n");
 323       printf("before:\n");
 324       printf("K=%d\n", K[j]);
 325       printf("h:\n");
 326       printVector(h[j], N);
 327       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 328       printf("G:\n");
 329       printMatrix(G[j], N, N);
 330       printf("GBar:\n");
 331       printMatrix(GBar[j], N, N);
 332       printf("Q=%d\n", Q[j]);
 333       printf("D:\n");
 334       printVector(D[j], K[j]);
 335       printf("J:\n");
 336       printMatrix(J[j], K[j], K[j]);*/
 337       //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
 338       /*printf("\nafter:\n");
 339       printf("K=%d\n", K[j]);
 340       printf("h:\n");
 341       printVector(h[j], N);
 342       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 343       printf("G:\n");
 344       printMatrix(G[j], N, N);
 345       printf("GBar:\n");
 346       printMatrix(GBar[j], N, N);
 347       printf("Q=%d\n", Q[j]);
 348       printf("D:\n");
 349       printVector(D[j], K[j]);
 350       printf("J:\n");
 351       printMatrix(J[j], K[j], K[j]);*/
 352
 353     //}
 354
 355     Paulicounter++;
 356   }
 357
 358   double complex amplitude = 0.0 + 0.0*I;
 359   for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
 360     //printf("i=%d\n", i);
 361
 362     randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
 363
 364     origGamma = 1.0 + 0.0*I;
 365
 366     for(k=0; k<Paulicounter; k++) {
 367       origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
 368       //printf("k=%d\n", k);
 369   }
 370     /*printf("origK=%d\n", origK);
 371     printf("origG:\n");
 372     printMatrix(origG, N, N);
 373     printf("origGBar:\n");
 374     printMatrix(origGBar, N, N);
 375     printf("origh:\n");
 376     printVector(origh, N);*/
 377
 378     double complex stabstateaverage = 0.0 + 0.0*I;
 379
 380     for(j=0; j<pow(7,N/6); j++) {
 381       //printf("j=%d\n", j);
 382       double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK, origh, origG, origGBar, origQ, origD, origJ);
 383       stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
 384     }
 385     amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
 386
 387     deallocate_mem(&origG, N);
 388     deallocate_mem(&origGBar, N);
 389     free(origh);
 390     deallocate_mem(&origJ, origK);
 391     free(origD);
 392   }
 393
 394   printf("amplitude:\n");
 395   if(creal(amplitude+ZEROTHRESHOLD)>0)
 396     printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
 397   else
 398     printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
 399
 400
 401
 402   for(i=0; i<PdN; i++)
 403     free(Pd[i]);
 404   free(Pd);
 405
 406   return 0;
 407 }
 408
 409 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 410 {
 411
 412   int newomega, newalpha, newbeta, newgamma, newdelta;
 413   int i;
 414
 415   if(scanf("%d", &newomega) != EOF) {
 416     *omega = newomega;
 417     for(i=0; i<numqubits; i++) {
 418       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
 419         printf("Error: Too few input coeffs!\n");
 420         exit(0);
 421       }
 422       if(newalpha+newbeta+newgamma+newdelta > 1) {
 423         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 424         exit(0);
 425       }
 426       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 427     }
 428     return 1;
 429   } else
 430     return 0;
 431
 432 }