X-Git-Url: https://s3miclassical.com/gitweb/?p=weak_simulation_stab_extent.git;a=blobdiff_plain;f=multipauli.c;fp=multipauli.c;h=2161b72083f676d3055279186392eb0f2d5707c4;hp=0000000000000000000000000000000000000000;hb=7bc25656311c391af6a4d5814cc82080981c63c4;hpb=d914e7c183347b1340988e0869a6e053dd573991

diff --git a/multipauli.c b/multipauli.c
new file mode 100644
index 0000000..2161b72
--- /dev/null
+++ b/multipauli.c
@@ -0,0 +1,301 @@
+#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* omega, 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, j;
+  
+  int N;              // number of qubits
+  scanf("%d", &N);
+
+  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)
+  scanf("%d", &K);  
+
+  complex double **IN;  // N-qubit identity matrix
+  IN = I2;
+  for(i=1; i<N; i++) {
+    IN = outerMatrix(IN,I2,pow(2,i),pow(2,i),2,2);
+  }
+
+  complex double (*(psiT[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
+  complex double (*(psi0[])) = { (double complex[]) {1.0+0.0*I}, (double complex[]) {0.0*I}}; // T gate magic state
+
+  complex double **psiN;
+  if(K > 0) {
+    psiN = psiT;
+    for(i=1; i<K; i++) 
+      psiN = outerMatrix(psiN, psiT, pow(2,i), 1, 2, 1);
+    for(i=K; i<N; i++) 
+      psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+  } else {
+    psiN = psi0;
+    for(i=1; i<N; i++) 
+      psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+  }
+      
+  
+  int omega, alpha[N], beta[N], gamma[N], delta[N];
+
+  int Paulicounter = 0;
+
+  complex double **fullP;  // full product: \prod_i 1/2*(1+P_i)
+  complex double **P;      // P (P_i above) is made up of products of one-qubit Paulis, P1
+  complex double **P1[N];  // one-qubit Paulis
+
+  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(&omega, alpha, beta, gamma, delta, N)) { // go over the product of 1/2*(I+Paulis) that makes up the full projector
+
+    Paulicounter++;
+    if(Paulicounter > N) {
+      printf("Error: Number of Paulis is greater than N!\n");
+      return 1;
+    }
+      
+
+    printf("%d\n", omega);
+    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);
+    }
+    
+    P = P1[0];
+    for(i=1; i<N; i++)
+      P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
+    P = scalarmultMatrix(cpow(I,omega),P,pow(2,N),pow(2,N));
+    
+    if(Paulicounter == 1)
+      fullP = scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N));
+    else {
+      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));
+    }
+    deallocate_mem(&P, pow(2,N));
+
+  }
+
+  complex double **psiNfinal = multMatrix(fullP,psiN,pow(2,N),pow(2,N),pow(2,N),1);
+
+  printf("fullP:\n");
+  for(i=0; i<pow(2,N); i++) {
+    for(j=0; j<pow(2,N); j++) {
+      printf("%lf+%lfI ", creal(fullP[i][j]), cimag(fullP[i][j]));
+    }
+    printf("\n");
+  }
+  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));
+  }
+
+  if(creal(tr+0.00000001)>0)
+    printf("%.10lf %c %.10lf I\n", cabs(creal(tr)), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
+  else
+    printf("%.10lf %c %.10lf I\n", creal(tr), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
+  //printf("%lf %c %lf I\n", creal(tr), cimag(tr)>0?'+':'-' , cabs(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(&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 *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+  int i;
+
+  if(scanf("%d", omega) != EOF) {
+    for(i=0; i<numqubits; i++) {
+      if(scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]) == EOF) {
+	printf("Error: Too few input coeffs!\n");
+	exit(0);
+      }
+      if(alpha[i]+beta[i]+gamma[i]+delta[i] > 1) {
+	printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+	exit(0);
+      }
+    }
+    return 1;
+  } else
+    return 0;
+
+}
+
+
+