7 #include "supplement.h"
8 #include "supplement2.h"
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
23 * For phi=pi/4 this should give you a uniform distro over [0,2^T-1].
24 ****************************************************************************/
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) {
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));
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);
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;
46 printf("delta = %lf\n", delta);
48 printf("L1Norm = %lf\n", L1Norm);
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;
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);
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));
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");
87 } else if(coherentSampling==0) {
88 //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
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
99 printf("Error: Unknown value for \"coherent sampling\" flag");
104 printf("numStabStates = %d\n", *numStabStates);
107 double cumL1Norm, sum;
108 long index, tmpIndex;
109 int shiftCounter, digit;
111 *stabStateIndices = calloc(*numStabStates,sizeof(long));
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
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
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));
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);
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));
135 for(k=0; k<(T-j); k++) {
136 digit = (int)((((double)rand())/((double)(RAND_MAX)))*(T-k));
137 //printf("digit=%d\n", digit);
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
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)
149 while(tmpIndex%2 == 1) {
150 tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
154 //printf("shiftCounter = %d\n", shiftCounter);
155 index += pow(2,shiftCounter);
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);
169 //printf("%lf\n", sum);
171 //printf("index:%ld\n ", index);
172 //printf("(*stabStateIndices)[%d]:%ld\n ", i, (*stabStateIndices)[i]);
180 // binomial coefficient
181 double binomcoeff(int m, int n) {
182 //((double)fac(T))/((double)(fac(j)*fac(T-j)))
186 for(i=n+1;i<=m;i++) {
187 coeff *= (double)i/((double)i-n);
188 //printf("%lf ", coeff);