X-Git-Url: https://s3miclassical.com/gitweb/?p=weak_simulation_stab_extent.git;a=blobdiff_plain;f=sparsify.c;h=268dd983683bf249e57ccbe5d6eb21df6ad19e49;hp=8e164fc35f669969c9b7c20419c5bade94c19b10;hb=7bc25656311c391af6a4d5814cc82080981c63c4;hpb=d914e7c183347b1340988e0869a6e053dd573991 diff --git a/sparsify.c b/sparsify.c index 8e164fc..268dd98 100644 --- a/sparsify.c +++ b/sparsify.c @@ -5,6 +5,7 @@ #include "matrix.h" #include "supplement.h" +#include "supplement2.h" #include "sparsify.h" /**************************************************************************** @@ -24,32 +25,79 @@ // Note: Make sure you seed the random number generator before calling this! // i.e. srand((unsigned)time(NULL)); -int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling) { +int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor) { int i, j, k, l; double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T); //printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0)); + // fraction of optimal T supplemental states + //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/((double)T); + //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/(4.0+pow(delta,2.0))/((double)T); + double alpha; + + // fraction of T supplemental states that we actually need for optimal scaling + // if first term is kept from Seddon et al. + //double fractionSupplement = alpha*0.25*delta*pow(L1Norm,2.0)/((double)T); + // if first term isn't kept + double fractionSupplement; + printf("delta = %lf\n", delta); - //printf("L1Norm = %lf\n", L1Norm); + printf("L1Norm = %lf\n", L1Norm); double gamma; int increment; - if(coherentSampling) { + if(coherentSampling!=0) { + if(coherentSampling==2) { + alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(2*T-1)); + fractionSupplement = alpha*1.0; + *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(2*T-1))+1)))*(floor(fractionSupplement*(2*T-1))+1); + increment = floor(fractionSupplement*(2*T-1))+1; + } else if(coherentSampling==1) { + alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(T)); + fractionSupplement = alpha*1.0; + *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(T))+1)))*(floor(fractionSupplement*(T))+1); + increment = floor(fractionSupplement*(T))+1; + } + printf("Performing coherent sampling!\n"); - gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T; - // since we will be taking a multiple of (T+1) samples, we round numStabStates up to the nearest multiple of (T+1) - *numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1); - increment = T+1; - } else { + printf("Using an optimal fraction of T supplemental states!\n"); + printf("Optimal fraction = %lf\n", fractionSupplement); + // since we will be taking a multiple of (fractionSupplement*T+1) samples, we round numStabStates up to the nearest multiple of (fractionSupplement*T+1) + // if first term is kept from Seddon et al. + //*numStabStates = ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + // if first term isn't kept + //*numStabStates = ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + //*numStabStates = ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + + printf("increment= %d\n", increment); + //printf("prebinned numStabStates=%d\n", (int)ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0))); + //printf("prebinned numStabStates=%d\n", (int)ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0)))); + //printf("prebinned numStabStates=%d\n", (int)ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0)))); + printf("prebinned numStabStates=%d\n", (int)ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta)); + fflush(stdout); + if(fractionSupplement > 1.0) { + printf("Too many supplemental states requested!\n"); + printf("Note: the correlated sampling is currently only implemented for a maximum of (2t-1) supplemental correlated states!\n"); + return 1; + } + } else if(coherentSampling==0) { //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T; + // SPARSIFY //*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1); - gamma = 1.0; - *numStabStates = ceil((pow(L1Norm,2.0)-gamma)/pow(delta,2.0)); + // if first term is kept from Seddon et al. + //*numStabStates = ceil((2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0)); + // if first term isn't kept + *numStabStates = ceil(samplePrefactor*(1.0*2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0)); + // multiply the factor of 2.0 by 1.0 to get the prior state-of-the-art + // multiply the factor of 2.0 by 0.0 to get the tighter bound on the prior state-of-the-art increment = 1; + } else { + printf("Error: Unknown value for \"coherent sampling\" flag"); + return 1; } @@ -108,10 +156,14 @@ int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, dou } (*stabStateIndices)[i] = index; //printf("first index = %d\n", index); - if(coherentSampling) - for(k=1; k<=T; k++) + if(coherentSampling==2) + //for(k=1; k<=T; k++) + //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index + supplement2(&(*stabStateIndices)[i], T, fractionSupplement); + else if(coherentSampling==1) + //for(k=1; k<=T; k++) //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index - supplement(&(*stabStateIndices)[i], T); + supplement(&(*stabStateIndices)[i], T, fractionSupplement); break; } //printf("%lf\n", sum);