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);
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 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)
35 complex double **IN; // N-qubit identity matrix
38 IN = outerMatrix(IN,I2,pow(2,i),pow(2,i),2,2);
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
44 complex double **psiN;
48 psiN = outerMatrix(psiN, psiT, pow(2,i), 1, 2, 1);
50 psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
54 psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
58 int omega, alpha[N], beta[N], gamma[N], delta[N];
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
69 for(i=0; i<pow(2,N); i++) {
70 printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0]));
74 while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) { // go over the product of 1/2*(I+Paulis) that makes up the full projector
77 if(Paulicounter > N) {
78 printf("Error: Number of Paulis is greater than N!\n");
83 printf("%d\n", omega);
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);
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));
95 fullP = scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N));
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));
99 deallocate_mem(&P, pow(2,N));
103 complex double **psiNfinal = multMatrix(fullP,psiN,pow(2,N),pow(2,N),pow(2,N),1);
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]));
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]));
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));
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)));
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 */
132 deallocate_mem(&psiNfinal, pow(2,N));
135 // deallocate_mem(&psiN, pow(2,N));
141 complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
147 C = calloc(cols, sizeof(complex double*));
148 for(i=0; i<cols; i++)
149 C[i] = calloc(rows, sizeof(complex double));
151 for(i=0; i<rows; i++)
152 for(j=0; j<cols; j++)
153 C[i][j] = A[i][j] + B[i][j];
158 complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
164 C = calloc(cols, sizeof(complex double*));
165 for(i=0; i<cols; i++)
166 C[i] = calloc(rows, sizeof(complex double));
168 for(i=0; i<rows; i++)
169 for(j=0; j<cols; j++)
170 C[i][j] = scalar*a[i][j];
175 complex double trace(complex double **a, int rows, int cols)
178 complex double tr = 0.0*I;
180 for(i=0; i<rows; i++)
186 complex double** transp(complex double **a, int rows, int cols)
191 C = calloc(cols, sizeof(complex double*));
192 for(i=0; i<cols; i++)
193 C[i] = calloc(rows, sizeof(complex double));
195 for(i=0; i<cols; i++)
196 for(j=0; j<rows; j++) {
203 complex double** conjtransp(complex double **a, int rows, int cols)
208 C = calloc(cols, sizeof(complex double*));
209 for(i=0; i<cols; i++)
210 C[i] = calloc(rows, sizeof(complex double));
212 for(i=0; i<cols; i++)
213 for(j=0; j<rows; j++) {
214 C[i][j] = conj(a[j][i]);
220 void printMatrix(complex double** a, int rows, int cols)
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]));
232 complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
236 C = calloc(ro1, sizeof(complex double*));
238 C[i] = calloc(co2, sizeof(complex double));
240 for(i=0; i<ro1; i++) {
241 for(j=0; j<co2; j++) {
244 C[i][j] += A[i][k] * B[k][j];
251 complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
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));
263 C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
270 void deallocate_mem(complex double ***arr, int rows)
273 for(i=0; i<rows; i++)
278 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
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");
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);