weaksim.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 "sparsify.h"
  14
  15 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
  16
  17 // order of matrix elements is [row][column]!!!
  18
  19 int main(int argc, char *argv[])
  20 {
  21
  22   if(argc != 4) {
  23     printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
  24     exit(0);
  25   }
  26
  27   double additiveErrorDelta = atof(argv[1]);        // additive error delta
  28   double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
  29   int coherentSampling = atoi(argv[3]);           // perform coherent sampling (true=1 or false=0)
  30
  31   int N;              // number of qubits
  32   scanf("%d", &N);
  33
  34   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)
  35   scanf("%d", &T);
  36
  37   printf("phi = %lf\n", phi);
  38
  39   int omega;
  40   int alpha[N], beta[N], gamma[N], delta[N];
  41   int Paulicounter = 0;
  42
  43   int i, j, k, l, m;
  44
  45   //|T> = e^(pi*i/8)/2*sqrt(4-2*sqrt(2))* (sqrt(-i)*(|0>+i|1>)/sqrt(2) + (|0>+|1>)/sqrt(2))
  46
  47   //double additiveErrorDelta = 0.1;
  48
  49   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>
  50   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>
  51   // alternative coefficient to use instead of coeffb to get overall entangled state
  52   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))));
  53
  54   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} };
  55   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} };
  56
  57   long* stabStateIndices;
  58   int numStabStates;
  59
  60   srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
  61
  62   if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
  63     return 1;
  64
  65   //printf("checking: numStabStateIndices:\n");
  66   //for(i=0; i<numStabStates; i++)
  67   //  printf("%d ", stabStateIndices[i]);
  68   //printf("\n");
  69   //fflush(stdout);
  70
  71   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
  72   double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
  73   G = calloc(numStabStates,sizeof(int*)); GBar = calloc(numStabStates,sizeof(int*));
  74   h = calloc(numStabStates,sizeof(int*));
  75
  76   J = calloc(numStabStates,sizeof(int*)); D = calloc(numStabStates,sizeof(int*)); Q = calloc(numStabStates,sizeof(int));
  77
  78   K = calloc(numStabStates, sizeof(int));
  79
  80   double complex origGamma[(int)numStabStates];
  81   int *origK, *origQ, **origD, ***origJ;
  82   int ***origG, ***origGBar, **origh;
  83
  84   origG = calloc(numStabStates,sizeof(int*)); origGBar = calloc(numStabStates,sizeof(int*));
  85   origh = calloc(numStabStates,sizeof(int*));
  86
  87   origJ = calloc(numStabStates,sizeof(int*)); origD = calloc(numStabStates,sizeof(int*)); origQ = calloc(numStabStates,sizeof(int));
  88
  89   origK = calloc(numStabStates, sizeof(int));
  90
  91   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
  92
  93   double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
  94
  95   for(j=0; j<numStabStates; j++) {
  96
  97     combination = stabStateIndices[j];
  98
  99     K[j] = 0.0;
 100
 101     for(k=0; k<N; k++) {
 102       K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 103       combination /= 2;
 104     }
 105     combination = stabStateIndices[j];
 106     origK[j] = K[j];
 107
 108     Gamma[j] = 1.0;
 109     Gamma[j] *= L1Norm/((double)numStabStates);
 110
 111     // the coefficients which are a product of 'coeffa', 'coeffb', 'coeffbb' (that are subsequently multiplied into Gamma[j]) is multiplied by 'norm'
 112     //Gamma[j] *= norm;
 113
 114     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
 115     h[j] = calloc(N, sizeof(int));
 116
 117     if(K[j] > 0) {
 118       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
 119       for(k=0; k<K[j]; k++)
 120         J[j][k] = calloc(K[j], sizeof(int));
 121     }
 122
 123     origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
 124     origh[j] = calloc(N, sizeof(int));
 125
 126     if(K[j] > 0) {
 127       origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
 128       for(k=0; k<K[j]; k++)
 129         origJ[j][k] = calloc(K[j], sizeof(int));
 130     }
 131
 132     for(k=0; k<N; k++) {
 133       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
 134       origG[j][k] = calloc(N, sizeof(int)); origGBar[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 3) of 'j' in that we go through with 'k'
 138     int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
 139
 140     // 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
 141     for(k=0; k<N; k++) {
 142       if(combination%2==1) {
 143         Gamma[j] *= coeffb/coeffbb;
 144         break; // break out of loop
 145       }
 146       combination /= 2; // shift to the right by one (in base-2 arithmetic)
 147     }
 148     combination = stabStateIndices[j];
 149
 150     for(k=0; k<N; k++) {
 151
 152       Q[j] += ((combination%2)==1)*Q2 + ((combination%2)==0)*Q1;
 153
 154       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
 155
 156       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 157       for(l=0; l<Kcombo; l++) {
 158           // 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
 159           switch(combination%2) {
 160           case 0:
 161             D[j][Kcounter+l] = D1[l];
 162             break;
 163           case 1:
 164             D[j][Kcounter+l] = D2[l];
 165             break;
 166           default:
 167             printf("error");
 168             return 1;
 169           }
 170         for(m=0; m<Kcombo; m++) {
 171           // 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
 172           switch(combination%2) {
 173           case 0:
 174             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
 175             break;
 176           case 1:
 177             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
 178             break;
 179           default:
 180             printf("error");
 181             return 1;
 182           }
 183         }
 184       }
 185
 186       for(l=0; l<n1; l++) { // assuming n1=n2
 187         h[j][k*n1+l] = ((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l];
 188       }
 189       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
 190       for(l=0; l<Kcombo; l++) {
 191         for(m=0; m<n1; m++) { // assuming n1=n2
 192           G[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
 193           GBar[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
 194         }
 195       }
 196       Kcounter = Kcounter + Kcombo;
 197
 198       /* printf("intermediate G[%d]:\n", j); */
 199       /* printMatrix(G[j], N, N); */
 200       /* printf("intermediate GBar[%d]:\n", j); */
 201       /* printMatrix(GBar[j], N, N); */
 202       //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 203
 204       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 205
 206       combination /= 2; // shift to the right by one (in base-2 arithmetic)
 207     }
 208     //printf("!\n");
 209
 210     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
 211     combination = stabStateIndices[j];
 212     for(k=0; k<(N); k++) {
 213       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
 214       //printf("Kcounter=%d\n", Kcounter);
 215       // G and GBar rows that are outside the first 'k' spanning basis states
 216       for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
 217         //printf("l=%d\n", l);
 218         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 219           /* printf("m=%d\n", m); */
 220           /* printf("Kcounter+l=%d\n", Kcounter+l); */
 221           /* printf("k*n1+m=%d\n", k*n1+m); */
 222           G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
 223           GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
 224         }
 225       }
 226       Kcounter = Kcounter + (n1-Kcombo);
 227
 228       /* printf("intermediate G[%d]:\n", j); */
 229       /* printMatrix(G[j], N, N); */
 230       /* printf("intermediate GBar[%d]:\n", j); */
 231       /* printMatrix(GBar[j], N, N); */
 232
 233       combination /= 2;
 234     }
 235     for(k=0; k<N; k++) {
 236       memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 237     }
 238     for(k=0; k<K[j]; k++) {
 239       memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 240     }
 241
 242     memcpy(origh[j], h[j], N*sizeof(int));
 243     memcpy(origD[j], D[j], K[j]*sizeof(int));
 244
 245   }
 246   //exit(0);
 247   memcpy(origGamma, Gamma, numStabStates*sizeof(double complex));
 248
 249   memcpy(origQ, Q, numStabStates*sizeof(int));
 250
 251   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
 252
 253     Paulicounter++;
 254     if(Paulicounter > N) {
 255       printf("Error: Number of Paulis is greater than N!\n");
 256       return 1;
 257     }
 258
 259     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 260     // Y_i = -I*Z*X
 261     for(i=0; i<N; i++) {
 262       if(delta[i]){
 263         omega += 3; // -I = I^3
 264         beta[i] = delta[i];
 265         gamma[i] = delta[i];
 266       }
 267     }
 268
 269
 270     for(j=0; j<numStabStates; j++) { // the kets
 271
 272       Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
 273
 274     }
 275
 276   }
 277
 278   double complex amplitude = 0.0 + 0.0*I;
 279   for(i=0; i<numStabStates; i++) { // the bras
 280     for(j=0; j<numStabStates; j++) {
 281       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]);
 282       amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
 283     }
 284   }
 285
 286   //printf("numStabStates=%d\n", numStabStates);
 287   printf("L1Norm=%lf\n", L1Norm);
 288
 289   printf("\namplitude:\n");
 290   if(creal(amplitude+0.00000001)>0)
 291     printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
 292   else
 293     printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
 294   //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
 295   printf("\nabs(amplitude):\n");
 296   printf("%lf\n", cabs(amplitude));
 297
 298   return 0;
 299
 300 }
 301
 302
 303
 304 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 305 {
 306
 307   int newomega, newalpha, newbeta, newgamma, newdelta;
 308   int i;
 309
 310   if(scanf("%d", &newomega) != EOF) {
 311     *omega = newomega;
 312     for(i=0; i<numqubits; i++) {
 313       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
 314         printf("Error: Too few input coeffs!\n");
 315         exit(0);
 316       }
 317       if(newalpha+newbeta+newgamma+newdelta > 1) {
 318         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 319         exit(0);
 320       }
 321       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 322     }
 323     return 1;
 324   } else
 325     return 0;
 326
 327 }