+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+#include <time.h>
+#include "matrix.h"
+#include "exponentialsum.h"
+#include "shrinkstar.h"
+#include "extend.h"
+#include "measurepauli.h"
+#include "innerproduct.h"
+#include "randomstabilizerstate.h"
+#include "sparsify.h"
+
+#define ZEROTHRESHOLD (0.00000001)
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+int main(int argc, char *argv[])
+{
+
+ if(argc != 5) {
+ printf("weaksim_rellerr argument: \"number of stabilizer state samples\" \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+ exit(0);
+ }
+
+ int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples
+ double additiveErrorDelta = atof(argv[2]); // additive error delta
+ double phi = PI*atof(argv[3]);//PI/4.0; // PI/4.0 is the T gate magic state
+ int coherentSampling = atoi(argv[4]); // perform coherent sampling (true=1 or false=0)
+
+ int N; // number of qubits
+ scanf("%d", &N);
+
+ 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)
+ scanf("%d", &T);
+
+ printf("phi = %lf\n", phi);
+
+ int omega[N]; // max of N measurements
+ int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
+ int Paulicounter = 0;
+
+ int i, j, k, l, m;
+
+ FILE *fp;
+ char buff[255];
+
+ srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
+
+ fp = fopen("Pd.txt", "r");
+
+ if(fscanf(fp, "%s", buff) == EOF) {
+ printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
+ return 1;
+ }
+
+ double** Pd;
+
+ int PdN = atoi(buff);
+
+ Pd = calloc(PdN, sizeof(double*));
+ for(i=0; i<PdN; i++)
+ Pd[i] = calloc(PdN+1, sizeof(double));
+
+ double tmp;
+
+ for(i=1; i<PdN; i++) {
+ tmp = 0.0;
+ for(j=0; j<=i; j++) {
+ if(fscanf(fp, "%s", buff) == EOF) {
+ printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
+ return 1;
+ }
+ Pd[i][j] = atof(buff);
+ //printf("%e ", Pd[i][j]);
+ tmp += Pd[i][j];
+ }
+ //printf("\n");
+ //printf("total=%f\n", tmp);
+ }
+
+ double complex amplitude;
+
+ double complex coeffa = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*cexp(-PI*I*0.25)*I/sqrt(2.0)*(-I+cexp(-0.25*PI*I))*(-I+cexp(I*phi)))); // factor of cexp(PI*I/8.0)*cexp(-PI*I*0.25) comes from converting (|0>+|1>)/sqrt(2) under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ double complex coeffb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*phi)))); // factor of cexp(PI*I/8.0) comes from converting |0> under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ // alternative coefficient to use instead of coeffb to get overall entangled state
+ double complex coeffbb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*0.25*PI))));
+
+ int n1 = 1; int k1 = 1; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
+ int n2 = 1; int k2 = 1; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {0}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
+
+ long* stabStateIndices;
+ int numStabStates;
+
+ srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+
+
+ if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+ return 1;
+
+ //printf("checking: numStabStateIndices:\n");
+ //for(i=0; i<numStabStates; i++)
+ // printf("%ld ", stabStateIndices[i]);
+ //printf("\n");
+ //fflush(stdout);
+
+ int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
+ double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
+ G = calloc(numStabStates,sizeof(int*)); GBar = calloc(numStabStates,sizeof(int*));
+ h = calloc(numStabStates,sizeof(int*));
+
+ J = calloc(numStabStates,sizeof(int*)); D = calloc(numStabStates,sizeof(int*)); Q = calloc(numStabStates,sizeof(int));
+
+ K = calloc(numStabStates, sizeof(int));
+
+ int origK, origQ, *origD;
+ int **origJ;
+ int **origG, **origGBar;
+ int *origh;
+ double complex origGamma;
+
+ long combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+
+ double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
+
+ for(j=0; j<numStabStates; j++) {
+
+ combination = stabStateIndices[j];
+
+ K[j] = 0.0;
+
+ for(k=0; k<N; k++) {
+ K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ combination /= 2;
+ }
+ combination = j;
+
+ Gamma[j] = 1.0;
+ Gamma[j] *= L1Norm/((double)numStabStates);
+
+ // the coefficients which are a product of 'coeffa', 'coeffb', 'coeffbb' (that are subsequently multiplied into Gamma[j]) is multiplied by 'norm'
+ //Gamma[j] *= norm;
+
+ G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
+ h[j] = calloc(N, sizeof(int));
+
+ if(K[j] > 0) {
+ J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
+ for(k=0; k<K[j]; k++)
+ J[j][k] = calloc(K[j], sizeof(int));
+ }
+
+ for(k=0; k<N; k++) {
+ G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
+ }
+
+ 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'
+ int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
+
+ // if combination contains at least one instance of the second state, i.e. contains the 0 digit in binary, then we want to have it have one instance of coeffb instead of coeffbb
+ for(k=0; k<N; k++) {
+ if(combination%2==1) {
+ Gamma[j] *= coeffb/coeffbb;
+ break; // break out of loop
+ }
+ combination /= 2; // shift to the right by one (in base-2 arithmetic)
+ }
+ combination = stabStateIndices[j];
+
+ for(k=0; k<N; k++) {
+
+ Q[j] += (((combination%2)==1)*Q2 + ((combination%2)==0)*Q1);
+
+
+ Gamma[j] *= (((combination%2)==1)*coeffbb + ((combination%2)==0)*coeffa); // only assign coeffbb instead of coeffb; coeffb replaces one instance of coeffbb before this loop
+
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ for(l=0; l<Kcombo; l++) {
+ // D1 has a different number of rows 'l' than D2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ D[j][Kcounter+l] = D1[l];
+ break;
+ case 1:
+ D[j][Kcounter+l] = D2[l];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ for(m=0; m<Kcombo; m++) {
+ // J1 has a different number of rows 'l' than J2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ J[j][Kcounter+l][Kcounter+m] = J1[l][m];
+ break;
+ case 1:
+ J[j][Kcounter+l][Kcounter+m] = J2[l][m];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ }
+ }
+
+ for(l=0; l<n1; l++) { // assuming n1=n2
+ h[j][k*n1+l] = (((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l]);
+ }
+ // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
+ for(l=0; l<Kcombo; l++) {
+ for(m=0; m<n1; m++) { // assuming n1=n2
+ G[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
+ GBar[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+ }
+ }
+ Kcounter = Kcounter + Kcombo;
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+ //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
+
+ //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
+
+ combination /= 2; // shift to the right by one (in base-7 arithmetic)
+ }
+ //printf("!\n");
+
+ // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
+ combination = j;
+ for(k=0; k<(N); k++) {
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ //printf("Kcounter=%d\n", Kcounter);
+ // G and GBar rows that are outside the first 'k' spanning basis states
+ for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
+ //printf("l=%d\n", l);
+ for(m=0; m<n1; m++) { // assuming n1=n2=n3
+ /* printf("m=%d\n", m); */
+ /* printf("Kcounter+l=%d\n", Kcounter+l); */
+ /* printf("k*n1+m=%d\n", k*n1+m); */
+ G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
+ GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+ }
+ }
+ Kcounter = Kcounter + (n1-Kcombo);
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+
+ combination /= 2;
+ }
+
+ /*printf("G[%d]:\n", j);
+ printMatrix(G[j], N, N);
+ printf("GBar[%d]:\n", j);
+ printMatrix(GBar[j], N, N);
+
+ printf("h[%d]:\n", j);
+ printVector(h[j], N);
+
+ printf("J[%d]:\n", j);
+ printMatrix(J[j], K[j], K[j]);
+
+ printf("D[%d]:\n", j);
+ printVector(D[j], K[j]);
+
+ printf("Q[%d]=%d\n", j, Q[j]);*/
+
+ }
+ //exit(0);
+
+ while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
+
+ if((Paulicounter+1) > N) {
+ printf("Error: Number of Paulis is greater than N!\n");
+ return 1;
+ }
+
+ // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
+ // Y_i = -I*Z*X
+ for(i=0; i<N; i++) {
+ if(delta[Paulicounter][i]){
+ omega[Paulicounter] += 3; // -I = I^3
+ beta[Paulicounter][i] = delta[Paulicounter][i];
+ gamma[Paulicounter][i] = delta[Paulicounter][i];
+ }
+ }
+
+ /*printf("*******\n");
+ printf("*******\n");
+ printf("omega=%d\n", omega);
+ printf("X:\n");
+ printVector(gamma, N);
+ printf("Z:\n");
+ printVector(beta, N);
+ printf("*******\n");
+ printf("*******\n");*/
+
+ //for(j=0; j<numStabStates; j++) { // the kets
+
+ /*printf("========\n");
+ printf("before:\n");
+ printf("K=%d\n", K[j]);
+ printf("h:\n");
+ printVector(h[j], N);
+ printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
+ printf("G:\n");
+ printMatrix(G[j], N, N);
+ printf("GBar:\n");
+ printMatrix(GBar[j], N, N);
+ printf("Q=%d\n", Q[j]);
+ printf("D:\n");
+ printVector(D[j], K[j]);
+ printf("J:\n");
+ printMatrix(J[j], K[j], K[j]);*/
+ //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+ /*printf("\nafter:\n");
+ printf("K=%d\n", K[j]);
+ printf("h:\n");
+ printVector(h[j], N);
+ printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
+ printf("G:\n");
+ printMatrix(G[j], N, N);
+ printf("GBar:\n");
+ printMatrix(GBar[j], N, N);
+ printf("Q=%d\n", Q[j]);
+ printf("D:\n");
+ printVector(D[j], K[j]);
+ printf("J:\n");
+ printMatrix(J[j], K[j], K[j]);*/
+
+ //}
+
+ Paulicounter++;
+ }
+
+ amplitude = 0.0 + 0.0*I;
+ for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
+ //printf("i=%d\n", i);
+
+ randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
+
+ origGamma = 1.0 + 0.0*I;
+
+ for(k=0; k<Paulicounter; k++) {
+ origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
+ //printf("k=%d\n", k);
+ }
+ /*printf("origK=%d\n", origK);
+ printf("origG:\n");
+ printMatrix(origG, N, N);
+ printf("origGBar:\n");
+ printMatrix(origGBar, N, N);
+ printf("origh:\n");
+ printVector(origh, N);*/
+
+ double complex stabstateaverage = 0.0 + 0.0*I;
+
+ for(j=0; j<numStabStates; j++) {
+ //printf("j=%d\n", j);
+ 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);
+ stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
+ }
+ amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+
+ deallocate_mem(&origG, N);
+ deallocate_mem(&origGBar, N);
+ free(origh);
+ deallocate_mem(&origJ, origK);
+ free(origD);
+ }
+
+ //printf("numStabStates=%d\n", numStabStates);
+ printf("L1Norm=%lf\n", L1Norm);
+
+ printf("\namplitude:\n");
+ if(creal(amplitude+ZEROTHRESHOLD)>0)
+ printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+ else
+ printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+ printf("\nabs(amplitude):\n");
+ printf("%lf\n", cabs(amplitude));
+
+
+ for(i=0; i<PdN; i++)
+ free(Pd[i]);
+ free(Pd);
+
+ return 0;
+}
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int newomega, newalpha, newbeta, newgamma, newdelta;
+ int i;
+
+ if(scanf("%d", &newomega) != EOF) {
+ *omega = newomega;
+ for(i=0; i<numqubits; i++) {
+ if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
+ printf("Error: Too few input coeffs!\n");
+ exit(0);
+ }
+ if(newalpha+newbeta+newgamma+newdelta > 1) {
+ printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+ exit(0);
+ }
+ alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
+ }
+ return 1;
+ } else
+ return 0;
+
+}