strongsim1_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("strongsim2_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%2 != 0) {
  35   //  printf("'N' needs to be a multiple of 2 for a k=2 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=0; i<PdN; i++) {
  71     tmp = 0.0;
  72     for(j=0; j<=(i+1); 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 coeffa = (1.0/sqrt(2.0))*cexp(0.0*PI*I);
  87   double complex coeffb = (1.0/sqrt(2.0))*cexp(0.25*PI*I);
  88
  89   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} };
  90   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} };
  91
  92   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
  93   double complex Gamma[(int)pow(2,N)]; // prefactor in front of resultant state
  94   G = calloc(pow(2,N),sizeof(int*)); GBar = calloc(pow(2,N),sizeof(int*));
  95   h = calloc(pow(2,N),sizeof(int*));
  96
  97   J = calloc(pow(2,N),sizeof(int*)); D = calloc(pow(2,N),sizeof(int*)); Q = calloc(pow(2,N),sizeof(int));
  98
  99   K = calloc(pow(2,N), sizeof(int));
 100
 101   int origK, origQ, *origD;
 102   int **origJ;
 103   int **origG, **origGBar;
 104   int *origh;
 105   double complex origGamma;
 106
 107   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
 108
 109
 110   for(j=0; j<pow(2,N); j++) { // there will be 2^(N) combinations when using k=12 tensor factors
 111
 112     combination = j;
 113
 114     K[j] = 0.0;
 115
 116     for(k=0; k<N; k++) {
 117       K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 118       combination /= 2;
 119     }
 120     combination = j;
 121
 122     Gamma[j] = 1.0;
 123
 124     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
 125     h[j] = calloc(N, sizeof(int));
 126
 127     if(K[j] > 0) {
 128       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
 129       for(k=0; k<K[j]; k++)
 130         J[j][k] = calloc(K[j], sizeof(int));
 131     }
 132
 133     for(k=0; k<N; k++) {
 134       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
 135     }
 136
 137     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 2) of 'j' in that we go through with 'k'
 138     int Kcombo; // Kcombo stores the k<(n1=n2) dimension of the member of the combination that we are currently adding
 139     for(k=0; k<N; k++) {
 140
 141       Q[j] += (((combination%2)==1)*Q2 + ((combination%2)==0)*Q1);
 142
 143
 144       Gamma[j] *= (((combination%2)==1)*coeffb + ((combination%2)==0)*coeffa);
 145
 146       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 147       for(l=0; l<Kcombo; l++) {
 148         // D1 has a different number of rows 'l' than D2 and D3 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
 149         switch(combination%2) {
 150         case 0:
 151           D[j][Kcounter+l] = D1[l];
 152           break;
 153         case 1:
 154           D[j][Kcounter+l] = D2[l];
 155           break;
 156         default:
 157           printf("error");
 158           return 1;
 159           }
 160         for(m=0; m<Kcombo; m++) {
 161           // J1 has a different number of rows 'l' than J2 and J3 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
 162           switch(combination%2) {
 163           case 0:
 164             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
 165             break;
 166           case 1:
 167             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
 168             break;
 169           default:
 170             printf("error");
 171             return 1;
 172           }
 173         }
 174       }
 175
 176       for(l=0; l<n1; l++) { // assuming n1=n2=n3
 177         h[j][k*n1+l] = (((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l]);
 178       }
 179       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
 180       for(l=0; l<Kcombo; l++) {
 181         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 182           G[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
 183           GBar[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
 184         }
 185       }
 186       Kcounter = Kcounter + Kcombo;
 187
 188       /* printf("intermediate G[%d]:\n", j); */
 189       /* printMatrix(G[j], N, N); */
 190       /* printf("intermediate GBar[%d]:\n", j); */
 191       /* printMatrix(GBar[j], N, N); */
 192       //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 193
 194       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 195
 196       combination /= 2; // shift to the right by one (in base-7 arithmetic)
 197     }
 198     //printf("!\n");
 199
 200     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
 201     combination = j;
 202     for(k=0; k<(N); k++) {
 203       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 204       //printf("Kcounter=%d\n", Kcounter);
 205       // G and GBar rows that are outside the first 'k' spanning basis states
 206       for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
 207         //printf("l=%d\n", l);
 208         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 209           /* printf("m=%d\n", m); */
 210           /* printf("Kcounter+l=%d\n", Kcounter+l); */
 211           /* printf("k*n1+m=%d\n", k*n1+m); */
 212           G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
 213           GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
 214         }
 215       }
 216       Kcounter = Kcounter + (n1-Kcombo);
 217
 218       /* printf("intermediate G[%d]:\n", j); */
 219       /* printMatrix(G[j], N, N); */
 220       /* printf("intermediate GBar[%d]:\n", j); */
 221       /* printMatrix(GBar[j], N, N); */
 222
 223       combination /= 2;
 224     }
 225
 226     /*printf("G[%d]:\n", j);
 227     printMatrix(G[j], N, N);
 228     printf("GBar[%d]:\n", j);
 229     printMatrix(GBar[j], N, N);
 230
 231     printf("h[%d]:\n", j);
 232     printVector(h[j], N);
 233
 234     printf("J[%d]:\n", j);
 235     printMatrix(J[j], K[j], K[j]);
 236
 237     printf("D[%d]:\n", j);
 238     printVector(D[j], K[j]);
 239
 240     printf("Q[%d]=%d\n", j, Q[j]);*/
 241
 242   }
 243   //exit(0);
 244
 245   while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
 246
 247     if((Paulicounter+1) > N) {
 248       printf("Error: Number of Paulis is greater than N!\n");
 249       return 1;
 250     }
 251
 252     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 253     // Y_i = -I*Z*X
 254     for(i=0; i<N; i++) {
 255       if(delta[Paulicounter][i]){
 256         omega[Paulicounter] += 3; // -I = I^3
 257         beta[Paulicounter][i] = delta[Paulicounter][i];
 258         gamma[Paulicounter][i] = delta[Paulicounter][i];
 259       }
 260     }
 261
 262     /*printf("*******\n");
 263     printf("*******\n");
 264     printf("omega=%d\n", omega);
 265     printf("X:\n");
 266     printVector(gamma, N);
 267     printf("Z:\n");
 268     printVector(beta, N);
 269     printf("*******\n");
 270     printf("*******\n");*/
 271
 272     //for(j=0; j<pow(2,N); j++) { // the kets
 273
 274       /*printf("========\n");
 275       printf("before:\n");
 276       printf("K=%d\n", K[j]);
 277       printf("h:\n");
 278       printVector(h[j], N);
 279       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 280       printf("G:\n");
 281       printMatrix(G[j], N, N);
 282       printf("GBar:\n");
 283       printMatrix(GBar[j], N, N);
 284       printf("Q=%d\n", Q[j]);
 285       printf("D:\n");
 286       printVector(D[j], K[j]);
 287       printf("J:\n");
 288       printMatrix(J[j], K[j], K[j]);*/
 289       //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega[Paulicounter], gamma[Paulicounter], beta[Paulicounter]);
 290       /*printf("\nafter:\n");
 291       printf("K=%d\n", K[j]);
 292       printf("h:\n");
 293       printVector(h[j], N);
 294       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
 295       printf("G:\n");
 296       printMatrix(G[j], N, N);
 297       printf("GBar:\n");
 298       printMatrix(GBar[j], N, N);
 299       printf("Q=%d\n", Q[j]);
 300       printf("D:\n");
 301       printVector(D[j], K[j]);
 302       printf("J:\n");
 303       printMatrix(J[j], K[j], K[j]);*/
 304
 305     //}
 306
 307     Paulicounter++;
 308   }
 309
 310   double complex amplitude = 0.0 + 0.0*I;
 311   for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
 312     //printf("i=%d\n", i);
 313
 314     randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
 315     /*printf("origK=%d\n", origK);
 316     printf("origG:\n");
 317     printMatrix(origG, N, N);
 318     printf("origGBar:\n");
 319     printMatrix(origGBar, N, N);
 320     printf("origh:\n");
 321     printVector(origh, N);
 322     printf("origQ=%d\n", origQ);
 323     printf("origD:\n");
 324     printVector(origD, origK);
 325     printf("origJ:\n");
 326     printMatrix(origJ, origK, origK);*/
 327
 328     origGamma = 1.0 + 0.0*I;
 329
 330     for(k=Paulicounter-1; k>=0; k--) { // go backwards through the measurements since acting on the bra (doesn't matter if they commute though)
 331       origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
 332       //printf("k=%d\n", k);
 333     }
 334     /*printf("origK=%d\n", origK);
 335     printf("origG:\n");
 336     printMatrix(origG, N, N);
 337     printf("origGBar:\n");
 338     printMatrix(origGBar, N, N);
 339     printf("origh:\n");
 340     printVector(origh, N);*/
 341
 342     double complex stabstateaverage = 0.0 + 0.0*I;
 343
 344     for(j=0; j<pow(2,N); j++) {
 345       //printf("j=%d\n", j);
 346       /*  printf("K=%d\n", K[j]);
 347     printf("G:\n");
 348     printMatrix(G[j], N, N);
 349     printf("h:\n");
 350     printVector(h[j], N);
 351     printf("Q=%d\n", Q[j]);
 352     printf("D:\n");
 353     printVector(D[j], K[j]);
 354     printf("J:\n");
 355     printMatrix(J[j], K[j], K[j]);*/
 356       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);
 357       /*printf("newamplitude = %lf + %lf i\n", creal(newamplitude), cimag(newamplitude));
 358       printf("origGamma = %lf + %lf i\n", creal(origGamma), cimag(origGamma));
 359       printf("Gamma = %lf + %lf i\n", creal(Gamma[j]), cimag(Gamma[j]));*/
 360       stabstateaverage = stabstateaverage + conj(origGamma)*Gamma[j]*newamplitude;
 361     }
 362     amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double complex)(NUMSTABSTATESAMPLES))*cpow(2.0,T);
 363
 364     deallocate_mem(&origG, N);
 365     deallocate_mem(&origGBar, N);
 366     free(origh);
 367     deallocate_mem(&origJ, origK);
 368     free(origD);
 369   }
 370
 371   printf("amplitude:\n");
 372   if(creal(amplitude+ZEROTHRESHOLD)>0)
 373     printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
 374   else
 375     printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
 376
 377
 378
 379   for(i=0; i<PdN; i++)
 380     free(Pd[i]);
 381   free(Pd);
 382
 383   return 0;
 384 }
 385
 386 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 387 {
 388
 389   int newomega, newalpha, newbeta, newgamma, newdelta;
 390   int i;
 391
 392   if(scanf("%d", &newomega) != EOF) {
 393     *omega = newomega;
 394     for(i=0; i<numqubits; i++) {
 395       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
 396         printf("Error: Too few input coeffs!\n");
 397         exit(0);
 398       }
 399       if(newalpha+newbeta+newgamma+newdelta > 1) {
 400         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 401         exit(0);
 402       }
 403       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 404     }
 405     return 1;
 406   } else
 407     return 0;
 408
 409 }