X-Git-Url: https://s3miclassical.com/gitweb/?p=strong_simulation_stabilizer_rank.git;a=blobdiff_plain;f=strongsim2.c;fp=strongsim2.c;h=3b9b36b7c50e01578935beca832fc61a510546e8;hp=0000000000000000000000000000000000000000;hb=f36d60f8c85e0879a7d2b52e816638d2b4fdf86a;hpb=a85d4d567fc589503c5e263e1520295c8d433a45

diff --git a/strongsim2.c b/strongsim2.c
new file mode 100644
index 0000000..3b9b36b
--- /dev/null
+++ b/strongsim2.c
@@ -0,0 +1,210 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+#include "matrix.h"
+#include "exponentialsum.h"
+#include "shrink.h"
+#include "extend.h"
+#include "measurepauli.h"
+#include "innerproduct.h"
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+int main()
+{
+
+  int N;              // number of qubits
+  scanf("%d", &N);
+
+  if(N%2 != 0) {
+    printf("'N' needs to be an even number for a k=2 tensor factor decomposition!\n");
+    return 1;
+  }
+
+  int T;              // 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", &T);  
+
+  int omega;
+  int alpha[N], beta[N], gamma[N], delta[N];
+  int Paulicounter = 0;
+
+  int i, j, k, l, m;
+
+
+  int n1 = 2; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1}, (int[]) {0, 1} }; int (*(GBar1[])) = { (int[]) {1, 0}, (int[]) {1, 1} }; int h1[] = {0, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
+  int n2 = 2; int k2 = 1; int (*(G2[])) = { (int[]) {1, 1}, (int[]) {0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0}, (int[]) {1, 1} }; int h2[] = {0, 1}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
+
+  int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
+  double complex Gamma[(int)pow(2.0,N/2)]; // prefactor in front of resultant state
+  G = calloc(pow(2,N/2),sizeof(int*)); GBar = calloc(pow(2,N/2),sizeof(int*));
+  h = calloc(pow(2,N/2),sizeof(int*));
+  
+  J = calloc(pow(2,N/2),sizeof(int*)); D = calloc(pow(2,N/2),sizeof(int*)); Q = calloc(pow(2,N/2),sizeof(int));
+
+  K = calloc(pow(2,N/2), sizeof(int));
+
+  double complex origGamma[(int)pow(2,N/2)];
+  int *origK, *origQ, **origD, ***origJ;
+  int ***origG, ***origGBar, **origh;
+
+  origG = calloc(pow(2,N/2),sizeof(int*)); origGBar = calloc(pow(2,N/2),sizeof(int*));
+  origh = calloc(pow(2,N/2),sizeof(int*));
+  
+  origJ = calloc(pow(2,N/2),sizeof(int*)); origD = calloc(pow(2,N/2),sizeof(int*)); origQ = calloc(pow(2,N/2),sizeof(int));
+
+  origK = calloc(pow(2,N/2), sizeof(int));
+
+  int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+  
+
+  for(j=0; j<pow(2,N/2); j++) { // there will be 2^(N/2) combinations when using k=2 tensor factors
+
+    K[j] = k1*N/2; // assuming k1=k2
+    origK[j] = K[j];
+
+    combination = j;
+
+    Gamma[j] = 1.0;
+    
+    G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
+    h[j] = calloc(N, sizeof(int));
+
+    if(K[j] > 0) {
+      J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
+      for(k=0; k<K[j]; k++)
+	J[j][k] = calloc(K[j], sizeof(int));
+    }
+
+    origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
+    origh[j] = calloc(N, sizeof(int));
+    
+    if(K[j] > 0) {
+      origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
+      for(k=0; k<K[j]; k++)
+	origJ[j][k] = calloc(K[j], sizeof(int));
+    }
+
+    for(k=0; k<N/2; k++) {
+
+      Q[j] += (combination & 0x1)*Q2 + (~combination & 0x1)*Q1;
+      
+      Gamma[j] *= (1.0/sqrt(2.0))*((combination & 0x1)*cexp(0.25*PI*I) + (~combination & 0x1)*cexp(0.0*PI*I));
+      
+      G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
+      G[j][N/2+k] = calloc(N, sizeof(int)); GBar[j][N/2+k] = calloc(N, sizeof(int));
+
+      origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
+      origG[j][N/2+k] = calloc(N, sizeof(int)); origGBar[j][N/2+k] = calloc(N, sizeof(int));
+      
+      for(l=0; l<k1; l++) { // assuming k1=k2
+	D[j][k*k1+l] = (combination & 0x1)*D2[l] + (~combination & 0x1)*D1[l];
+	for(m=0; m<k1; m++) // assuming k1=k2
+	  J[j][k*k1+l][k*k1+m] = (combination & 0x1)*J2[l][m] + (~combination & 0x1)*J1[l][m];
+      }
+	    
+      for(l=0; l<n1; l++) { // assuming n1=n2
+	h[j][k*n1+l] = (combination & 0x1)*h2[l] + (~combination & 0x1)*h1[l];
+      	G[j][k][k*n1+l] = (combination & 0x1)*G2[0][l] + (~combination & 0x1)*G1[0][l];
+      	G[j][N/2+k][k*n1+l] = (combination & 0x1)*G2[1][l] + (~combination & 0x1)*G1[1][l]; // basis outside of support
+      	GBar[j][k][k*n1+l] = (combination & 0x1)*GBar2[0][l] + (~combination & 0x1)*GBar1[0][l];
+      	GBar[j][N/2+k][k*n1+l] = (combination & 0x1)*GBar2[1][l] + (~combination & 0x1)*GBar1[1][l]; // basis outside of support
+      }
+      memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
+      memcpy(origG[j][N/2+k], G[j][N/2+k], N*sizeof(int)); memcpy(origGBar[j][N/2+k], GBar[j][N/2+k], N*sizeof(int));
+
+      memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
+      
+      combination /= 2; // shift to the right by one (in base-2 arithmetic)
+    }
+
+    memcpy(origh[j], h[j], N*sizeof(int));
+    memcpy(origD[j], D[j], K[j]*sizeof(int));
+
+  }
+  memcpy(origGamma, Gamma, pow(2,N/2)*sizeof(double complex));
+
+  memcpy(origQ, Q, pow(2,N/2)*sizeof(int));
+
+  while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
+
+    Paulicounter++;
+    if(Paulicounter > N) {
+      printf("Error: Number of Paulis is greater than N!\n");
+      return 1;
+    }
+    
+    // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
+    // Y_i = -I*Z*X
+    for(i=0; i<N; i++) {
+      if(delta[i]){
+	omega += 3; // -I = I^3
+	beta[i] = delta[i];
+	gamma[i] = delta[i];
+      }
+    }
+
+    for(j=0; j<pow(2,N/2); j++) { // the kets
+
+      Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+
+    }
+
+  }
+
+  double complex amplitude = 0.0 + 0.0*I;
+  for(i=0; i<pow(2,N/2); i++) { // the bras
+    for(j=0; j<pow(2,N/2); j++) {
+      // check to see if second arguments are modified!!! They shouldn't be!
+      double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+      amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+    }
+  }
+
+  printf("amplitude:\n");
+  if(creal(amplitude+0.00000001)>0)
+    printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
+  else
+    printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
+  //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
+
+  return 0;
+
+}
+
+
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+    
+  int newomega, newalpha, newbeta, newgamma, newdelta;
+  int i;
+
+  if(scanf("%d", &newomega) != EOF) {
+    *omega = newomega;
+    for(i=0; i<numqubits; i++) {
+      if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
+	printf("Error: Too few input coeffs!\n");
+	exit(0);
+      }
+      if(newalpha+newbeta+newgamma+newdelta > 1) {
+	printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+	exit(0);
+      }
+      alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
+    }
+    return 1;
+  } else
+    return 0;
+    
+}
+
+
+
+
+
+
+