From: Lucas <lkocia@s1040071.ca.sandia.gov>
Date: Wed, 2 Feb 2022 00:36:46 +0000 (-0800)
Subject: final v1.0 files
X-Git-Url: https://s3miclassical.com/gitweb/?a=commitdiff_plain;h=HEAD;p=weak_simulation_stab_extent.git

final v1.0 files
---

diff --git a/weaksim.c b/weaksim.c
index 8ec6222..b2c93f4 100644
--- a/weaksim.c
+++ b/weaksim.c
@@ -11,6 +11,7 @@
 #include "measurepauli.h"
 #include "innerproduct.h"
 #include "sparsify.h"
+#include "omp.h"
 
 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
 
@@ -19,14 +20,18 @@ int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, i
 int main(int argc, char *argv[])
 {
 
-  if(argc != 4) {
-    printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+  if(argc != 5 && argc != 6) {
+    printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\" \"seed\" \"(optional) sample number prefactor\"\n");
     exit(0);
   }
 
   double additiveErrorDelta = atof(argv[1]);        // additive error delta
   double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
-  int coherentSampling = atoi(argv[3]);           // perform coherent sampling (true=1 or false=0)
+  int coherentSampling = atoi(argv[3]);           // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
+  int seed = atoi(argv[4]); // random number seed
+  double samplePrefactor = 1.0;
+  if(argc == 6)
+    samplePrefactor = atof(argv[5]); // factor that multiplies k, the number of samples used
 
   int N;              // number of qubits
   scanf("%d", &N);
@@ -57,16 +62,17 @@ int main(int argc, char *argv[])
   long* stabStateIndices;
   int numStabStates;
 
-  srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+  //srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+  srand((unsigned)seed); // seeding the random number generator for sparsify()
 
-  if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+  if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling, samplePrefactor))
     return 1;
 
-  //printf("checking: numStabStateIndices:\n");
-  //for(i=0; i<numStabStates; i++)
-  //  printf("%d ", stabStateIndices[i]);
-  //printf("\n");
-  //fflush(stdout);
+  /* printf("checking: numStabStateIndices:\n"); */
+  /* for(i=0; i<numStabStates; i++) */
+  /*   printf("%ld ", stabStateIndices[i]); */
+  /* printf("\n"); */
+  /* fflush(stdout); */
 
   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
   double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
@@ -248,6 +254,34 @@ int main(int argc, char *argv[])
 
   memcpy(origQ, Q, numStabStates*sizeof(int));
 
+  double complex amplitude = 0.0 + 0.0*I;
+  double complex lastamplitude = 0.0 + 0.0*I;
+
+  double complex probability = 1.0 + 0.0*I;
+
+  // the first measurement should be the identity
+  omega = 0;
+  for(i=0; i<N; i++) {
+    alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+  }
+
+  for(j=0; j<numStabStates; j++) { // the kets
+    Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+  }
+
+  double complex newamplitude;
+  
+  #pragma omp parallel
+  {
+  #pragma omp for schedule(static) collapse(2)
+  for(i=0; i<numStabStates; i++) { // the bras
+    for(j=0; j<numStabStates; j++) {
+      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]);
+      amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+    }
+  }
+  }
+
   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
   
     Paulicounter++;
@@ -273,14 +307,26 @@ int main(int argc, char *argv[])
 
     }
 
-  }
-
-  double complex amplitude = 0.0 + 0.0*I;
-  for(i=0; i<numStabStates; i++) { // the bras
-    for(j=0; j<numStabStates; j++) {
-      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]);
-      amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+    lastamplitude = amplitude;
+    amplitude = 0.0 + 0.0*I;
+    #pragma omp parallel
+    {
+    #pragma omp for schedule(static) collapse(2)
+    for(i=0; i<numStabStates; i++) { // the bras
+      for(j=0; j<numStabStates; j++) {
+	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]);
+	amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+      }
+    }
     }
+    printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
+    printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+
+    probability *= amplitude/lastamplitude;
+    //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
+
+    printf("probability = %lf\n", cabs(probability));
+    printf("!\n");
   }
 
   //printf("numStabStates=%d\n", numStabStates);
@@ -295,6 +341,8 @@ int main(int argc, char *argv[])
   printf("\nabs(amplitude):\n");
   printf("%lf\n", cabs(amplitude));
 
+  printf("probability = %lf\n", cabs(probability));
+
   return 0;
 
 }
diff --git a/weaksim_relerr.c b/weaksim_relerr.c
index 2cc45fb..09bf002 100644
--- a/weaksim_relerr.c
+++ b/weaksim_relerr.c
@@ -6,7 +6,7 @@
 #include <time.h>
 #include "matrix.h"
 #include "exponentialsum.h"
-#include "shrinkstar.h"
+#include "shrink.h"
 #include "extend.h"
 #include "measurepauli.h"
 #include "innerproduct.h"
@@ -23,25 +23,25 @@ int main(int argc, char *argv[])
 {
 
   if(argc != 5) {
-    printf("weaksim_rellerr argument: \"number of stabilizer state samples\" \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+    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");
     exit(0);
   }
 
   int NUMSTABSTATESAMPLES = atoi(argv[1]);        // number of stabilizer state samples
   double additiveErrorDelta = atof(argv[2]);        // additive error delta
   double phi = PI*atof(argv[3]);//PI/4.0; // PI/4.0 is the T gate magic state
-  int coherentSampling = atoi(argv[4]);           // perform coherent sampling (true=1 or false=0)
+  int coherentSampling = atoi(argv[4]);           // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
 
-  int N;                  // number of qubits
+  int N;              // number of qubits
   scanf("%d", &N);
 
   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)
-  scanf("%d", &T);
+  scanf("%d", &T);  
 
   printf("phi = %lf\n", phi);
 
-  int omega[N]; // max of N measurements
-  int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
+  int omega;
+  int alpha[N], beta[N], gamma[N], delta[N];
   int Paulicounter = 0;
 
   int i, j, k, l, m;
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
   FILE *fp;
   char buff[255];
 
-  srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
+  srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate _equatorial()
   
   fp = fopen("Pd.txt", "r");
 
@@ -63,7 +63,7 @@ int main(int argc, char *argv[])
   int PdN = atoi(buff);
 
   Pd = calloc(PdN, sizeof(double*));
-  for(i=0; i<PdN; i++) 
+  for(i=0; i<PdN; i++)
     Pd[i] = calloc(PdN+1, sizeof(double));
 
   double tmp;
@@ -72,8 +72,8 @@ int main(int argc, char *argv[])
     tmp = 0.0;
     for(j=0; j<=i; j++) {
       if(fscanf(fp, "%s", buff) == EOF) {
-	printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
-	return 1;
+  	printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
+  	return 1;
       }
       Pd[i][j] = atof(buff);
       //printf("%e ", Pd[i][j]);
@@ -83,8 +83,10 @@ int main(int argc, char *argv[])
     //printf("total=%f\n", tmp);
   }
 
-  double complex amplitude;
-  
+  //|T> = e^(pi*i/8)/2*sqrt(4-2*sqrt(2))* (sqrt(-i)*(|0>+i|1>)/sqrt(2) + (|0>+|1>)/sqrt(2))
+
+  //double additiveErrorDelta = 0.1;
+
   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>
   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>
   // alternative coefficient to use instead of coeffb to get overall entangled state
@@ -98,15 +100,14 @@ int main(int argc, char *argv[])
 
   srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
 
-
   if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
     return 1;
 
-  //printf("checking: numStabStateIndices:\n");
-  //for(i=0; i<numStabStates; i++)
-  //  printf("%ld ", stabStateIndices[i]);
-  //printf("\n");
-  //fflush(stdout);
+  /* printf("checking: numStabStateIndices:\n"); */
+  /* for(i=0; i<numStabStates; i++) */
+  /*   printf("%ld ", stabStateIndices[i]); */
+  /* printf("\n"); */
+  /* fflush(stdout); */
 
   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
   double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
@@ -117,16 +118,33 @@ int main(int argc, char *argv[])
 
   K = calloc(numStabStates, sizeof(int));
 
-  int origK, origQ, *origD;
-  int **origJ;
-  int **origG, **origGBar;
-  int *origh;
-  double complex origGamma;
+  double complex origGamma[(int)NUMSTABSTATESAMPLES];
+  int *origK, *origQ, **origD, ***origJ;
+  int ***origG, ***origGBar, **origh;
+
+  origG = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origGBar = calloc(NUMSTABSTATESAMPLES,sizeof(int*));
+  origh = calloc(NUMSTABSTATESAMPLES,sizeof(int*));
+  
+  origJ = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origD = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origQ = calloc(NUMSTABSTATESAMPLES,sizeof(int));
+
+  origK = calloc(NUMSTABSTATESAMPLES, sizeof(int));
 
-  long combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+  int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
   
   double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
 
+  printf("!\n");
+  fflush(stdout);
+
+  for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+    origGamma[j] = 1.0+0.0*I;
+    randomstabilizerstate(N, &origK[j], &origh[j], &origG[j], &origGBar[j], &origQ[j], &origD[j], &origJ[j], Pd);
+    //randomstabilizerstate_equatorial(N, &origK[j], &origh[j], &origG[j], &origGBar[j], &origQ[j], &origD[j], &origJ[j]);
+  }
+
+  printf("NUMSTABSTATESAMPLES=%d stabilizer states allocated and initialized!\n", NUMSTABSTATESAMPLES);
+  fflush(stdout);
+
   for(j=0; j<numStabStates; j++) {
 
     combination = stabStateIndices[j];
@@ -137,7 +155,8 @@ int main(int argc, char *argv[])
       K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
       combination /= 2;
     }
-    combination = j;
+    combination = stabStateIndices[j];
+    //origK[j] = K[j];
 
     Gamma[j] = 1.0;
     Gamma[j] *= L1Norm/((double)numStabStates);
@@ -154,8 +173,18 @@ int main(int argc, char *argv[])
 	J[j][k] = calloc(K[j], sizeof(int));
     }
 
+    //origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
+    //origh[j] = calloc(N, sizeof(int));
+
+    //if(K[j] > 0) {
+    //  origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
+    //  for(k=0; k<K[j]; k++)
+    //	origJ[j][k] = calloc(K[j], sizeof(int));
+    //}
+
     for(k=0; k<N; k++) {
       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
+      //origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
     }
 
     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'
@@ -173,27 +202,26 @@ int main(int argc, char *argv[])
 
     for(k=0; k<N; k++) {
 
-      Q[j] += (((combination%2)==1)*Q2 + ((combination%2)==0)*Q1);
+      Q[j] += ((combination%2)==1)*Q2 + ((combination%2)==0)*Q1;
       
-
       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
 
       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
       for(l=0; l<Kcombo; l++) {
-	// D1 has 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
-	switch(combination%2) {
-	case 0:
-	  D[j][Kcounter+l] = D1[l];
-	  break;
-	case 1:
-	  D[j][Kcounter+l] = D2[l];
-	  break;
-	default:
-	  printf("error");
-	  return 1;
+	  // 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
+	  switch(combination%2) {
+	  case 0:
+	    D[j][Kcounter+l] = D1[l];
+	    break;
+	  case 1:
+	    D[j][Kcounter+l] = D2[l];
+	    break;
+	  default:
+	    printf("error");
+	    return 1;
 	  }
 	for(m=0; m<Kcombo; m++) {
-	  // J1 has 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
+	  // 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
 	  switch(combination%2) {
 	  case 0:
 	    J[j][Kcounter+l][Kcounter+m] = J1[l][m];
@@ -209,13 +237,13 @@ int main(int argc, char *argv[])
       }
 
       for(l=0; l<n1; l++) { // assuming n1=n2
-	h[j][k*n1+l] = (((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l]);
+	h[j][k*n1+l] = ((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l];
       }
       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
       for(l=0; l<Kcombo; l++) {
 	for(m=0; m<n1; m++) { // assuming n1=n2
-	  G[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
-	  GBar[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+	  G[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
+	  GBar[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
 	}
       }
       Kcounter = Kcounter + Kcombo;
@@ -228,12 +256,12 @@ int main(int argc, char *argv[])
 
       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
       
-      combination /= 2; // shift to the right by one (in base-7 arithmetic)
+      combination /= 2; // shift to the right by one (in base-2 arithmetic)
     }
     //printf("!\n");
 
     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
-    combination = j;
+    combination = stabStateIndices[j];
     for(k=0; k<(N); k++) {
       Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
       //printf("Kcounter=%d\n", Kcounter);
@@ -244,8 +272,8 @@ int main(int argc, char *argv[])
 	  /* printf("m=%d\n", m); */
 	  /* printf("Kcounter+l=%d\n", Kcounter+l); */
 	  /* printf("k*n1+m=%d\n", k*n1+m); */
-      	  G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
-      	  GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+      	  G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
+      	  GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
       	}
       }
       Kcounter = Kcounter + (n1-Kcombo);
@@ -257,29 +285,56 @@ int main(int argc, char *argv[])
       
       combination /= 2;
     }
+    /* for(k=0; k<N; k++) { */
+    /*   memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int)); */
+    /* } */
+    /* for(k=0; k<K[j]; k++) { */
+    /*   memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));       */
+    /* } */
+
+    /* memcpy(origh[j], h[j], N*sizeof(int)); */
+    /* memcpy(origD[j], D[j], K[j]*sizeof(int)); */
 
-    /*printf("G[%d]:\n", j);
-    printMatrix(G[j], N, N);
-    printf("GBar[%d]:\n", j);
-    printMatrix(GBar[j], N, N);
+  }
+  //exit(0);
+  /* memcpy(origGamma, Gamma, numStabStates*sizeof(double complex)); */
 
-    printf("h[%d]:\n", j);
-    printVector(h[j], N);
+  /* memcpy(origQ, Q, numStabStates*sizeof(int)); */
 
-    printf("J[%d]:\n", j);
-    printMatrix(J[j], K[j], K[j]);
-    
-    printf("D[%d]:\n", j);
-    printVector(D[j], K[j]);
+  printf("numStabStates=%d stabilizer states allocated and initialized!\n", numStabStates);
+  fflush(stdout);
+
+  double complex amplitude = 0.0 + 0.0*I;
+  double complex lastamplitude = 0.0 + 0.0*I;
 
-    printf("Q[%d]=%d\n", j, Q[j]);*/
+  double complex probability = 1.0 + 0.0*I;
 
+  // the first measurement should be the identity
+  omega = 0;
+  for(i=0; i<N; i++) {
+    alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+  }
+
+  for(j=0; j<numStabStates; j++) { // the kets
+    Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
   }
-  //exit(0);
 
-  while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
+  double complex newamplitude;
+  double complex stabstateaverage;
+  for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+    stabstateaverage = 0.0 + 0.0*I;
+    for(i=0; i<numStabStates; i++) { // the bras
+      //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]);
+      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]);
+      stabstateaverage = stabstateaverage + origGamma[j]*Gamma[i]*newamplitude;
+    }
+    amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+  }
 
-    if((Paulicounter+1) > N) {
+  while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
+  
+    Paulicounter++;
+    if(Paulicounter > N) {
       printf("Error: Number of Paulis is greater than N!\n");
       return 1;
     }
@@ -287,116 +342,61 @@ int main(int argc, char *argv[])
     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
     // Y_i = -I*Z*X
     for(i=0; i<N; i++) {
-      if(delta[Paulicounter][i]){
-	omega[Paulicounter] += 3; // -I = I^3
-	beta[Paulicounter][i] = delta[Paulicounter][i];
-	gamma[Paulicounter][i] = delta[Paulicounter][i];
+      if(delta[i]){
+	omega += 3; // -I = I^3
+	beta[i] = delta[i];
+	gamma[i] = delta[i];
       }
     }
 
-    /*printf("*******\n");
-    printf("*******\n");
-    printf("omega=%d\n", omega);
-    printf("X:\n");
-    printVector(gamma, N);
-    printf("Z:\n");
-    printVector(beta, N);
-    printf("*******\n");
-    printf("*******\n");*/
-
-    //for(j=0; j<numStabStates; j++) { // the kets
-
-      /*printf("========\n");
-      printf("before:\n");
-      printf("K=%d\n", K[j]);
-      printf("h:\n");
-      printVector(h[j], N);
-      printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
-      printf("G:\n");
-      printMatrix(G[j], N, N);
-      printf("GBar:\n");
-      printMatrix(GBar[j], N, N);
-      printf("Q=%d\n", Q[j]);
-      printf("D:\n");
-      printVector(D[j], K[j]);
-      printf("J:\n");
-      printMatrix(J[j], K[j], K[j]);*/
-      //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
-      /*printf("\nafter:\n");
-      printf("K=%d\n", K[j]);
-      printf("h:\n");
-      printVector(h[j], N);
-      printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
-      printf("G:\n");
-      printMatrix(G[j], N, N);
-      printf("GBar:\n");
-      printMatrix(GBar[j], N, N);
-      printf("Q=%d\n", Q[j]);
-      printf("D:\n");
-      printVector(D[j], K[j]);
-      printf("J:\n");
-      printMatrix(J[j], K[j], K[j]);*/
 
-    //}
+    for(j=0; j<numStabStates; j++) { // the kets
 
-    Paulicounter++;
-  }
-
-  amplitude = 0.0 + 0.0*I;
-  for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
-    //printf("i=%d\n", i);
+      Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
 
-    randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
+    }
 
-    origGamma = 1.0 + 0.0*I;
-    
-    for(k=0; k<Paulicounter; k++) {
-      origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
-      //printf("k=%d\n", k);
-  }
-    /*printf("origK=%d\n", origK);
-    printf("origG:\n");
-    printMatrix(origG, N, N);
-    printf("origGBar:\n");
-    printMatrix(origGBar, N, N);
-    printf("origh:\n");
-    printVector(origh, N);*/
-
-    double complex stabstateaverage = 0.0 + 0.0*I;
-    
-    for(j=0; j<numStabStates; j++) {
-      //printf("j=%d\n", j);
-      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);
-      stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
+    lastamplitude = amplitude;
+    amplitude = 0.0 + 0.0*I;
+    for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+      stabstateaverage = 0.0 + 0.0*I;
+      for(i=0; i<numStabStates; i++) { // the bras
+	//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]);
+	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]);
+	stabstateaverage = stabstateaverage + origGamma[j]*Gamma[i]*newamplitude;
+      }
+      amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
     }
-    amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+    printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
+    printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+
+    probability *= amplitude/lastamplitude;
+    //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
 
-    deallocate_mem(&origG, N);
-    deallocate_mem(&origGBar, N);
-    free(origh);
-    deallocate_mem(&origJ, origK);
-    free(origD);
+    printf("probability = %lf\n", cabs(probability));
+    printf("!\n");
   }
 
   //printf("numStabStates=%d\n", numStabStates);
   printf("L1Norm=%lf\n", L1Norm);
 
   printf("\namplitude:\n");
-  if(creal(amplitude+ZEROTHRESHOLD)>0)
-    printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+  if(creal(amplitude+0.00000001)>0)
+    printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
   else
-    printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+    printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+  //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
   printf("\nabs(amplitude):\n");
   printf("%lf\n", cabs(amplitude));
 
-
-  for(i=0; i<PdN; i++) 
-    free(Pd[i]);
-  free(Pd);
+  printf("probability = %lf\n", cabs(probability));
 
   return 0;
+
 }
 
+
+
 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 {