8 #include "exponentialsum.h"
9 #include "shrinkstar.h"
11 #include "measurepauli.h"
12 #include "innerproduct.h"
13 #include "randomstabilizerstate.h"
15 #define ZEROTHRESHOLD (0.00000001)
17 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
19 // order of matrix elements is [row][column]!!!
21 int main(int argc, char *argv[])
25 printf("strongsim6_rellerr argument: \"number of stabilizer state samples\"\n");
29 int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples
31 int N; // number of qubits
35 printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n");
39 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)
42 int omega[N]; // max of N measurements
43 int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
51 srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
53 fp = fopen("Pd.txt", "r");
55 if(fscanf(fp, "%s", buff) == EOF) {
56 printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
64 Pd = calloc(PdN, sizeof(double*));
66 Pd[i] = calloc(PdN+1, sizeof(double));
70 for(i=1; i<PdN; i++) {
73 if(fscanf(fp, "%s", buff) == EOF) {
74 printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
77 Pd[i][j] = atof(buff);
78 //printf("%e ", Pd[i][j]);
82 //printf("total=%f\n", tmp);
86 double complex coeffb60 = (-16.0+12.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
87 double complex coeffb66 = (96.0-68.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
88 double complex coeffe6 = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)*cpow(2.0,2.5);
89 double complex coeffo6 = (-14.0+10.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-14.0*PI*I/8.0)*cpow(2.0,2.5);
90 double complex coeffk6 = (7.0-5.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-8.0*PI*I/8.0)*4.0*csqrt(2.0)*cpow(2.0,0.5);
91 double complex coeffphiprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
92 double complex coeffphidprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
95 int n1 = 6; int k1 = 6; int (*(G1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h1[] = {0, 0, 0, 0, 0, 0}; int Q1 = 0; int D1[] = {0, 0, 0, 0, 0, 0}; int (*(J1[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
97 int n2 = 6; int k2 = 6; int (*(G2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h2[] = {0, 0, 0, 0, 0, 0}; int Q2 = 4; int D2[] = {4, 4, 4, 4, 4, 4}; int (*(J2[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
99 int n3 = 6; int k3 = 5; int (*(G3[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar3[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h3[] = {1, 0, 0, 0, 0, 0}; int Q3 = 4; int D3[] = {0, 0, 0, 0, 0}; int (*(J3[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
101 int n4 = 6; int k4 = 5; int (*(G4[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar4[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h4[] = {0, 0, 0, 0, 0, 0}; int Q4 = 4; int D4[] = {4, 4, 4, 4, 4}; int (*(J4[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
103 int n5 = 6; int k5 = 1; int (*(G5[])) = { (int[]) {1, 1, 1, 1, 1, 1}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1} }; int (*(GBar5[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1} }; int h5[] = {1, 1, 1, 1, 1, 1}; int Q5 = 6; int D5[] = {2}; int (*(J5[])) = { (int[]) {4} };
105 int n6 = 6; int k6 = 5; int (*(G6[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar6[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h6[] = {1, 0, 0, 0, 0, 0}; int Q6 = 0; int D6[] = {0, 0, 0, 0, 0}; int (*(J6[])) = { (int[]) {0, 4, 0, 0, 4}, (int[]) {4, 0, 4, 0, 0}, (int[]) {0, 4, 0, 4, 0}, (int[]) {0, 0, 4, 0, 4}, (int[]) {4, 0, 0, 4, 0} };
107 int n7 = 6; int k7 = 5; int (*(G7[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar7[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h7[] = {1, 0, 0, 0, 0, 0}; int Q7 = 0; int D7[] = {0, 0, 0, 0, 0}; int (*(J7[])) = { (int[]) {0, 0, 4, 4, 0}, (int[]) {0, 0, 0, 4, 4}, (int[]) {4, 0, 0, 0, 4}, (int[]) {4, 4, 0, 0, 0}, (int[]) {0, 4, 4, 0, 0} };
110 int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
111 double complex Gamma[(int)pow(7,N/6)]; // prefactor in front of resultant state
112 G = calloc(pow(7,N/6),sizeof(int*)); GBar = calloc(pow(7,N/6),sizeof(int*));
113 h = calloc(pow(7,N/6),sizeof(int*));
115 J = calloc(pow(7,N/6),sizeof(int*)); D = calloc(pow(7,N/6),sizeof(int*)); Q = calloc(pow(7,N/6),sizeof(int));
117 K = calloc(pow(7,N/6), sizeof(int));
119 int origK, origQ, *origD;
121 int **origG, **origGBar;
123 double complex origGamma;
125 int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
128 for(j=0; j<pow(7,N/6); j++) { // there will be 7^(N/6) combinations when using k=12 tensor factors
134 for(k=0; k<N/6; k++) {
135 K[j] += (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
142 G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
143 h[j] = calloc(N, sizeof(int));
146 J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
147 for(k=0; k<K[j]; k++)
148 J[j][k] = calloc(K[j], sizeof(int));
152 G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
155 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'
156 int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
157 for(k=0; k<N/6; k++) {
159 Q[j] += (((combination%7)==6)*Q7 + ((combination%7)==5)*Q6 + ((combination%7)==4)*Q5 + ((combination%7)==3)*Q4 + ((combination%7)==2)*Q3 + ((combination%7)==1)*Q2 + ((combination%7)==0)*Q1);
162 Gamma[j] *= (((combination%7)==6)*coeffphidprime + ((combination%7)==5)*coeffphiprime + ((combination%7)==4)*coeffk6 + ((combination%7)==3)*coeffo6 + ((combination%7)==2)*coeffe6 + ((combination%7)==1)*coeffb66 + ((combination%7)==0)*coeffb60);
164 Kcombo = (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
165 for(l=0; l<Kcombo; l++) {
166 // 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
167 switch(combination%7) {
169 D[j][Kcounter+l] = D1[l];
172 D[j][Kcounter+l] = D2[l];
175 D[j][Kcounter+l] = D3[l];
178 D[j][Kcounter+l] = D4[l];
181 D[j][Kcounter+l] = D5[l];
184 D[j][Kcounter+l] = D6[l];
187 D[j][Kcounter+l] = D7[l];
193 for(m=0; m<Kcombo; m++) {
194 // 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
195 switch(combination%7) {
197 J[j][Kcounter+l][Kcounter+m] = J1[l][m];
200 J[j][Kcounter+l][Kcounter+m] = J2[l][m];
203 J[j][Kcounter+l][Kcounter+m] = J3[l][m];
206 J[j][Kcounter+l][Kcounter+m] = J4[l][m];
209 J[j][Kcounter+l][Kcounter+m] = J5[l][m];
212 J[j][Kcounter+l][Kcounter+m] = J6[l][m];
215 J[j][Kcounter+l][Kcounter+m] = J7[l][m];
224 for(l=0; l<n1; l++) { // assuming n1=n2=n3
225 h[j][k*n1+l] = (((combination%7)==6)*h7[l] + ((combination%7)==5)*h6[l] + ((combination%7)==4)*h5[l] + ((combination%7)==3)*h4[l] + ((combination%7)==2)*h3[l] + ((combination%7)==1)*h2[l] + ((combination%7)==0)*h1[l]);
227 // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
228 for(l=0; l<Kcombo; l++) {
229 for(m=0; m<n1; m++) { // assuming n1=n2=n3
230 G[j][Kcounter+l][k*n1+m] = (((combination%7)==6)*G7[l][m] + ((combination%7)==5)*G6[l][m] + ((combination%7)==4)*G5[l][m] + ((combination%7)==3)*G4[l][m] + ((combination%7)==2)*G3[l][m] + ((combination%7)==1)*G2[l][m] + ((combination%7)==0)*G1[l][m]);
231 GBar[j][Kcounter+l][k*n1+m] = (((combination%7)==6)*GBar7[l][m] + ((combination%7)==5)*GBar6[l][m] + ((combination%7)==4)*GBar5[l][m] + ((combination%7)==3)*GBar4[l][m] + ((combination%7)==2)*GBar3[l][m] + ((combination%7)==1)*GBar2[l][m] + ((combination%7)==0)*GBar1[l][m]);
234 Kcounter = Kcounter + Kcombo;
236 /* printf("intermediate G[%d]:\n", j); */
237 /* printMatrix(G[j], N, N); */
238 /* printf("intermediate GBar[%d]:\n", j); */
239 /* printMatrix(GBar[j], N, N); */
240 //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
242 //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
244 combination /= 7; // shift to the right by one (in base-7 arithmetic)
248 // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
250 for(k=0; k<(N/6); k++) {
251 Kcombo = (((combination%7)==6)*k7 + ((combination%7)==5)*k6 + ((combination%7)==4)*k5 + ((combination%7)==3)*k4 + ((combination%7)==2)*k3 + ((combination%7)==1)*k2 + ((combination%7)==0)*k1);
252 //printf("Kcounter=%d\n", Kcounter);
253 // G and GBar rows that are outside the first 'k' spanning basis states
254 for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
255 //printf("l=%d\n", l);
256 for(m=0; m<n1; m++) { // assuming n1=n2=n3
257 /* printf("m=%d\n", m); */
258 /* printf("Kcounter+l=%d\n", Kcounter+l); */
259 /* printf("k*n1+m=%d\n", k*n1+m); */
260 G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%7)==6)*G7[l][m] + ((combination%7)==5)*G6[l][m] + ((combination%7)==4)*G5[l][m] + ((combination%7)==3)*G4[l][m] + ((combination%7)==2)*G3[l][m] + ((combination%7)==1)*G2[l][m] + ((combination%7)==0)*G1[l][m]);
261 GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%7)==6)*GBar7[l][m] + ((combination%7)==5)*GBar6[l][m] + ((combination%7)==4)*GBar5[l][m] + ((combination%7)==3)*GBar4[l][m] + ((combination%7)==2)*GBar3[l][m] + ((combination%7)==1)*GBar2[l][m] + ((combination%7)==0)*GBar1[l][m]);
264 Kcounter = Kcounter + (n1-Kcombo);
266 /* printf("intermediate G[%d]:\n", j); */
267 /* printMatrix(G[j], N, N); */
268 /* printf("intermediate GBar[%d]:\n", j); */
269 /* printMatrix(GBar[j], N, N); */
274 /*printf("G[%d]:\n", j);
275 printMatrix(G[j], N, N);
276 printf("GBar[%d]:\n", j);
277 printMatrix(GBar[j], N, N);
279 printf("h[%d]:\n", j);
280 printVector(h[j], N);
282 printf("J[%d]:\n", j);
283 printMatrix(J[j], K[j], K[j]);
285 printf("D[%d]:\n", j);
286 printVector(D[j], K[j]);
288 printf("Q[%d]=%d\n", j, Q[j]);*/
293 while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
295 if((Paulicounter+1) > N) {
296 printf("Error: Number of Paulis is greater than N!\n");
300 // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
303 if(delta[Paulicounter][i]){
304 omega[Paulicounter] += 3; // -I = I^3
305 beta[Paulicounter][i] = delta[Paulicounter][i];
306 gamma[Paulicounter][i] = delta[Paulicounter][i];
310 /*printf("*******\n");
312 printf("omega=%d\n", omega);
314 printVector(gamma, N);
316 printVector(beta, N);
318 printf("*******\n");*/
320 //for(j=0; j<pow(7,N/6); j++) { // the kets
322 /*printf("========\n");
324 printf("K=%d\n", K[j]);
326 printVector(h[j], N);
327 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
329 printMatrix(G[j], N, N);
331 printMatrix(GBar[j], N, N);
332 printf("Q=%d\n", Q[j]);
334 printVector(D[j], K[j]);
336 printMatrix(J[j], K[j], K[j]);*/
337 //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
338 /*printf("\nafter:\n");
339 printf("K=%d\n", K[j]);
341 printVector(h[j], N);
342 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
344 printMatrix(G[j], N, N);
346 printMatrix(GBar[j], N, N);
347 printf("Q=%d\n", Q[j]);
349 printVector(D[j], K[j]);
351 printMatrix(J[j], K[j], K[j]);*/
358 double complex amplitude = 0.0 + 0.0*I;
359 for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
360 //printf("i=%d\n", i);
362 randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
364 origGamma = 1.0 + 0.0*I;
366 for(k=0; k<Paulicounter; k++) {
367 origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
368 //printf("k=%d\n", k);
370 /*printf("origK=%d\n", origK);
372 printMatrix(origG, N, N);
373 printf("origGBar:\n");
374 printMatrix(origGBar, N, N);
376 printVector(origh, N);*/
378 double complex stabstateaverage = 0.0 + 0.0*I;
380 for(j=0; j<pow(7,N/6); j++) {
381 //printf("j=%d\n", j);
382 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);
383 stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
385 amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
387 deallocate_mem(&origG, N);
388 deallocate_mem(&origGBar, N);
390 deallocate_mem(&origJ, origK);
394 printf("amplitude:\n");
395 if(creal(amplitude+ZEROTHRESHOLD)>0)
396 printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
398 printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
409 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
412 int newomega, newalpha, newbeta, newgamma, newdelta;
415 if(scanf("%d", &newomega) != EOF) {
417 for(i=0; i<numqubits; i++) {
418 if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
419 printf("Error: Too few input coeffs!\n");
422 if(newalpha+newbeta+newgamma+newdelta > 1) {
423 printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
426 alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;