hilbertspace_vector.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* 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;
  28
  29   int N;              // number of qubits
  30   scanf("%d", &N);
  31
  32   complex double (*(psi[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
  33
  34   complex double **psiN = psi;
  35   for(i=1; i<N; i++)
  36     psiN = outerMatrix(psiN, psi, pow(2,i), 1, 2, 1);
  37
  38   int alpha[N], beta[N], gamma[N], delta[N];
  39
  40   complex double **P1[N];
  41
  42   complex double tr;
  43
  44   /* printf("psiN:\n"); */
  45   /* for(i=0; i<pow(2,N); i++) { */
  46   /*   printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0])); */
  47   /* } */
  48
  49   while(readPaulicoeffs(alpha, beta, gamma, delta, N)) {
  50
  51     for(i=0; i<N; i++) {
  52       //printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
  53       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);
  54     }
  55
  56     complex double **P = P1[0];
  57     for(i=1; i<N; i++)
  58       P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
  59
  60     complex double **psiNfinal = multMatrix(P,psiN,pow(2,N),pow(2,N),pow(2,N),1);
  61
  62     /* printf("psiNfinal:\n"); */
  63     /* for(i=0; i<pow(2,N); i++) { */
  64     /*   printf("%d: %lf+%lfI\n", i, creal(psiNfinal[i][0]), cimag(psiNfinal[i][0])); */
  65     /* } */
  66
  67     tr = 0.0 + 0.0*I;
  68     //printf("tr:\n");
  69     for(i=0; i<pow(2,N); i++) {
  70       tr += conj(psiN[i][0])*psiNfinal[i][0];
  71       //printf("%d: %lf+%lfI\n", i, creal(tr), cimag(tr));
  72     }
  73
  74     //    printf("%lf+%lfI\n", cabs(creal(tr)), cimag(tr));
  75     printf("%lf\n", cabs(creal(tr)));
  76     /* 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 */
  77
  78     deallocate_mem(&P, pow(2,N));
  79     deallocate_mem(&psiNfinal, pow(2,N));
  80
  81   }
  82
  83   //  deallocate_mem(&psiN, pow(2,N));
  84
  85   return 0;
  86 }
  87
  88
  89 complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
  90 {
  91   int i, j;
  92
  93   complex double** C;
  94
  95   C = calloc(cols, sizeof(complex double*));
  96   for(i=0; i<cols; i++)
  97     C[i] = calloc(rows, sizeof(complex double));
  98
  99   for(i=0; i<rows; i++)
 100     for(j=0; j<cols; j++)
 101       C[i][j] = A[i][j] + B[i][j];
 102
 103   return C;
 104 }
 105
 106 complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
 107 {
 108   int i, j;
 109
 110   complex double** C;
 111
 112   C = calloc(cols, sizeof(complex double*));
 113   for(i=0; i<cols; i++)
 114     C[i] = calloc(rows, sizeof(complex double));
 115
 116   for(i=0; i<rows; i++)
 117     for(j=0; j<cols; j++)
 118       C[i][j] = scalar*a[i][j];
 119
 120   return C;
 121 }
 122
 123 complex double trace(complex double **a, int rows, int cols)
 124 {
 125   int i;
 126   complex double tr = 0.0*I;
 127
 128   for(i=0; i<rows; i++)
 129     tr += a[i][i];
 130
 131   return tr;
 132 }
 133
 134 complex double** transp(complex double **a, int rows, int cols)
 135 {
 136   int i, j;
 137   complex double** C;
 138
 139   C = calloc(cols, sizeof(complex double*));
 140   for(i=0; i<cols; i++)
 141     C[i] = calloc(rows, sizeof(complex double));
 142
 143   for(i=0; i<cols; i++)
 144     for(j=0; j<rows; j++) {
 145       C[i][j] = a[j][i];
 146     }
 147
 148   return C;
 149 }
 150
 151 complex double** conjtransp(complex double **a, int rows, int cols)
 152 {
 153   int i, j;
 154   complex double** C;
 155
 156   C = calloc(cols, sizeof(complex double*));
 157   for(i=0; i<cols; i++)
 158     C[i] = calloc(rows, sizeof(complex double));
 159
 160   for(i=0; i<cols; i++)
 161     for(j=0; j<rows; j++) {
 162       C[i][j] = conj(a[j][i]);
 163     }
 164
 165   return C;
 166 }
 167
 168 void printMatrix(complex double** a, int rows, int cols)
 169 {
 170   int i, j;
 171   printf("Matrix[%d][%d]\n", rows, cols);
 172   for(i=0; i<rows; i++) {
 173     for(j=0; j<cols; j++) {
 174       printf("%lf+%lfI ", creal(a[i][j]), cimag(a[i][j]));
 175     }
 176     printf("\n");
 177   }
 178 }
 179
 180 complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
 181 {
 182   int i, j, k;
 183   complex double **C;
 184   C = calloc(ro1, sizeof(complex double*));
 185   for(i=0; i<ro1; i++)
 186     C[i] = calloc(co2, sizeof(complex double));
 187
 188   for(i=0; i<ro1; i++) {
 189     for(j=0; j<co2; j++) {
 190       C[i][j] = 0;
 191       for(k=0; k<co1; k++)
 192         C[i][j] += A[i][k] * B[k][j];
 193     }
 194   }
 195
 196   return C;
 197 }
 198
 199 complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
 200 {
 201   int i, j, k, l;
 202   complex double **C;
 203   C = calloc(ro1*ro2, sizeof(complex double*));
 204   for(i=0; i<ro1*ro2; i++)
 205     C[i] = calloc(co1*co2, sizeof(complex double));
 206
 207   for(i=0; i<ro1; i++)
 208     for(j=0; j<ro2; j++)
 209       for(k=0; k<co1; k++)
 210         for(l=0; l<co2; l++) {
 211           C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
 212         }
 213
 214   return C;
 215 }
 216
 217
 218
 219 void deallocate_mem(complex double ***arr, int rows)
 220 {
 221   int i;
 222   for(i=0; i<rows; i++)
 223     free((*arr)[i]);
 224   free(*arr);
 225 }
 226
 227 int readPaulicoeffs(int *alpha, int *beta, int *gamma, int *delta, int numqubits)
 228 {
 229
 230   int i;
 231
 232   if(scanf("%d %d %d %d", &alpha[0], &beta[0], &gamma[0], &delta[0]) != EOF) {
 233     for(i=1; i<numqubits; i++) {
 234       scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]);
 235     }
 236     return 1;
 237   } else
 238     return 0;
 239
 240 }
 241
 242
 243