--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <complex.h>
+#include <math.h>
+
+void deallocate_mem(complex double ***arr, int rows);
+void printMatrix(complex double** a, int rows, int cols);
+complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
+complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
+complex double** transp(complex double **a, int rows, int cols);
+complex double** conjtransp(complex double **a, int rows, int cols);
+complex double trace(complex double **a, int rows, int cols);
+complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols);
+complex double** addMatrix(complex double **A, complex double **B, int rows, int cols);
+int readPaulicoeffs(int* alpha, int* beta, int* gamma, int* delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+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} };
+static complex double (*(X[])) = { (double complex[]) {0.0*I, 1.0+0.0*I}, (double complex[]) {1.0+0.0*I, 0.0*I} };
+static complex double (*(Y[])) = { (double complex[]) {0.0*I, 0.0-1.0*I}, (double complex[]) {0.0+1.0*I, 0.0*I} };
+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} };
+
+int main()
+{
+
+ int i;
+
+ int N; // number of qubits
+ scanf("%d", &N);
+
+ complex double (*(psi[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
+
+ complex double **psiN = psi;
+ for(i=1; i<N; i++)
+ psiN = outerMatrix(psiN, psi, pow(2,i), 1, 2, 1);
+
+ int alpha[N], beta[N], gamma[N], delta[N];
+
+ complex double **P1[N];
+
+ complex double tr;
+
+ /* printf("psiN:\n"); */
+ /* for(i=0; i<pow(2,N); i++) { */
+ /* printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0])); */
+ /* } */
+
+ while(readPaulicoeffs(alpha, beta, gamma, delta, N)) {
+
+ for(i=0; i<N; i++) {
+ //printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
+ 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);
+ }
+
+ complex double **P = P1[0];
+ for(i=1; i<N; i++)
+ P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
+
+ complex double **psiNfinal = multMatrix(P,psiN,pow(2,N),pow(2,N),pow(2,N),1);
+
+ /* printf("psiNfinal:\n"); */
+ /* for(i=0; i<pow(2,N); i++) { */
+ /* printf("%d: %lf+%lfI\n", i, creal(psiNfinal[i][0]), cimag(psiNfinal[i][0])); */
+ /* } */
+
+ tr = 0.0 + 0.0*I;
+ //printf("tr:\n");
+ for(i=0; i<pow(2,N); i++) {
+ tr += conj(psiN[i][0])*psiNfinal[i][0];
+ //printf("%d: %lf+%lfI\n", i, creal(tr), cimag(tr));
+ }
+
+ // printf("%lf+%lfI\n", cabs(creal(tr)), cimag(tr));
+ printf("%lf\n", cabs(creal(tr)));
+ /* 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 */
+
+ deallocate_mem(&P, pow(2,N));
+ deallocate_mem(&psiNfinal, pow(2,N));
+
+ }
+
+ // deallocate_mem(&psiN, pow(2,N));
+
+ return 0;
+}
+
+
+complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
+{
+ int i, j;
+
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = A[i][j] + B[i][j];
+
+ return C;
+}
+
+complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
+{
+ int i, j;
+
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = scalar*a[i][j];
+
+ return C;
+}
+
+complex double trace(complex double **a, int rows, int cols)
+{
+ int i;
+ complex double tr = 0.0*I;
+
+ for(i=0; i<rows; i++)
+ tr += a[i][i];
+
+ return tr;
+}
+
+complex double** transp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = a[j][i];
+ }
+
+ return C;
+}
+
+complex double** conjtransp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = conj(a[j][i]);
+ }
+
+ return C;
+}
+
+void printMatrix(complex double** a, int rows, int cols)
+{
+ int i, j;
+ printf("Matrix[%d][%d]\n", rows, cols);
+ for(i=0; i<rows; i++) {
+ for(j=0; j<cols; j++) {
+ printf("%lf+%lfI ", creal(a[i][j]), cimag(a[i][j]));
+ }
+ printf("\n");
+ }
+}
+
+complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k;
+ complex double **C;
+ C = calloc(ro1, sizeof(complex double*));
+ for(i=0; i<ro1; i++)
+ C[i] = calloc(co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++) {
+ for(j=0; j<co2; j++) {
+ C[i][j] = 0;
+ for(k=0; k<co1; k++)
+ C[i][j] += A[i][k] * B[k][j];
+ }
+ }
+
+ return C;
+}
+
+complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k, l;
+ complex double **C;
+ C = calloc(ro1*ro2, sizeof(complex double*));
+ for(i=0; i<ro1*ro2; i++)
+ C[i] = calloc(co1*co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ for(k=0; k<co1; k++)
+ for(l=0; l<co2; l++) {
+ C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
+ }
+
+ return C;
+}
+
+
+
+void deallocate_mem(complex double ***arr, int rows)
+{
+ int i;
+ for(i=0; i<rows; i++)
+ free((*arr)[i]);
+ free(*arr);
+}
+
+int readPaulicoeffs(int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int i;
+
+ if(scanf("%d %d %d %d", &alpha[0], &beta[0], &gamma[0], &delta[0]) != EOF) {
+ for(i=1; i<numqubits; i++) {
+ scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]);
+ }
+ return 1;
+ } else
+ return 0;
+
+}
+
+
+
projects / strong_simulation_gauss_sum_rank.git / blobdiff