From: Lucas K <lkocia@s3miclassical.com>
Date: Tue, 22 Dec 2020 22:56:20 +0000 (-0800)
Subject: delete specific instances of strongsim.c and strongsim_relerr.c
X-Git-Url: https://s3miclassical.com/gitweb/?p=strong_simulation_stabilizer_rank.git;a=commitdiff_plain;h=00f8f5c3f25b6137aec8180d0da144608ae55555

delete specific instances of strongsim.c and strongsim_relerr.c
---

diff --git a/strongsim.c b/strongsim.c
deleted file mode 100644
index 389a737..0000000
--- a/strongsim.c
+++ /dev/null
@@ -1,286 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include <complex.h>
-#include "matrix.h"
-#include "exponentialsum.h"
-#include "shrink.h"
-#include "extend.h"
-#include "measurepauli.h"
-#include "innerproduct.h"
-
-int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
-
-// order of matrix elements is [row][column]!!!
-
-int main()
-{
-
-  int N;              // number of qubits
-  scanf("%d", &N);
-
-  if(N%3 != 0) {
-    printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n");
-    return 1;
-  }
-
-  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);  
-
-  int omega;
-  int alpha[N], beta[N], gamma[N], delta[N];
-  int Paulicounter = 0;
-
-  int i, j, k, l, m;
-
-  double complex coeffa = -0.25*(1.0-I)*(-1.0-I+sqrt(2.0))*csqrt(-I);
-  double complex coeffb = 0.25*(-1.0-I)*(1.0-I+sqrt(2.0))*csqrt(I);
-  double complex coeffc = 0.25*(-1.0-I)*(-1.0+I+sqrt(2.0))*csqrt(I);
-
-  int n1 = 3; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1, 1}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0}, (int[]) {1, 1, 0}, (int[]) {1, 0, 1}}; int h1[] = {1, 1, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
-  int n2 = 3; int k2 = 3; int (*(G2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h2[] = {0, 0, 0}; int Q2 = 2; int D2[] = {2, 2, 0}; int (*(J2[])) = { (int[]) {4, 0, 0}, (int[]) {0, 4, 0}, (int[]) {0, 0, 0} };
-  int n3 = 3; int k3 = 3; int (*(G3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h3[] = {0, 0, 0}; int Q3 = 2; int D3[] = {6, 6, 0}; int (*(J3[])) = { (int[]) {4, 4, 4}, (int[]) {4, 4, 4}, (int[]) {4, 4, 0} };
-
-  int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
-  double complex Gamma[(int)pow(3,N/3)]; // prefactor in front of resultant state
-  G = calloc(pow(3,N/3),sizeof(int*)); GBar = calloc(pow(3,N/3),sizeof(int*));
-  h = calloc(pow(3,N/3),sizeof(int*));
-  
-  J = calloc(pow(3,N/3),sizeof(int*)); D = calloc(pow(3,N/3),sizeof(int*)); Q = calloc(pow(3,N/3),sizeof(int));
-
-  K = calloc(pow(3,N/3), sizeof(int));
-
-  double complex origGamma[(int)pow(3,N/3)];
-  int *origK, *origQ, **origD, ***origJ;
-  int ***origG, ***origGBar, **origh;
-
-  origG = calloc(pow(3,N/3),sizeof(int*)); origGBar = calloc(pow(3,N/3),sizeof(int*));
-  origh = calloc(pow(3,N/3),sizeof(int*));
-  
-  origJ = calloc(pow(3,N/3),sizeof(int*)); origD = calloc(pow(3,N/3),sizeof(int*)); origQ = calloc(pow(3,N/3),sizeof(int));
-
-  origK = calloc(pow(3,N/3), sizeof(int));
-
-  int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
-  
-
-  for(j=0; j<pow(3,N/3); j++) { // there will be 3^(N/3) combinations when using k=3 tensor factors
-
-    combination = j;
-
-    K[j] = 0.0;
-    
-    for(k=0; k<N/3; k++) {
-      K[j] += (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      combination /= 3;
-    }
-    combination = j;
-    origK[j] = K[j];
-
-    Gamma[j] = 1.0;
-
-    G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
-    h[j] = calloc(N, sizeof(int));
-
-    if(K[j] > 0) {
-      J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
-      for(k=0; k<K[j]; k++)
-	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'
-    int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
-    for(k=0; k<N/3; k++) {
-
-      Q[j] += ((combination%3)==2)*Q3 + ((combination%3)==1)*Q2 + ((combination%3)==0)*Q1;
-      
-      Gamma[j] *= (((combination%3)==2)*coeffc + ((combination%3)==1)*coeffb + ((combination%3)==0)*coeffa);
-
-      Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      for(l=0; l<Kcombo; l++) {
-	  // 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
-	  switch(combination%3) {
-	  case 0:
-	    D[j][Kcounter+l] = D1[l];
-	    break;
-	  case 1:
-	    D[j][Kcounter+l] = D2[l];
-	    break;
-	  case 2:
-	    D[j][Kcounter+l] = D3[l];
-	    break;
-	  default:
-	    printf("error");
-	    return 1;
-	  }
-	for(m=0; m<Kcombo; m++) {
-	  // 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
-	  switch(combination%3) {
-	  case 0:
-	    J[j][Kcounter+l][Kcounter+m] = J1[l][m];
-	    break;
-	  case 1:
-	    J[j][Kcounter+l][Kcounter+m] = J2[l][m];
-	    break;
-	  case 2:
-	    J[j][Kcounter+l][Kcounter+m] = J3[l][m];
-	    break;
-	  default:
-	    printf("error");
-	    return 1;
-	  }
-	}
-      }
-
-      for(l=0; l<n1; l++) { // assuming n1=n2=n3
-	h[j][k*n1+l] = ((combination%3)==2)*h3[l] + ((combination%3)==1)*h2[l] + ((combination%3)==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=n3
-	  G[j][Kcounter+l][k*n1+m] = ((combination%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m];
-	  GBar[j][Kcounter+l][k*n1+m] = ((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m];
-	}
-      }
-      Kcounter = Kcounter + Kcombo;
-      
-      /* printf("intermediate G[%d]:\n", j); */
-      /* printMatrix(G[j], N, N); */
-      /* printf("intermediate GBar[%d]:\n", j); */
-      /* printMatrix(GBar[j], N, N); */
-      //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
-
-      //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
-      
-      combination /= 3; // shift to the right by one (in base-3 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;
-    for(k=0; k<(N/3); k++) {
-      Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      //printf("Kcounter=%d\n", Kcounter);
-      // G and GBar rows that are outside the first 'k' spanning basis states
-      for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
-	//printf("l=%d\n", l);
-      	for(m=0; m<n1; m++) { // assuming n1=n2=n3
-	  /* 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%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m];
-      	  GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m];
-      	}
-      }
-      Kcounter = Kcounter + (n1-Kcombo);
-
-      /* printf("intermediate G[%d]:\n", j); */
-      /* printMatrix(G[j], N, N); */
-      /* printf("intermediate GBar[%d]:\n", j); */
-      /* printMatrix(GBar[j], N, N); */
-      
-      combination /= 3;
-    }
-    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));
-
-  }
-  //exit(0);
-  memcpy(origGamma, Gamma, pow(3,N/3)*sizeof(double complex));
-
-  memcpy(origQ, Q, pow(3,N/3)*sizeof(int));
-
-  while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
-  
-    Paulicounter++;
-    if(Paulicounter > N) {
-      printf("Error: Number of Paulis is greater than N!\n");
-      return 1;
-    }
-    
-    // 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[i]){
-	omega += 3; // -I = I^3
-	beta[i] = delta[i];
-	gamma[i] = delta[i];
-      }
-    }
-
-
-    for(j=0; j<pow(3,N/3); 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 amplitude = 0.0 + 0.0*I;
-  for(i=0; i<pow(3,N/3); i++) { // the bras
-    for(j=0; j<pow(3,N/3); 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;
-    }
-  }
-
-  printf("amplitude:\n");
-  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("%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)));
-
-  return 0;
-
-}
-
-
-
-int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
-{
-    
-  int newomega, newalpha, newbeta, newgamma, newdelta;
-  int i;
-
-  if(scanf("%d", &newomega) != EOF) {
-    *omega = newomega;
-    for(i=0; i<numqubits; i++) {
-      if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
-	printf("Error: Too few input coeffs!\n");
-	exit(0);
-      }
-      if(newalpha+newbeta+newgamma+newdelta > 1) {
-	printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
-	exit(0);
-      }
-      alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
-    }
-    return 1;
-  } else
-    return 0;
-    
-}
diff --git a/strongsim_relerr.c b/strongsim_relerr.c
deleted file mode 100644
index 24ec5b3..0000000
--- a/strongsim_relerr.c
+++ /dev/null
@@ -1,393 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include <complex.h>
-#include <time.h>
-#include "matrix.h"
-#include "exponentialsum.h"
-#include "shrinkstar.h"
-#include "extend.h"
-#include "measurepauli.h"
-#include "innerproduct.h"
-#include "randomstabilizerstate.h"
-
-#define ZEROTHRESHOLD (0.00000001)
-
-int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
-
-// order of matrix elements is [row][column]!!!
-
-int main(int argc, char *argv[])
-{
-
-  if(argc != 2) {
-    printf("strongsim3_rellerr argument: \"number of stabilizer state samples\"\n");
-    exit(0);
-  }
-
-  int NUMSTABSTATESAMPLES = atoi(argv[1]);        // number of stabilizer state samples
-
-  int N;                  // number of qubits
-  scanf("%d", &N);
-
-  if(N%3 != 0) {
-    printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n");
-    return 1;
-  }
-
-  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);
-
-  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 Paulicounter = 0;
-
-  int i, j, k, l, m;
-
-  FILE *fp;
-  char buff[255];
-
-  srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
-  
-  fp = fopen("Pd.txt", "r");
-
-  if(fscanf(fp, "%s", buff) == EOF) {
-    printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
-    return 1;
-  }
-  
-  double** Pd;
-
-  int PdN = atoi(buff);
-
-  Pd = calloc(PdN, sizeof(double*));
-  for(i=0; i<PdN; i++) 
-    Pd[i] = calloc(PdN+1, sizeof(double));
-
-  double tmp;
-  
-  for(i=1; i<PdN; i++) {
-    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;
-      }
-      Pd[i][j] = atof(buff);
-      //printf("%e ", Pd[i][j]);
-      tmp += Pd[i][j];
-    }
-    //printf("\n");
-    //printf("total=%f\n", tmp);
-  }
-
-
-  double complex coeffa = -0.25*(1.0-I)*(-1.0-I+sqrt(2.0))*csqrt(-I);
-  double complex coeffb = 0.25*(-1.0-I)*(1.0-I+sqrt(2.0))*csqrt(I);
-  double complex coeffc = 0.25*(-1.0-I)*(-1.0+I+sqrt(2.0))*csqrt(I);
-
-  int n1 = 3; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1, 1}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0}, (int[]) {1, 1, 0}, (int[]) {1, 0, 1}}; int h1[] = {1, 1, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
-  int n2 = 3; int k2 = 3; int (*(G2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h2[] = {0, 0, 0}; int Q2 = 2; int D2[] = {2, 2, 0}; int (*(J2[])) = { (int[]) {4, 0, 0}, (int[]) {0, 4, 0}, (int[]) {0, 0, 0} };
-  int n3 = 3; int k3 = 3; int (*(G3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h3[] = {0, 0, 0}; int Q3 = 2; int D3[] = {6, 6, 0}; int (*(J3[])) = { (int[]) {4, 4, 4}, (int[]) {4, 4, 4}, (int[]) {4, 4, 0} };
-
-
-  int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
-  double complex Gamma[(int)pow(3,N/3)]; // prefactor in front of resultant state
-  G = calloc(pow(3,N/3),sizeof(int*)); GBar = calloc(pow(3,N/3),sizeof(int*));
-  h = calloc(pow(3,N/3),sizeof(int*));
-  
-  J = calloc(pow(3,N/3),sizeof(int*)); D = calloc(pow(3,N/3),sizeof(int*)); Q = calloc(pow(3,N/3),sizeof(int));
-
-  K = calloc(pow(3,N/3), sizeof(int));
-
-  int origK, origQ, *origD;
-  int **origJ;
-  int **origG, **origGBar;
-  int *origh;
-  double complex origGamma;
-
-  int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
-  
-
-  for(j=0; j<pow(3,N/3); j++) { // there will be 3^(N/3) combinations when using k=12 tensor factors
-
-    combination = j;
-
-    K[j] = 0.0;
-    
-    for(k=0; k<N/3; k++) {
-      K[j] += (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      combination /= 3;
-    }
-    combination = j;
-
-    Gamma[j] = 1.0;
-
-    G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
-    h[j] = calloc(N, sizeof(int));
-
-    if(K[j] > 0) {
-      J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
-      for(k=0; k<K[j]; k++)
-	J[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));
-    }
-
-    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'
-    int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
-    for(k=0; k<N/3; k++) {
-
-      Q[j] += (((combination%3)==2)*Q3 + ((combination%3)==1)*Q2 + ((combination%3)==0)*Q1);
-      
-
-      Gamma[j] *= (((combination%3)==2)*coeffc + ((combination%3)==1)*coeffb + ((combination%3)==0)*coeffa);
-
-      Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      for(l=0; l<Kcombo; l++) {
-	// 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
-	switch(combination%3) {
-	case 0:
-	  D[j][Kcounter+l] = D1[l];
-	  break;
-	case 1:
-	  D[j][Kcounter+l] = D2[l];
-	  break;
-	case 2:
-	  D[j][Kcounter+l] = D3[l];
-	  break;
-	default:
-	  printf("error");
-	  return 1;
-	  }
-	for(m=0; m<Kcombo; m++) {
-	  // 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
-	  switch(combination%3) {
-	  case 0:
-	    J[j][Kcounter+l][Kcounter+m] = J1[l][m];
-	    break;
-	  case 1:
-	    J[j][Kcounter+l][Kcounter+m] = J2[l][m];
-	    break;
-	  case 2:
-	    J[j][Kcounter+l][Kcounter+m] = J3[l][m];
-	    break;
-	  default:
-	    printf("error");
-	    return 1;
-	  }
-	}
-      }
-
-      for(l=0; l<n1; l++) { // assuming n1=n2=n3
-	h[j][k*n1+l] = (((combination%3)==2)*h3[l] + ((combination%3)==1)*h2[l] + ((combination%3)==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=n3
-	  G[j][Kcounter+l][k*n1+m] = (((combination%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m]);
-	  GBar[j][Kcounter+l][k*n1+m] = (((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m]);
-	}
-      }
-      Kcounter = Kcounter + Kcombo;
-      
-      /* printf("intermediate G[%d]:\n", j); */
-      /* printMatrix(G[j], N, N); */
-      /* printf("intermediate GBar[%d]:\n", j); */
-      /* printMatrix(GBar[j], N, N); */
-      //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
-
-      //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
-      
-      combination /= 3; // shift to the right by one (in base-7 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;
-    for(k=0; k<(N/3); k++) {
-      Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
-      //printf("Kcounter=%d\n", Kcounter);
-      // G and GBar rows that are outside the first 'k' spanning basis states
-      for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
-	//printf("l=%d\n", l);
-      	for(m=0; m<n1; m++) { // assuming n1=n2=n3
-	  /* 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%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m]);
-      	  GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m]);
-      	}
-      }
-      Kcounter = Kcounter + (n1-Kcombo);
-
-      /* printf("intermediate G[%d]:\n", j); */
-      /* printMatrix(G[j], N, N); */
-      /* printf("intermediate GBar[%d]:\n", j); */
-      /* printMatrix(GBar[j], N, N); */
-      
-      combination /= 3;
-    }
-
-    /*printf("G[%d]:\n", j);
-    printMatrix(G[j], N, N);
-    printf("GBar[%d]:\n", j);
-    printMatrix(GBar[j], N, N);
-
-    printf("h[%d]:\n", j);
-    printVector(h[j], N);
-
-    printf("J[%d]:\n", j);
-    printMatrix(J[j], K[j], K[j]);
-    
-    printf("D[%d]:\n", j);
-    printVector(D[j], K[j]);
-
-    printf("Q[%d]=%d\n", j, Q[j]);*/
-
-  }
-  //exit(0);
-
-  while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
-
-    if((Paulicounter+1) > N) {
-      printf("Error: Number of Paulis is greater than N!\n");
-      return 1;
-    }
-    
-    // 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];
-      }
-    }
-
-    /*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<pow(3,N/3); 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]);*/
-
-    //}
-
-    Paulicounter++;
-  }
-
-  double complex amplitude = 0.0 + 0.0*I;
-  for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
-    //printf("i=%d\n", i);
-
-    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<pow(3,N/3); 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;
-    }
-    amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
-
-    deallocate_mem(&origG, N);
-    deallocate_mem(&origGBar, N);
-    free(origh);
-    deallocate_mem(&origJ, origK);
-    free(origD);
-  }
-
-  printf("amplitude:\n");
-  if(creal(amplitude+ZEROTHRESHOLD)>0)
-    printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
-  else
-    printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
-  
-  
-
-  for(i=0; i<PdN; i++) 
-    free(Pd[i]);
-  free(Pd);
-
-  return 0;
-}
-
-int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
-{
-    
-  int newomega, newalpha, newbeta, newgamma, newdelta;
-  int i;
-
-  if(scanf("%d", &newomega) != EOF) {
-    *omega = newomega;
-    for(i=0; i<numqubits; i++) {
-      if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
-	printf("Error: Too few input coeffs!\n");
-	exit(0);
-      }
-      if(newalpha+newbeta+newgamma+newdelta > 1) {
-	printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
-	exit(0);
-      }
-      alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
-    }
-    return 1;
-  } else
-    return 0;
-    
-}