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