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 "sparsify.h"
   9
  10 /****************************************************************************
  11  * SPARSIFY procedure introduced by Bravyi et al.
  12  * (S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and
  13  * M. Howard, Simulation of Quantum Circuits by Low-Rank Sta-
  14  * bilizer Decompositions, Quantum 3, 181 (2019).)
  15  * Outputs stabilizer state indices in stabStateIndices
  16  * as understood by weaksim.c in the tilde'd "computational" basis.
  17  ***Needs address of stabStateIndices since it allocates it.***
  18  * 'numStabStates' is the number of stabilizer state indices,
  19  * 'T' is the number of T-gate magic states, 'phi' is the angle of the
  20  * intermediate diagonal magic state, and 'delta' is the additive
  21  * error of the state.
  22  * For phi=pi/4 this should give you a uniform distro over [0,2^T-1].
  23  ****************************************************************************/
  24
  25 // Note: Make sure you seed the random number generator before calling this!
  26 // i.e. srand((unsigned)time(NULL));
  27 int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling) {
  28
  29   int i, j, k, l;
  30
  31   double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
  32   //printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0));
  33
  34   printf("delta = %lf\n", delta);
  35
  36   //printf("L1Norm = %lf\n", L1Norm);
  37
  38   double gamma;
  39   int increment;
  40
  41   if(coherentSampling) {
  42     printf("Performing coherent sampling!\n");
  43     gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
  44     // since we will be taking a multiple of (T+1) samples, we round numStabStates up to the nearest multiple of (T+1)
  45     *numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
  46     increment = T+1;
  47   } else {
  48     //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
  49     //*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
  50     gamma = 1.0;
  51     *numStabStates = ceil((pow(L1Norm,2.0)-gamma)/pow(delta,2.0));
  52     increment = 1;
  53   }
  54
  55
  56   printf("numStabStates = %d\n", *numStabStates);
  57   fflush(stdout);
  58
  59   double cumL1Norm, sum;
  60   long index, tmpIndex;
  61   int shiftCounter, digit;
  62
  63   *stabStateIndices = calloc(*numStabStates,sizeof(long));
  64
  65   for(i=0; i<*numStabStates; i=i+increment) {
  66     // We will pick our stabilizer index by first selecting how many tilde'd 1s it will have
  67     // by randomly sampling weighted by the binomial expansion of that in the L1norm,
  68     // and then uniformly randomly picking which digits to make that many tilde 1's
  69
  70     // Get random number and scale by the L_1 norm of the t-qubit full stabilizer decomposition
  71     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
  72     sum = 0.0;
  73     //printf("%d cumL1Norm = %lf\n", i, cumL1Norm);
  74     for(j=T; j>=0; j--) {
  75       // add cumulative L1 norm to sum until you are greater than cumL1Norm of random state
  76       //printf("binomcoeff(%d,%d)=%lf\n", T, j, binomcoeff(T,j));
  77       if(j < T)
  78         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);
  79       else
  80         sum += pow((sqrt(1.0-sin(phi))),T);
  81       //printf("%d sum=%lf\n", j, sum);
  82       if(cumL1Norm <= sum) {
  83         // the randomly chosen tilde'd computational state will have j tilde'd 1s
  84         // Now uniformly randomly choose a state with j tilde'd 1s
  85         //printf("%d tilde'd ones\n", (T-j));
  86         index = 0;
  87         for(k=0; k<(T-j); k++) {
  88           digit = (int)((((double)rand())/((double)(RAND_MAX)))*(T-k));
  89           //printf("digit=%d\n", digit);
  90           tmpIndex = index;
  91           shiftCounter = 0;
  92           // first shift past leading 1's
  93           while(tmpIndex%2 == 1) {
  94             tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
  95             shiftCounter++;
  96           }
  97           // now shift past 'digit' 0's
  98           for(l=0; l<digit; l++){
  99             tmpIndex /= 2; // shift to the left by one (base 2 arithmetic)
 100             shiftCounter++;
 101             while(tmpIndex%2 == 1) {
 102               tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
 103               shiftCounter++;
 104             }
 105           }
 106           //printf("shiftCounter = %d\n", shiftCounter);
 107           index += pow(2,shiftCounter);
 108         }
 109         (*stabStateIndices)[i] = index;
 110         //printf("first index = %d\n", index);
 111         if(coherentSampling)
 112           for(k=1; k<=T; k++)
 113             //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
 114             supplement(&(*stabStateIndices)[i], T);
 115         break;
 116       }
 117       //printf("%lf\n", sum);
 118     }
 119     //printf("index:%ld\n ", index);
 120     //printf("(*stabStateIndices)[%d]:%ld\n ", i, (*stabStateIndices)[i]);
 121   }
 122   //printf("\n");
 123
 124   return 0;
 125
 126 }
 127
 128 // binomial coefficient
 129 double binomcoeff(int m, int n) {
 130   //((double)fac(T))/((double)(fac(j)*fac(T-j)))
 131   int i;
 132   double coeff = 1.0;
 133   if (n<=m) {
 134     for(i=n+1;i<=m;i++) {
 135       coeff *= (double)i/((double)i-n);
 136       //printf("%lf ", coeff);
 137     }
 138     //printf("\n");
 139     return coeff;
 140   }
 141   else
 142     return 1;
 143 }