/* When the user clicks on the button, toggle between hiding and showing the dropdown content */



added code, README and Makefile
[strong_simulation_gauss_sum_rank.git] / 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