strongsim1.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <complex.h>
   6 #include "matrix.h"
   7 #include "exponentialsum.h"
   8 #include "shrink.h"
   9 #include "extend.h"
  10 #include "measurepauli.h"
  11 #include "innerproduct.h"
  12
  13 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
  14
  15 // order of matrix elements is [row][column]!!!
  16
  17 int main()
  18 {
  19
  20   int N;              // number of qubits
  21   scanf("%d", &N);
  22
  23   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)
  24   scanf("%d", &T);
  25
  26   int omega;
  27   int alpha[N], beta[N], gamma[N], delta[N];
  28   int Paulicounter = 0;
  29
  30   int i, j, k, l;
  31
  32
  33   int n1 = 1; int k1 = 0; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {0}; int (*(J1[])) = { (int[]) {0} };
  34   int n2 = 1; int k2 = 0; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {1}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
  35
  36   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
  37   double complex Gamma[(int)pow(2,N)]; // prefactor in front of resultant state
  38   G = calloc(pow(2,N),sizeof(int*)); GBar = calloc(pow(2,N),sizeof(int*));
  39   h = calloc(pow(2,N),sizeof(int*));
  40
  41   J = calloc(pow(2,N),sizeof(int*)); D = calloc(pow(2,N),sizeof(int*)); Q = calloc(pow(2,N),sizeof(int));
  42
  43   K = calloc(pow(2,N), sizeof(int));
  44
  45   double complex origGamma[(int)pow(2,N)];
  46   int *origK, *origQ, **origD, ***origJ;
  47   int ***origG, ***origGBar, **origh;
  48
  49   origG = calloc(pow(2,N),sizeof(int*)); origGBar = calloc(pow(2,N),sizeof(int*));
  50   origh = calloc(pow(2,N),sizeof(int*));
  51
  52   origJ = calloc(pow(2,N),sizeof(int*)); origD = calloc(pow(2,N),sizeof(int*)); origQ = calloc(pow(2,N),sizeof(int));
  53
  54   origK = calloc(pow(2,N), sizeof(int));
  55
  56   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
  57
  58
  59   for(j=0; j<pow(2,N); j++) { // there will be 2^N combinations when using k=1 tensor factors
  60     combination = j;
  61
  62     Gamma[j] = 1.0;
  63
  64     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
  65     h[j] = calloc(N, sizeof(int));
  66
  67     if(K[j] > 0) {
  68
  69       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
  70     }
  71
  72     origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
  73     origh[j] = calloc(N, sizeof(int));
  74
  75     if(K[j] > 0) {
  76       origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
  77     }
  78
  79     for(k=0; k<N; k++) {
  80
  81       Gamma[j] *= (1.0/sqrt(2.0))*((combination & 0x1)*cexp(0.25*PI*I) + (~combination & 0x1)*cexp(0.0*PI*I));
  82
  83       G[j][k] = calloc(N, sizeof(int*)); GBar[j][k] = calloc(N, sizeof(int*));
  84
  85       if(K[j] > 0)
  86         J[j] = calloc(K[j], sizeof(int*));
  87
  88       origG[j][k] = calloc(N, sizeof(int*)); origGBar[j][k] = calloc(N, sizeof(int*));
  89
  90       if(K[j] > 0)
  91         origJ[j] = calloc(K[j], sizeof(int*));
  92
  93       h[j][k] = (combination & 0x1)*h2[0] + (~combination & 0x1)*h1[0];
  94       for(l=0; l<n1; l++) { // assuming n1=n2
  95         G[j][k][k*n1+l] = (combination & 0x1)*G2[0][l] + (~combination & 0x1)*G1[0][l];
  96         GBar[j][k][k*n1+l] = (combination & 0x1)*GBar2[0][l] + (~combination & 0x1)*GBar1[0][l];
  97       }
  98       memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
  99
 100       combination /= 2; // shift to the right by one (in base-2 arithmetic)
 101     }
 102
 103     memcpy(origh[j], h[j], N*sizeof(int));
 104
 105   }
 106
 107   memcpy(origGamma, Gamma, pow(2,N)*sizeof(double complex));
 108
 109   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
 110
 111     Paulicounter++;
 112     if(Paulicounter > N) {
 113       printf("Error: Number of Paulis is greater than N!\n");
 114       return 1;
 115     }
 116
 117     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 118     // Y_i = -I*Z*X
 119     for(i=0; i<N; i++) {
 120       if(delta[i]){
 121         omega += 3; // -I = I^3
 122         beta[i] = delta[i];
 123         gamma[i] = delta[i];
 124       }
 125     }
 126
 127     for(j=0; j<pow(2,N); j++) { // the kets
 128
 129       Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
 130
 131     }
 132
 133   }
 134
 135   double complex amplitude = 0.0 + 0.0*I;
 136   for(i=0; i<pow(2,N); i++) { // the bras
 137     for(j=0; j<pow(2,N); j++) {
 138       // check to see if second arguments are modified!!! They shouldn't be!
 139       double complex newamplitude = conj(origGamma[i])*Gamma[j]*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]);
 140       amplitude = amplitude + newamplitude;
 141     }
 142   }
 143
 144   printf("amplitude:\n");
 145   if(creal(amplitude+0.00000001)>0)
 146     printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
 147   else
 148     printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
 149
 150   return 0;
 151
 152 }
 153
 154
 155
 156 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 157 {
 158
 159   int newomega, newalpha, newbeta, newgamma, newdelta;
 160   int i;
 161
 162   if(scanf("%d", &newomega) != EOF) {
 163     *omega = newomega;
 164     for(i=0; i<numqubits; i++) {
 165       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
 166         printf("Error: Too few input coeffs!\n");
 167         exit(0);
 168       }
 169       if(newalpha+newbeta+newgamma+newdelta > 1) {
 170         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 171         exit(0);
 172       }
 173       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 174     }
 175     return 1;
 176   } else
 177     return 0;
 178
 179 }
 180
 181
 182
 183
 184
 185
 186