sparsify.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <complex.h>
   5
   6 #include "matrix.h"
   7 #include "supplement.h"
   8 #include "supplement2.h"
   9 #include "sparsify.h"
  10
  11 /****************************************************************************
  12  * SPARSIFY procedure introduced by Bravyi et al.
  13  * (S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and
  14  * M. Howard, Simulation of Quantum Circuits by Low-Rank Sta-
  15  * bilizer Decompositions, Quantum 3, 181 (2019).)
  16  * Outputs stabilizer state indices in stabStateIndices
  17  * as understood by weaksim.c in the tilde'd "computational" basis.
  18  ***Needs address of stabStateIndices since it allocates it.***
  19  * 'numStabStates' is the number of stabilizer state indices,
  20  * 'T' is the number of T-gate magic states, 'phi' is the angle of the
  21  * intermediate diagonal magic state, and 'delta' is the additive
  22  * error of the state.
  23  * For phi=pi/4 this should give you a uniform distro over [0,2^T-1].
  24  ****************************************************************************/
  25
  26 // Note: Make sure you seed the random number generator before calling this!
  27 // i.e. srand((unsigned)time(NULL));
  28 int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor) {
  29
  30   int i, j, k, l;
  31
  32   double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
  33   //printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0));
  34
  35   // fraction of optimal T supplemental states
  36   //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/((double)T);
  37   //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/(4.0+pow(delta,2.0))/((double)T);
  38   double alpha;
  39
  40   // fraction of T supplemental states that we actually need for optimal scaling
  41   // if first term is kept from Seddon et al.
  42   //double fractionSupplement = alpha*0.25*delta*pow(L1Norm,2.0)/((double)T);
  43   // if first term isn't kept
  44   double fractionSupplement;
  45
  46   printf("delta = %lf\n", delta);
  47
  48   printf("L1Norm = %lf\n", L1Norm);
  49
  50   double gamma;
  51   int increment;
  52
  53   if(coherentSampling!=0) {
  54     if(coherentSampling==2) {
  55       alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(2*T-1));
  56       fractionSupplement = alpha*1.0;
  57       *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(2*T-1))+1)))*(floor(fractionSupplement*(2*T-1))+1);
  58       increment = floor(fractionSupplement*(2*T-1))+1;
  59     } else if(coherentSampling==1) {
  60       alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(T));
  61       fractionSupplement = alpha*1.0;
  62       *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(T))+1)))*(floor(fractionSupplement*(T))+1);
  63       increment = floor(fractionSupplement*(T))+1;
  64     }
  65
  66     printf("Performing coherent sampling!\n");
  67     printf("Using an optimal fraction of T supplemental states!\n");
  68     printf("Optimal fraction = %lf\n", fractionSupplement);
  69     // since we will be taking a multiple of (fractionSupplement*T+1) samples, we round numStabStates up to the nearest multiple of (fractionSupplement*T+1)
  70     // if first term is kept from Seddon et al.
  71     //*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);
  72     // if first term isn't kept
  73     //*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);
  74     //*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);
  75
  76     printf("increment= %d\n", increment);
  77     //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)));
  78     //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))));
  79     //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))));
  80     printf("prebinned numStabStates=%d\n", (int)ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta));
  81     fflush(stdout);
  82     if(fractionSupplement > 1.0) {
  83       printf("Too many supplemental states requested!\n");
  84       printf("Note: the correlated sampling is currently only implemented for a maximum of (2t-1) supplemental correlated states!\n");
  85       return 1;
  86     }
  87   } else if(coherentSampling==0) {
  88     //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
  89     // SPARSIFY
  90     //*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
  91     // if first term is kept from Seddon et al.
  92     //*numStabStates = ceil((2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
  93     // if first term isn't kept
  94     *numStabStates = ceil(samplePrefactor*(1.0*2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
  95     // multiply the factor of 2.0 by 1.0 to get the prior state-of-the-art
  96     // multiply the factor of 2.0 by 0.0 to get the tighter bound on the prior state-of-the-art
  97     increment = 1;
  98   } else {
  99     printf("Error: Unknown value for \"coherent sampling\" flag");
 100     return 1;
 101   }
 102
 103
 104   printf("numStabStates = %d\n", *numStabStates);
 105   fflush(stdout);
 106
 107   double cumL1Norm, sum;
 108   long index, tmpIndex;
 109   int shiftCounter, digit;
 110
 111   *stabStateIndices = calloc(*numStabStates,sizeof(long));
 112
 113   for(i=0; i<*numStabStates; i=i+increment) {
 114     // We will pick our stabilizer index by first selecting how many tilde'd 1s it will have
 115     // by randomly sampling weighted by the binomial expansion of that in the L1norm,
 116     // and then uniformly randomly picking which digits to make that many tilde 1's
 117
 118     // Get random number and scale by the L_1 norm of the t-qubit full stabilizer decomposition
 119     cumL1Norm = ((((double)rand())/((double)(RAND_MAX)))*L1Norm); // this is the "cumulative" L1 norm of the randomly selected state if you index by the tilde'd computational basis states truncated to an integer
 120     sum = 0.0;
 121     //printf("%d cumL1Norm = %lf\n", i, cumL1Norm);
 122     for(j=T; j>=0; j--) {
 123       // add cumulative L1 norm to sum until you are greater than cumL1Norm of random state
 124       //printf("binomcoeff(%d,%d)=%lf\n", T, j, binomcoeff(T,j));
 125       if(j < T)
 126         sum += binomcoeff(T,j)*(sqrt(1.0-cos(phi)))*pow(sqrt(1.0 - cos(0.25*PI)),T-1-j)*pow(sqrt(1.0-sin(phi)),j);
 127       else
 128         sum += pow((sqrt(1.0-sin(phi))),T);
 129       //printf("%d sum=%lf\n", j, sum);
 130       if(cumL1Norm <= sum) {
 131         // the randomly chosen tilde'd computational state will have j tilde'd 1s
 132         // Now uniformly randomly choose a state with j tilde'd 1s
 133         //printf("%d tilde'd ones\n", (T-j));
 134         index = 0;
 135         for(k=0; k<(T-j); k++) {
 136           digit = (int)((((double)rand())/((double)(RAND_MAX)))*(T-k));
 137           //printf("digit=%d\n", digit);
 138           tmpIndex = index;
 139           shiftCounter = 0;
 140           // first shift past leading 1's
 141           while(tmpIndex%2 == 1) {
 142             tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
 143             shiftCounter++;
 144           }
 145           // now shift past 'digit' 0's
 146           for(l=0; l<digit; l++){
 147             tmpIndex /= 2; // shift to the left by one (base 2 arithmetic)
 148             shiftCounter++;
 149             while(tmpIndex%2 == 1) {
 150               tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
 151               shiftCounter++;
 152             }
 153           }
 154           //printf("shiftCounter = %d\n", shiftCounter);
 155           index += pow(2,shiftCounter);
 156         }
 157         (*stabStateIndices)[i] = index;
 158         //printf("first index = %d\n", index);
 159         if(coherentSampling==2)
 160           //for(k=1; k<=T; k++)
 161             //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
 162           supplement2(&(*stabStateIndices)[i], T, fractionSupplement);
 163         else if(coherentSampling==1)
 164           //for(k=1; k<=T; k++)
 165             //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
 166           supplement(&(*stabStateIndices)[i], T, fractionSupplement);
 167         break;
 168       }
 169       //printf("%lf\n", sum);
 170     }
 171     //printf("index:%ld\n ", index);
 172     //printf("(*stabStateIndices)[%d]:%ld\n ", i, (*stabStateIndices)[i]);
 173   }
 174   //printf("\n");
 175
 176   return 0;
 177
 178 }
 179
 180 // binomial coefficient
 181 double binomcoeff(int m, int n) {
 182   //((double)fac(T))/((double)(fac(j)*fac(T-j)))
 183   int i;
 184   double coeff = 1.0;
 185   if (n<=m) {
 186     for(i=n+1;i<=m;i++) {
 187       coeff *= (double)i/((double)i-n);
 188       //printf("%lf ", coeff);
 189     }
 190     //printf("\n");
 191     return coeff;
 192   }
 193   else
 194     return 1;
 195 }