strongsim3.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   if(N%3 != 0) {
  24     printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n");
  25     return 1;
  26   }
  27
  28   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)
  29   scanf("%d", &T);
  30
  31   int omega;
  32   int alpha[N], beta[N], gamma[N], delta[N];
  33   int Paulicounter = 0;
  34
  35   int i, j, k, l, m;
  36
  37   double complex coeffa = -0.25*(1.0-I)*(-1.0-I+sqrt(2.0))*csqrt(-I);
  38   double complex coeffb = 0.25*(-1.0-I)*(1.0-I+sqrt(2.0))*csqrt(I);
  39   double complex coeffc = 0.25*(-1.0-I)*(-1.0+I+sqrt(2.0))*csqrt(I);
  40
  41   int n1 = 3; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1, 1}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0}, (int[]) {1, 1, 0}, (int[]) {1, 0, 1}}; int h1[] = {1, 1, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
  42   int n2 = 3; int k2 = 3; int (*(G2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h2[] = {0, 0, 0}; int Q2 = 2; int D2[] = {2, 2, 0}; int (*(J2[])) = { (int[]) {4, 0, 0}, (int[]) {0, 4, 0}, (int[]) {0, 0, 0} };
  43   int n3 = 3; int k3 = 3; int (*(G3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h3[] = {0, 0, 0}; int Q3 = 2; int D3[] = {6, 6, 0}; int (*(J3[])) = { (int[]) {4, 4, 4}, (int[]) {4, 4, 4}, (int[]) {4, 4, 0} };
  44
  45   int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
  46   double complex Gamma[(int)pow(3,N/3)]; // prefactor in front of resultant state
  47   G = calloc(pow(3,N/3),sizeof(int*)); GBar = calloc(pow(3,N/3),sizeof(int*));
  48   h = calloc(pow(3,N/3),sizeof(int*));
  49
  50   J = calloc(pow(3,N/3),sizeof(int*)); D = calloc(pow(3,N/3),sizeof(int*)); Q = calloc(pow(3,N/3),sizeof(int));
  51
  52   K = calloc(pow(3,N/3), sizeof(int));
  53
  54   double complex origGamma[(int)pow(3,N/3)];
  55   int *origK, *origQ, **origD, ***origJ;
  56   int ***origG, ***origGBar, **origh;
  57
  58   origG = calloc(pow(3,N/3),sizeof(int*)); origGBar = calloc(pow(3,N/3),sizeof(int*));
  59   origh = calloc(pow(3,N/3),sizeof(int*));
  60
  61   origJ = calloc(pow(3,N/3),sizeof(int*)); origD = calloc(pow(3,N/3),sizeof(int*)); origQ = calloc(pow(3,N/3),sizeof(int));
  62
  63   origK = calloc(pow(3,N/3), sizeof(int));
  64
  65   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
  66
  67
  68   for(j=0; j<pow(3,N/3); j++) { // there will be 3^(N/3) combinations when using k=3 tensor factors
  69
  70     combination = j;
  71
  72     K[j] = 0.0;
  73
  74     for(k=0; k<N/3; k++) {
  75       K[j] += (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
  76       combination /= 3;
  77     }
  78     combination = j;
  79     origK[j] = K[j];
  80
  81     Gamma[j] = 1.0;
  82
  83     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
  84     h[j] = calloc(N, sizeof(int));
  85
  86     if(K[j] > 0) {
  87       J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
  88       for(k=0; k<K[j]; k++)
  89         J[j][k] = calloc(K[j], sizeof(int));
  90     }
  91
  92     origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
  93     origh[j] = calloc(N, sizeof(int));
  94
  95     if(K[j] > 0) {
  96       origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
  97       for(k=0; k<K[j]; k++)
  98         origJ[j][k] = calloc(K[j], sizeof(int));
  99     }
 100
 101     for(k=0; k<N; k++) {
 102       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
 103       origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
 104     }
 105
 106     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'
 107     int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
 108     for(k=0; k<N/3; k++) {
 109
 110       Q[j] += ((combination%3)==2)*Q3 + ((combination%3)==1)*Q2 + ((combination%3)==0)*Q1;
 111
 112       Gamma[j] *= (((combination%3)==2)*coeffc + ((combination%3)==1)*coeffb + ((combination%3)==0)*coeffa);
 113
 114       Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
 115       for(l=0; l<Kcombo; l++) {
 116           // 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
 117           switch(combination%3) {
 118           case 0:
 119             D[j][Kcounter+l] = D1[l];
 120             break;
 121           case 1:
 122             D[j][Kcounter+l] = D2[l];
 123             break;
 124           case 2:
 125             D[j][Kcounter+l] = D3[l];
 126             break;
 127           default:
 128             printf("error");
 129             return 1;
 130           }
 131         for(m=0; m<Kcombo; m++) {
 132           // 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
 133           switch(combination%3) {
 134           case 0:
 135             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
 136             break;
 137           case 1:
 138             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
 139             break;
 140           case 2:
 141             J[j][Kcounter+l][Kcounter+m] = J3[l][m];
 142             break;
 143           default:
 144             printf("error");
 145             return 1;
 146           }
 147         }
 148       }
 149
 150       for(l=0; l<n1; l++) { // assuming n1=n2=n3
 151         h[j][k*n1+l] = ((combination%3)==2)*h3[l] + ((combination%3)==1)*h2[l] + ((combination%3)==0)*h1[l];
 152       }
 153       // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
 154       for(l=0; l<Kcombo; l++) {
 155         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 156           G[j][Kcounter+l][k*n1+m] = ((combination%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m];
 157           GBar[j][Kcounter+l][k*n1+m] = ((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m];
 158         }
 159       }
 160       Kcounter = Kcounter + Kcombo;
 161
 162       /* printf("intermediate G[%d]:\n", j); */
 163       /* printMatrix(G[j], N, N); */
 164       /* printf("intermediate GBar[%d]:\n", j); */
 165       /* printMatrix(GBar[j], N, N); */
 166       //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 167
 168       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 169
 170       combination /= 3; // shift to the right by one (in base-3 arithmetic)
 171     }
 172     //printf("!\n");
 173
 174     // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
 175     combination = j;
 176     for(k=0; k<(N/3); k++) {
 177       Kcombo = (((combination%3)==2)*k3 + ((combination%3)==1)*k2 + ((combination%3)==0)*k1);
 178       //printf("Kcounter=%d\n", Kcounter);
 179       // G and GBar rows that are outside the first 'k' spanning basis states
 180       for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
 181         //printf("l=%d\n", l);
 182         for(m=0; m<n1; m++) { // assuming n1=n2=n3
 183           /* printf("m=%d\n", m); */
 184           /* printf("Kcounter+l=%d\n", Kcounter+l); */
 185           /* printf("k*n1+m=%d\n", k*n1+m); */
 186           G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%3)==2)*G3[l][m] + ((combination%3)==1)*G2[l][m] + ((combination%3)==0)*G1[l][m];
 187           GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%3)==2)*GBar3[l][m] + ((combination%3)==1)*GBar2[l][m] + ((combination%3)==0)*GBar1[l][m];
 188         }
 189       }
 190       Kcounter = Kcounter + (n1-Kcombo);
 191
 192       /* printf("intermediate G[%d]:\n", j); */
 193       /* printMatrix(G[j], N, N); */
 194       /* printf("intermediate GBar[%d]:\n", j); */
 195       /* printMatrix(GBar[j], N, N); */
 196
 197       combination /= 3;
 198     }
 199     for(k=0; k<N; k++) {
 200       memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
 201     }
 202     for(k=0; k<K[j]; k++) {
 203       memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
 204     }
 205
 206     memcpy(origh[j], h[j], N*sizeof(int));
 207     memcpy(origD[j], D[j], K[j]*sizeof(int));
 208
 209   }
 210   //exit(0);
 211   memcpy(origGamma, Gamma, pow(3,N/3)*sizeof(double complex));
 212
 213   memcpy(origQ, Q, pow(3,N/3)*sizeof(int));
 214
 215   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
 216
 217     Paulicounter++;
 218     if(Paulicounter > N) {
 219       printf("Error: Number of Paulis is greater than N!\n");
 220       return 1;
 221     }
 222
 223     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
 224     // Y_i = -I*Z*X
 225     for(i=0; i<N; i++) {
 226       if(delta[i]){
 227         omega += 3; // -I = I^3
 228         beta[i] = delta[i];
 229         gamma[i] = delta[i];
 230       }
 231     }
 232
 233
 234     for(j=0; j<pow(3,N/3); j++) { // the kets
 235
 236       Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
 237
 238     }
 239
 240   }
 241
 242   double complex amplitude = 0.0 + 0.0*I;
 243   for(i=0; i<pow(3,N/3); i++) { // the bras
 244     for(j=0; j<pow(3,N/3); j++) {
 245       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]);
 246       amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
 247     }
 248   }
 249
 250   printf("amplitude:\n");
 251   if(creal(amplitude+0.00000001)>0)
 252     printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
 253   else
 254     printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
 255   //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
 256
 257   return 0;
 258
 259 }
 260
 261
 262
 263 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 264 {
 265
 266   int newomega, newalpha, newbeta, newgamma, newdelta;
 267   int i;
 268
 269   if(scanf("%d", &newomega) != EOF) {
 270     *omega = newomega;
 271     for(i=0; i<numqubits; i++) {
 272       if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
 273         printf("Error: Too few input coeffs!\n");
 274         exit(0);
 275       }
 276       if(newalpha+newbeta+newgamma+newdelta > 1) {
 277         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 278         exit(0);
 279       }
 280       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
 281     }
 282     return 1;
 283   } else
 284     return 0;
 285
 286 }