#include "measurepauli.h"
#include "innerproduct.h"
#include "sparsify.h"
+#include "omp.h"
int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
int main(int argc, char *argv[])
{
- if(argc != 4) {
- printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+ if(argc != 5 && argc != 6) {
+ printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\" \"seed\" \"(optional) sample number prefactor\"\n");
exit(0);
}
double additiveErrorDelta = atof(argv[1]); // additive error delta
double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
- int coherentSampling = atoi(argv[3]); // perform coherent sampling (true=1 or false=0)
+ int coherentSampling = atoi(argv[3]); // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
+ int seed = atoi(argv[4]); // random number seed
+ double samplePrefactor = 1.0;
+ if(argc == 6)
+ samplePrefactor = atof(argv[5]); // factor that multiplies k, the number of samples used
int N; // number of qubits
scanf("%d", &N);
long* stabStateIndices;
int numStabStates;
- srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+ //srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+ srand((unsigned)seed); // seeding the random number generator for sparsify()
- if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+ if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling, samplePrefactor))
return 1;
- //printf("checking: numStabStateIndices:\n");
- //for(i=0; i<numStabStates; i++)
- // printf("%d ", stabStateIndices[i]);
- //printf("\n");
- //fflush(stdout);
+ /* 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
memcpy(origQ, Q, numStabStates*sizeof(int));
+ double complex amplitude = 0.0 + 0.0*I;
+ double complex lastamplitude = 0.0 + 0.0*I;
+
+ double complex probability = 1.0 + 0.0*I;
+
+ // the first measurement should be the identity
+ omega = 0;
+ for(i=0; i<N; i++) {
+ alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+ }
+
+ for(j=0; j<numStabStates; j++) { // the kets
+ Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+ }
+
+ double complex newamplitude;
+
+ #pragma omp parallel
+ {
+ #pragma omp for schedule(static) collapse(2)
+ for(i=0; i<numStabStates; i++) { // the bras
+ for(j=0; j<numStabStates; j++) {
+ newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+ amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+ }
+ }
+ }
+
while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
Paulicounter++;
}
- }
-
- double complex amplitude = 0.0 + 0.0*I;
- for(i=0; i<numStabStates; i++) { // the bras
- for(j=0; j<numStabStates; j++) {
- double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
- amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+ lastamplitude = amplitude;
+ amplitude = 0.0 + 0.0*I;
+ #pragma omp parallel
+ {
+ #pragma omp for schedule(static) collapse(2)
+ for(i=0; i<numStabStates; i++) { // the bras
+ for(j=0; j<numStabStates; j++) {
+ newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+ amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+ }
+ }
}
+ printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
+ printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+
+ probability *= amplitude/lastamplitude;
+ //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
+
+ printf("probability = %lf\n", cabs(probability));
+ printf("!\n");
}
//printf("numStabStates=%d\n", numStabStates);
printf("\nabs(amplitude):\n");
printf("%lf\n", cabs(amplitude));
+ printf("probability = %lf\n", cabs(probability));
+
return 0;
}
#include <time.h>
#include "matrix.h"
#include "exponentialsum.h"
-#include "shrinkstar.h"
+#include "shrink.h"
#include "extend.h"
#include "measurepauli.h"
#include "innerproduct.h"
{
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");
+ printf("weaksim_relerr argument: \"number of stabilizer state samples\" \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\\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 coherentSampling = atoi(argv[4]); // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3)
- int N; // number of qubits
+ 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);
+ 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 omega;
+ int alpha[N], beta[N], gamma[N], delta[N];
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()
+ srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate _equatorial()
fp = fopen("Pd.txt", "r");
int PdN = atoi(buff);
Pd = calloc(PdN, sizeof(double*));
- for(i=0; i<PdN; i++)
+ for(i=0; i<PdN; i++)
Pd[i] = calloc(PdN+1, sizeof(double));
double tmp;
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;
+ 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]);
//printf("total=%f\n", tmp);
}
- double complex amplitude;
-
+ //|T> = e^(pi*i/8)/2*sqrt(4-2*sqrt(2))* (sqrt(-i)*(|0>+i|1>)/sqrt(2) + (|0>+|1>)/sqrt(2))
+
+ //double additiveErrorDelta = 0.1;
+
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
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);
+ /* 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
K = calloc(numStabStates, sizeof(int));
- int origK, origQ, *origD;
- int **origJ;
- int **origG, **origGBar;
- int *origh;
- double complex origGamma;
+ double complex origGamma[(int)NUMSTABSTATESAMPLES];
+ int *origK, *origQ, **origD, ***origJ;
+ int ***origG, ***origGBar, **origh;
+
+ origG = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origGBar = calloc(NUMSTABSTATESAMPLES,sizeof(int*));
+ origh = calloc(NUMSTABSTATESAMPLES,sizeof(int*));
+
+ origJ = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origD = calloc(NUMSTABSTATESAMPLES,sizeof(int*)); origQ = calloc(NUMSTABSTATESAMPLES,sizeof(int));
+
+ origK = calloc(NUMSTABSTATESAMPLES, sizeof(int));
- long combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+ int 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);
+ printf("!\n");
+ fflush(stdout);
+
+ for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+ origGamma[j] = 1.0+0.0*I;
+ randomstabilizerstate(N, &origK[j], &origh[j], &origG[j], &origGBar[j], &origQ[j], &origD[j], &origJ[j], Pd);
+ //randomstabilizerstate_equatorial(N, &origK[j], &origh[j], &origG[j], &origGBar[j], &origQ[j], &origD[j], &origJ[j]);
+ }
+
+ printf("NUMSTABSTATESAMPLES=%d stabilizer states allocated and initialized!\n", NUMSTABSTATESAMPLES);
+ fflush(stdout);
+
for(j=0; j<numStabStates; j++) {
combination = stabStateIndices[j];
K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
combination /= 2;
}
- combination = j;
+ combination = stabStateIndices[j];
+ //origK[j] = K[j];
Gamma[j] = 1.0;
Gamma[j] *= L1Norm/((double)numStabStates);
J[j][k] = calloc(K[j], sizeof(int));
}
+ //origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
+ //origh[j] = calloc(N, sizeof(int));
+
+ //if(K[j] > 0) {
+ // origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
+ // for(k=0; k<K[j]; k++)
+ // origJ[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));
+ //origG[j][k] = calloc(N, sizeof(int)); origGBar[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'
for(k=0; k<N; k++) {
- Q[j] += (((combination%2)==1)*Q2 + ((combination%2)==0)*Q1);
+ 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;
+ // D1 may have 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
+ // J1 may have 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];
}
for(l=0; l<n1; l++) { // assuming n1=n2
- h[j][k*n1+l] = (((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l]);
+ 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]);
+ 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;
//memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
- combination /= 2; // shift to the right by one (in base-7 arithmetic)
+ combination /= 2; // shift to the right by one (in base-2 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;
+ combination = stabStateIndices[j];
for(k=0; k<(N); k++) {
Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
//printf("Kcounter=%d\n", Kcounter);
/* 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]);
+ 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);
combination /= 2;
}
+ /* for(k=0; k<N; k++) { */
+ /* memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int)); */
+ /* } */
+ /* for(k=0; k<K[j]; k++) { */
+ /* memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int)); */
+ /* } */
+
+ /* memcpy(origh[j], h[j], N*sizeof(int)); */
+ /* memcpy(origD[j], D[j], K[j]*sizeof(int)); */
- /*printf("G[%d]:\n", j);
- printMatrix(G[j], N, N);
- printf("GBar[%d]:\n", j);
- printMatrix(GBar[j], N, N);
+ }
+ //exit(0);
+ /* memcpy(origGamma, Gamma, numStabStates*sizeof(double complex)); */
- printf("h[%d]:\n", j);
- printVector(h[j], N);
+ /* memcpy(origQ, Q, numStabStates*sizeof(int)); */
- printf("J[%d]:\n", j);
- printMatrix(J[j], K[j], K[j]);
-
- printf("D[%d]:\n", j);
- printVector(D[j], K[j]);
+ printf("numStabStates=%d stabilizer states allocated and initialized!\n", numStabStates);
+ fflush(stdout);
+
+ double complex amplitude = 0.0 + 0.0*I;
+ double complex lastamplitude = 0.0 + 0.0*I;
- printf("Q[%d]=%d\n", j, Q[j]);*/
+ double complex probability = 1.0 + 0.0*I;
+ // the first measurement should be the identity
+ omega = 0;
+ for(i=0; i<N; i++) {
+ alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+ }
+
+ for(j=0; j<numStabStates; j++) { // the kets
+ Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
}
- //exit(0);
- while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
+ double complex newamplitude;
+ double complex stabstateaverage;
+ for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+ stabstateaverage = 0.0 + 0.0*I;
+ for(i=0; i<numStabStates; i++) { // the bras
+ //newamplitude = innerproduct_equatorial(N, K[i], h[i], G[i], GBar[i], Q[i], D[i], J[i], N, origK[j], origh[j], origG[j], origGBar[j], origQ[j], origD[j], origJ[j]);
+ newamplitude = innerproduct(N, K[i], h[i], G[i], GBar[i], Q[i], D[i], J[i], N, origK[j], origh[j], origG[j], origGBar[j], origQ[j], origD[j], origJ[j]);
+ stabstateaverage = stabstateaverage + origGamma[j]*Gamma[i]*newamplitude;
+ }
+ amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+ }
- if((Paulicounter+1) > N) {
+ while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
+
+ Paulicounter++;
+ if(Paulicounter > 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];
+ if(delta[i]){
+ omega += 3; // -I = I^3
+ beta[i] = delta[i];
+ gamma[i] = delta[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]);*/
- //}
+ for(j=0; j<numStabStates; j++) { // the kets
- Paulicounter++;
- }
-
- amplitude = 0.0 + 0.0*I;
- for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
- //printf("i=%d\n", i);
+ Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
- 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;
+ lastamplitude = amplitude;
+ amplitude = 0.0 + 0.0*I;
+ for(j=0; j<NUMSTABSTATESAMPLES; j++) {
+ stabstateaverage = 0.0 + 0.0*I;
+ for(i=0; i<numStabStates; i++) { // the bras
+ //newamplitude = innerproduct_equatorial(N, K[i], h[i], G[i], GBar[i], Q[i], D[i], J[i], N, origK[j], origh[j], origG[j], origGBar[j], origQ[j], origD[j], origJ[j]);
+ newamplitude = innerproduct(N, K[i], h[i], G[i], GBar[i], Q[i], D[i], J[i], N, origK[j], origh[j], origG[j], origGBar[j], origQ[j], origD[j], origJ[j]);
+ stabstateaverage = stabstateaverage + origGamma[j]*Gamma[i]*newamplitude;
+ }
+ amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
}
- amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+ printf("lastamplitude: %lf %c %lf I\n", cabs(creal(lastamplitude)), cimag(lastamplitude+0.00000001)>0?'+':'-' , cabs(cimag(lastamplitude)));
+ printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+
+ probability *= amplitude/lastamplitude;
+ //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement
- deallocate_mem(&origG, N);
- deallocate_mem(&origGBar, N);
- free(origh);
- deallocate_mem(&origJ, origK);
- free(origD);
+ printf("probability = %lf\n", cabs(probability));
+ printf("!\n");
}
//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)));
+ if(creal(amplitude+0.00000001)>0)
+ printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
else
- printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+ printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+ //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
printf("\nabs(amplitude):\n");
printf("%lf\n", cabs(amplitude));
-
- for(i=0; i<PdN; i++)
- free(Pd[i]);
- free(Pd);
+ printf("probability = %lf\n", cabs(probability));
return 0;
+
}
+
+
int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
{
projects / weak_simulation_stab_extent.git / commitdiff