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);
17 // order of matrix elements is [row][column]!!!
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} };
29 int N; // number of qubits
32 complex double (*(psi[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
34 complex double **psiN = psi;
36 psiN = outerMatrix(psiN, psi, pow(2,i), 1, 2, 1);
38 int alpha[N], beta[N], gamma[N], delta[N];
40 complex double **P1[N];
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])); */
49 while(readPaulicoeffs(alpha, beta, gamma, delta, N)) {
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);
56 complex double **P = P1[0];
58 P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
60 complex double **psiNfinal = multMatrix(P,psiN,pow(2,N),pow(2,N),pow(2,N),1);
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])); */
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));
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 */
78 deallocate_mem(&P, pow(2,N));
79 deallocate_mem(&psiNfinal, pow(2,N));
83 // deallocate_mem(&psiN, pow(2,N));
89 complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
95 C = calloc(cols, sizeof(complex double*));
97 C[i] = calloc(rows, sizeof(complex double));
100 for(j=0; j<cols; j++)
101 C[i][j] = A[i][j] + B[i][j];
106 complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
112 C = calloc(cols, sizeof(complex double*));
113 for(i=0; i<cols; i++)
114 C[i] = calloc(rows, sizeof(complex double));
116 for(i=0; i<rows; i++)
117 for(j=0; j<cols; j++)
118 C[i][j] = scalar*a[i][j];
123 complex double trace(complex double **a, int rows, int cols)
126 complex double tr = 0.0*I;
128 for(i=0; i<rows; i++)
134 complex double** transp(complex double **a, int rows, int cols)
139 C = calloc(cols, sizeof(complex double*));
140 for(i=0; i<cols; i++)
141 C[i] = calloc(rows, sizeof(complex double));
143 for(i=0; i<cols; i++)
144 for(j=0; j<rows; j++) {
151 complex double** conjtransp(complex double **a, int rows, int cols)
156 C = calloc(cols, sizeof(complex double*));
157 for(i=0; i<cols; i++)
158 C[i] = calloc(rows, sizeof(complex double));
160 for(i=0; i<cols; i++)
161 for(j=0; j<rows; j++) {
162 C[i][j] = conj(a[j][i]);
168 void printMatrix(complex double** a, int rows, int cols)
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]));
180 complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
184 C = calloc(ro1, sizeof(complex double*));
186 C[i] = calloc(co2, sizeof(complex double));
188 for(i=0; i<ro1; i++) {
189 for(j=0; j<co2; j++) {
192 C[i][j] += A[i][k] * B[k][j];
199 complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
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));
210 for(l=0; l<co2; l++) {
211 C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
219 void deallocate_mem(complex double ***arr, int rows)
222 for(i=0; i<rows; i++)
227 int readPaulicoeffs(int *alpha, int *beta, int *gamma, int *delta, int numqubits)
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]);