multipauli.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <complex.h>
   4 #include <math.h>
   5
   6 void deallocate_mem(complex double ***arr, int rows);
   7 void printMatrix(complex double** a, int rows, int cols);
   8 complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
   9 complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
  10 complex double** transp(complex double **a, int rows, int cols);
  11 complex double** conjtransp(complex double **a, int rows, int cols);
  12 complex double trace(complex double **a, int rows, int cols);
  13 complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols);
  14 complex double** addMatrix(complex double **A, complex double **B, int rows, int cols);
  15 int readPaulicoeffs(int* omega, int* alpha, int* beta, int* gamma, int* delta, int numqubits);
  16
  17 // order of matrix elements is [row][column]!!!
  18
  19 static complex double (*(I2[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, 1.0+0.0*I} };
  20 static complex double (*(X[])) = { (double complex[]) {0.0*I, 1.0+0.0*I}, (double complex[]) {1.0+0.0*I, 0.0*I} };
  21 static complex double (*(Y[])) = { (double complex[]) {0.0*I, 0.0-1.0*I}, (double complex[]) {0.0+1.0*I, 0.0*I} };
  22 static complex double (*(Z[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, -1.0+0.0*I} };
  23
  24 int main()
  25 {
  26
  27   int i, j;
  28
  29   int N;              // number of qubits
  30   scanf("%d", &N);
  31
  32   int K;              // 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)
  33   scanf("%d", &K);
  34
  35   complex double **IN;  // N-qubit identity matrix
  36   IN = I2;
  37   for(i=1; i<N; i++) {
  38     IN = outerMatrix(IN,I2,pow(2,i),pow(2,i),2,2);
  39   }
  40
  41   complex double (*(psiT[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
  42   complex double (*(psi0[])) = { (double complex[]) {1.0+0.0*I}, (double complex[]) {0.0*I}}; // T gate magic state
  43
  44   complex double **psiN;
  45   if(K > 0) {
  46     psiN = psiT;
  47     for(i=1; i<K; i++)
  48       psiN = outerMatrix(psiN, psiT, pow(2,i), 1, 2, 1);
  49     for(i=K; i<N; i++)
  50       psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
  51   } else {
  52     psiN = psi0;
  53     for(i=1; i<N; i++)
  54       psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
  55   }
  56
  57
  58   int omega, alpha[N], beta[N], gamma[N], delta[N];
  59
  60   int Paulicounter = 0;
  61
  62   complex double **fullP;  // full product: \prod_i 1/2*(1+P_i)
  63   complex double **P;      // P (P_i above) is made up of products of one-qubit Paulis, P1
  64   complex double **P1[N];  // one-qubit Paulis
  65
  66   complex double tr;
  67
  68   printf("psiN:\n");
  69   for(i=0; i<pow(2,N); i++) {
  70     printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0]));
  71   }
  72
  73
  74   while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) { // go over the product of 1/2*(I+Paulis) that makes up the full projector
  75
  76     Paulicounter++;
  77     if(Paulicounter > N) {
  78       printf("Error: Number of Paulis is greater than N!\n");
  79       return 1;
  80     }
  81
  82
  83     printf("%d\n", omega);
  84     for(i=0; i<N; i++) {
  85       printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
  86       P1[i] = addMatrix(addMatrix(addMatrix(scalarmultMatrix(alpha[i],I2,2,2),scalarmultMatrix(beta[i],Z,2,2),2,2),scalarmultMatrix(gamma[i],X,2,2),2,2),scalarmultMatrix(delta[i],Y,2,2),2,2);
  87     }
  88
  89     P = P1[0];
  90     for(i=1; i<N; i++)
  91       P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
  92     P = scalarmultMatrix(cpow(I,omega),P,pow(2,N),pow(2,N));
  93
  94     if(Paulicounter == 1)
  95       fullP = scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N));
  96     else {
  97       fullP = multMatrix(scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N)),fullP,pow(2,N),pow(2,N),pow(2,N),pow(2,N));
  98     }
  99     deallocate_mem(&P, pow(2,N));
 100
 101   }
 102
 103   complex double **psiNfinal = multMatrix(fullP,psiN,pow(2,N),pow(2,N),pow(2,N),1);
 104
 105   printf("fullP:\n");
 106   for(i=0; i<pow(2,N); i++) {
 107     for(j=0; j<pow(2,N); j++) {
 108       printf("%lf+%lfI ", creal(fullP[i][j]), cimag(fullP[i][j]));
 109     }
 110     printf("\n");
 111   }
 112   printf("psiNfinal:\n");
 113   for(i=0; i<pow(2,N); i++) {
 114     printf("%d: %lf+%lfI\n", i, creal(psiNfinal[i][0]), cimag(psiNfinal[i][0]));
 115   }
 116
 117   tr = 0.0 + 0.0*I;
 118   //printf("tr:\n");
 119   for(i=0; i<pow(2,N); i++) {
 120     tr += conj(psiN[i][0])*psiNfinal[i][0];
 121     //printf("%d: %lf+%lfI\n", i, creal(tr), cimag(tr));
 122   }
 123
 124   if(creal(tr+0.00000001)>0)
 125     printf("%.10lf %c %.10lf I\n", cabs(creal(tr)), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
 126   else
 127     printf("%.10lf %c %.10lf I\n", creal(tr), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
 128   //printf("%lf %c %lf I\n", creal(tr), cimag(tr)>0?'+':'-' , cabs(cimag(tr)));
 129   //printf("%lf\n", cabs(creal(tr)));
 130   /* cabs the creal part because it's always positive, but sometimes the 0.0 gets a minus sign which is annoying to see when comparing outputs */
 131
 132   deallocate_mem(&psiNfinal, pow(2,N));
 133
 134
 135   //  deallocate_mem(&psiN, pow(2,N));
 136
 137   return 0;
 138 }
 139
 140
 141 complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
 142 {
 143   int i, j;
 144
 145   complex double** C;
 146
 147   C = calloc(cols, sizeof(complex double*));
 148   for(i=0; i<cols; i++)
 149     C[i] = calloc(rows, sizeof(complex double));
 150
 151   for(i=0; i<rows; i++)
 152     for(j=0; j<cols; j++)
 153       C[i][j] = A[i][j] + B[i][j];
 154
 155   return C;
 156 }
 157
 158 complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
 159 {
 160   int i, j;
 161
 162   complex double** C;
 163
 164   C = calloc(cols, sizeof(complex double*));
 165   for(i=0; i<cols; i++)
 166     C[i] = calloc(rows, sizeof(complex double));
 167
 168   for(i=0; i<rows; i++)
 169     for(j=0; j<cols; j++)
 170       C[i][j] = scalar*a[i][j];
 171
 172   return C;
 173 }
 174
 175 complex double trace(complex double **a, int rows, int cols)
 176 {
 177   int i;
 178   complex double tr = 0.0*I;
 179
 180   for(i=0; i<rows; i++)
 181     tr += a[i][i];
 182
 183   return tr;
 184 }
 185
 186 complex double** transp(complex double **a, int rows, int cols)
 187 {
 188   int i, j;
 189   complex double** C;
 190
 191   C = calloc(cols, sizeof(complex double*));
 192   for(i=0; i<cols; i++)
 193     C[i] = calloc(rows, sizeof(complex double));
 194
 195   for(i=0; i<cols; i++)
 196     for(j=0; j<rows; j++) {
 197       C[i][j] = a[j][i];
 198     }
 199
 200   return C;
 201 }
 202
 203 complex double** conjtransp(complex double **a, int rows, int cols)
 204 {
 205   int i, j;
 206   complex double** C;
 207
 208   C = calloc(cols, sizeof(complex double*));
 209   for(i=0; i<cols; i++)
 210     C[i] = calloc(rows, sizeof(complex double));
 211
 212   for(i=0; i<cols; i++)
 213     for(j=0; j<rows; j++) {
 214       C[i][j] = conj(a[j][i]);
 215     }
 216
 217   return C;
 218 }
 219
 220 void printMatrix(complex double** a, int rows, int cols)
 221 {
 222   int i, j;
 223   printf("Matrix[%d][%d]\n", rows, cols);
 224   for(i=0; i<rows; i++) {
 225     for(j=0; j<cols; j++) {
 226       printf("%lf+%lfI ", creal(a[i][j]), cimag(a[i][j]));
 227     }
 228     printf("\n");
 229   }
 230 }
 231
 232 complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
 233 {
 234   int i, j, k;
 235   complex double **C;
 236   C = calloc(ro1, sizeof(complex double*));
 237   for(i=0; i<ro1; i++)
 238     C[i] = calloc(co2, sizeof(complex double));
 239
 240   for(i=0; i<ro1; i++) {
 241     for(j=0; j<co2; j++) {
 242       C[i][j] = 0;
 243       for(k=0; k<co1; k++)
 244         C[i][j] += A[i][k] * B[k][j];
 245     }
 246   }
 247
 248   return C;
 249 }
 250
 251 complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
 252 {
 253   int i, j, k, l;
 254   complex double **C;
 255   C = calloc(ro1*ro2, sizeof(complex double*));
 256   for(i=0; i<ro1*ro2; i++)
 257     C[i] = calloc(co1*co2, sizeof(complex double));
 258
 259   for(i=0; i<ro1; i++)
 260     for(j=0; j<ro2; j++)
 261       for(k=0; k<co1; k++)
 262         for(l=0; l<co2; l++)
 263           C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
 264
 265   return C;
 266 }
 267
 268
 269
 270 void deallocate_mem(complex double ***arr, int rows)
 271 {
 272   int i;
 273   for(i=0; i<rows; i++)
 274     free((*arr)[i]);
 275   free(*arr);
 276 }
 277
 278 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 279 {
 280
 281   int i;
 282
 283   if(scanf("%d", omega) != EOF) {
 284     for(i=0; i<numqubits; i++) {
 285       if(scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]) == EOF) {
 286         printf("Error: Too few input coeffs!\n");
 287         exit(0);
 288       }
 289       if(alpha[i]+beta[i]+gamma[i]+delta[i] > 1) {
 290         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
 291         exit(0);
 292       }
 293     }
 294     return 1;
 295   } else
 296     return 0;
 297
 298 }
 299
 300
 301