X-Git-Url: https://s3miclassical.com/gitweb/?p=strong_simulation_gauss_sum_rank.git;a=blobdiff_plain;f=hilbertspace_vector.c;fp=hilbertspace_vector.c;h=15d143d2ed00a6cb33c6e535c3745a29e7474e4d;hp=0000000000000000000000000000000000000000;hb=41acfb8cfb426c10f7a16e8cb04337012b25e3fb;hpb=13d8b4a4c252c341d36d3d90f66ed28600777224

diff --git a/hilbertspace_vector.c b/hilbertspace_vector.c
new file mode 100644
index 0000000..15d143d
--- /dev/null
+++ b/hilbertspace_vector.c
@@ -0,0 +1,243 @@
+#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;
+
+}
+
+
+