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

diff --git a/gausssums_multipleof6.c b/gausssums_multipleof6.c
new file mode 100644
index 0000000..b86aa4b
--- /dev/null
+++ b/gausssums_multipleof6.c
@@ -0,0 +1,176 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <complex.h>
+#include <math.h>
+
+int readPaulicoeffs(int* alpha, int* beta, int* gamma, int* delta, int numqubits);
+complex double Gausssum1d(int quadraticcoeff, int linearcoeff);
+complex double Kroneck(int arg);
+  
+// order of matrix elements is [row][column]!!!
+
+int main()
+{
+
+  int i;
+  int x, y, xp, yp;
+  
+  int N;              // number of qubits
+  scanf("%d", &N);
+  if(N%6!=0) {
+    printf("Error: N must be a multiple of 6!\n");
+    return 1;
+  }
+
+  int alpha[N], beta[N], gamma[N], delta[N];
+
+  double complex summand, sum;
+
+  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]); */
+
+    sum = 1.0+0.0*I;
+    for(i=0; i<N; i+=6) {
+      summand = 0.0*I;
+      for(y=0; y<=1; y++) {
+      	//for(x=0;x<=(y*(gamma[i]+gamma[i+1])+(y+1)*(gamma[i]+delta[i+1]))%2; x++) {
+      	for(x=y*(beta[i]+alpha[i+1]+beta[i+1])%2;x<=(y*(beta[i]+alpha[i+1]+beta[i+1])+y*(gamma[i]+gamma[i+1])+(y+1)*(gamma[i]+delta[i+1]))%2; x++) {
+	  for(yp=((gamma[i+3]+delta[i+4])*x)%2; yp<=(1+(gamma[i+3]+gamma[i+4])*x)%2; yp++) {
+	    //strange that I don't need the following modification to x
+	  //****for(yp=((gamma[i+3]+delta[i+4])*x*(1+alpha[i]+beta[i]))%2; yp<=(1+(gamma[i+3]+gamma[i+4])*x*(1+alpha[i]+beta[i]))%2; yp++) {
+	    for(xp=yp*(beta[i+3]+alpha[i+4]+beta[i+4])%2;xp<=(yp*(beta[i+3]+alpha[i+4]+beta[i+4])+yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	  //for(xp=0;xp<=(yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	      summand += 2.0*(0.125*pow(2.0,(1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*(1-yp*pow(gamma[i+3]-gamma[i+4],2)-pow(yp-1,2)*pow(gamma[i+3]-delta[i+4],2)))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*Gausssum1d(0.0,beta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(gamma[i+2]+delta[i+2]) 
+			      + 0.125*pow(2.0,(1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*(1-yp*pow(gamma[i+3]-gamma[i+4],2)-pow(yp-1,2)*pow(gamma[i+3]-delta[i+4],2)))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*csqrt(-I)*Gausssum1d(gamma[i+2]+delta[i+2],0)*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*Gausssum1d(0.0,beta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(gamma[i+5]+delta[i+5]) 
+		  + 0.125*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*csqrt(-I)*Gausssum1d(gamma[i+5]+delta[i+5],0)*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+	  for(yp=0; yp<=(1+alpha[i+5])%2; yp++) {
+	    for(xp=0;xp<=(1+beta[i+5])%2; xp++) {
+	      summand += (0.125*pow(2.0,1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*Gausssum1d(0.0,beta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(gamma[i+2]+delta[i+2]) 
+			  + 0.125*pow(2.0,1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*csqrt(-I)*Gausssum1d(gamma[i+2]+delta[i+2],0)*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+5]*yp+pow(xp+1,2)*(delta[i+3]+gamma[i+4]+2.0*beta[i+4])+xp*(delta[i+3]+delta[i+4])))*Gausssum1d(1,beta[i+3]+beta[i+4]+delta[i+3]+(xp+1)*gamma[i+4]+xp*delta[i+4])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(gamma[i+5]+delta[i+5]));
+	    }
+	  }
+	  //for(yp=x; yp<=1; yp++) {
+	  for(yp=x*(1+alpha[i]+beta[i])%2; yp<=1; yp++) {
+	    //for(xp=(yp*(delta[i+3]+delta[i+4])+x)%2;xp<=(1+yp*(gamma[i+3]+gamma[i+4])+x)%2; xp++) {
+	    for(xp=(yp*(delta[i+3]+delta[i+4])+x*(1+alpha[i]+beta[i]))%2;xp<=(1+yp*(gamma[i+3]+gamma[i+4])+x*(1+alpha[i]+beta[i]))%2; xp++) {
+	      summand += 2.0*(0.125*pow(2.0,(1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*(pow(yp-1,2)*(pow(xp,2)*(gamma[i+3]+delta[i+4])+pow(xp-1,2)*(gamma[i+4]+delta[i+3]))))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*Gausssum1d(0.0,beta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(gamma[i+2]+delta[i+2]) 
+			  + 0.125*pow(2.0,(1-y*pow(gamma[i]-gamma[i+1],2)-pow(y-1,2)*pow(gamma[i]-delta[i+1],2))*(pow(yp-1,2)*(pow(xp,2)*(gamma[i+3]+delta[i+4])+pow(xp-1,2)*(gamma[i+4]+delta[i+3]))))*cexp(M_PI*I/2.0*((delta[i+1]-gamma[i])*pow(y+1,2)+2.0*(beta[i]+beta[i+1]+gamma[i]+delta[i+1])*x*pow(y+1,2)+2.0*(delta[i]+delta[i+1]+beta[i]+beta[i+1])*x*y+(2.0*beta[i]+3.0*delta[i]+delta[i+1])*y))*csqrt(-I)*Gausssum1d(gamma[i+2]+delta[i+2],0)*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1])*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*(-I)*cexp(M_PI*I/2.0*(yp*(delta[i+3]+delta[i+4])+pow(yp+1,2)*(delta[i+3]+gamma[i+4]+2*beta[i+4])+pow(xp,2)+2*(beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp)*xp))*Gausssum1d(0,xp+beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp+gamma[i+5]+delta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+      	}
+      }
+
+      for(y=0; y<=(1+alpha[i+2])%2; y++) {
+	for(x=0; x<=(1+beta[i+2])%2; x++) {
+	  for(yp=0; yp<=1; yp++) {
+	    for(xp=yp*(beta[i+3]+alpha[i+4]+beta[i+4])%2;xp<=(yp*(beta[i+3]+alpha[i+4]+beta[i+4])+yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	      //for(xp=0;xp<=(yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	      summand += (0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+2]*y+pow(x+1,2)*(delta[i]+gamma[i+1]+2.0*beta[i+1])+x*(delta[i]+delta[i+1])))*Gausssum1d(1,beta[i]+beta[i+1]+delta[i]+(x+1)*gamma[i+1]+x*delta[i+1])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(gamma[i+2]+delta[i+2]))
+		*(0.125*pow(2.0,1-yp*pow(gamma[i+3]-gamma[i+4],2)-pow(yp-1,2)*pow(gamma[i+3]-delta[i+4],2))*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*Gausssum1d(0.0,beta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(gamma[i+5]+delta[i+5]) 
+		  + 0.125*pow(2.0,1-yp*pow(gamma[i+3]-gamma[i+4],2)-pow(yp-1,2)*pow(gamma[i+3]-delta[i+4],2))*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*csqrt(-I)*Gausssum1d(gamma[i+5]+delta[i+5],0)*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+	  for(yp=0; yp<=(1+alpha[i+5])%2; yp++) {
+	    for(xp=0;xp<=(1+beta[i+5])%2; xp++) {
+	      summand += (0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+2]*y+pow(x+1,2)*(delta[i]+gamma[i+1]+2.0*beta[i+1])+x*(delta[i]+delta[i+1])))*Gausssum1d(1,beta[i]+beta[i+1]+delta[i]+(x+1)*gamma[i+1]+x*delta[i+1])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(gamma[i+2]+delta[i+2]))
+		*(0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+5]*yp+pow(xp+1,2)*(delta[i+3]+gamma[i+4]+2.0*beta[i+4])+xp*(delta[i+3]+delta[i+4])))*Gausssum1d(1,beta[i+3]+beta[i+4]+delta[i+3]+(xp+1)*gamma[i+4]+xp*delta[i+4])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(gamma[i+5]+delta[i+5]));
+	    }
+	  }
+	  for(yp=0; yp<=1; yp++) {
+	    for(xp=(yp*(delta[i+3]+delta[i+4]))%2;xp<=(1+yp*(gamma[i+3]+gamma[i+4]))%2; xp++) {
+	      summand += (0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+2]*y+pow(x+1,2)*(delta[i]+gamma[i+1]+2.0*beta[i+1])+x*(delta[i]+delta[i+1])))*Gausssum1d(1,beta[i]+beta[i+1]+delta[i]+(x+1)*gamma[i+1]+x*delta[i+1])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(gamma[i+2]+delta[i+2]))
+		*(0.125*(-I)*pow(2.0,(pow(yp-1,2)*(pow(xp,2)*(gamma[i+3]+delta[i+4])+pow(xp-1,2)*(gamma[i+4]+delta[i+3]))))*cexp(M_PI*I/2.0*(yp*(delta[i+3]+delta[i+4])+pow(yp+1,2)*(delta[i+3]+gamma[i+4]+2*beta[i+4])+pow(xp,2)+2*(beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp)*xp))*Gausssum1d(0,xp+beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp+gamma[i+5]+delta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+      	}
+      }
+
+      for(y=0; y<=1; y++) {
+	for(x=(y*(delta[i]+delta[i+1]))%2; x<=(1+y*(gamma[i]+gamma[i+1]))%2; x++) {
+	  for(yp=((gamma[i+3]+delta[i+4])*y)%2; yp<=(1+(gamma[i+3]+gamma[i+4])*y)%2; yp++) {
+	    for(xp=yp*(beta[i+3]+alpha[i+4]+beta[i+4])%2;xp<=(yp*(beta[i+3]+alpha[i+4]+beta[i+4])+yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	      //for(xp=0;xp<=(yp*(gamma[i+3]+gamma[i+4])+(yp+1)*(gamma[i+3]+delta[i+4]))%2; xp++) {
+	      summand += 2.0*(0.125*pow(2.0,(pow(y-1,2)*(pow(x,2)*(gamma[i]+delta[i+1])+pow(x-1,2)*(gamma[i+1]+delta[i])))*(1-yp*pow(gamma[i+3]-gamma[i+4],2)-pow(yp-1,2)*pow(gamma[i+3]-delta[i+4],2)))*(-I)*cexp(M_PI*I/2.0*(y*(delta[i]+delta[i+1])+pow(y+1,2)*(delta[i]+gamma[i+1]+2*beta[i+1])+pow(x,2)+2*(beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y)*x))*Gausssum1d(0,x+beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y+gamma[i+2]+delta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*Gausssum1d(0.0,beta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(gamma[i+5]+delta[i+5]) 
+		  + 0.125*cexp(M_PI*I/2.0*((delta[i+4]-gamma[i+3])*pow(yp+1,2)+2.0*(beta[i+3]+beta[i+4]+gamma[i+3]+delta[i+4])*xp*pow(yp+1,2)+2.0*(delta[i+3]+delta[i+4]+beta[i+3]+beta[i+4])*xp*yp+(2.0*beta[i+3]+3.0*delta[i+3]+delta[i+4])*yp))*csqrt(-I)*Gausssum1d(gamma[i+5]+delta[i+5],0)*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4])*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+	  for(yp=0; yp<=(1+alpha[i+5])%2; yp++) {
+	    for(xp=0;xp<=(1+beta[i+5])%2; xp++) {
+	      summand += (0.125*pow(2.0,pow(y-1,2)*(pow(x,2)*(gamma[i]+delta[i+1])+pow(x-1,2)*(gamma[i+1]+delta[i])))*(-I)*cexp(M_PI*I/2.0*(y*(delta[i]+delta[i+1])+pow(y+1,2)*(delta[i]+gamma[i+1]+2*beta[i+1])+pow(x,2)+2*(beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y)*x))*Gausssum1d(0,x+beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y+gamma[i+2]+delta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*2.0*csqrt(-I)*cexp(M_PI*I/2.0*(2.0*beta[i+5]*yp+pow(xp+1,2)*(delta[i+3]+gamma[i+4]+2.0*beta[i+4])+xp*(delta[i+3]+delta[i+4])))*Gausssum1d(1,beta[i+3]+beta[i+4]+delta[i+3]+(xp+1)*gamma[i+4]+xp*delta[i+4])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(gamma[i+5]+delta[i+5]));
+	    }
+	  }
+	  for(yp=y; yp<=1; yp++) {
+	    for(xp=(yp*(delta[i+3]+delta[i+4])+y)%2;xp<=(1+yp*(gamma[i+3]+gamma[i+4])+y)%2; xp++) {
+	      summand += 2.0*(0.125*pow(2.0,(pow(y-1,2)*(pow(x,2)*(gamma[i]+delta[i+1])+pow(x-1,2)*(gamma[i+1]+delta[i])))*(pow(yp-1,2)*(pow(xp,2)*(gamma[i+3]+delta[i+4])+pow(xp-1,2)*(gamma[i+4]+delta[i+3]))))*(-I)*cexp(M_PI*I/2.0*(y*(delta[i]+delta[i+1])+pow(y+1,2)*(delta[i]+gamma[i+1]+2*beta[i+1])+pow(x,2)+2*(beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y)*x))*Gausssum1d(0,x+beta[i]+beta[i+1]+delta[i]+gamma[i+1]*(y+1)+delta[i+1]*y+gamma[i+2]+delta[i+2])*Kroneck(gamma[i]+delta[i]-gamma[i+1]-delta[i+1]+1)*Kroneck(alpha[i+2]+beta[i+2]))
+		*(0.125*(-I)*cexp(M_PI*I/2.0*(yp*(delta[i+3]+delta[i+4])+pow(yp+1,2)*(delta[i+3]+gamma[i+4]+2*beta[i+4])+pow(xp,2)+2*(beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp)*xp))*Gausssum1d(0,xp+beta[i+3]+beta[i+4]+delta[i+3]+gamma[i+4]*(yp+1)+delta[i+4]*yp+gamma[i+5]+delta[i+5])*Kroneck(gamma[i+3]+delta[i+3]-gamma[i+4]-delta[i+4]+1)*Kroneck(alpha[i+5]+beta[i+5]));
+	    }
+	  }
+      	}
+      }
+
+      sum *= summand;
+    }
+    //printf("%lf+%lfI\n", creal(sum), cimag(sum));
+    printf("%lf\n", cabs(creal(sum))); 
+
+
+  }
+
+  return 0;
+}
+
+complex double Kroneck(int arg)
+{
+  arg = (arg+1)%2; // output 1 if argument is 0 mod 2 and 0 otherwise
+  return ((complex double)arg);
+}
+
+complex double Gausssum1d(int quadraticcoeff, int linearcoeff)
+{
+  /*****************************************
+  /* NOTE! we assume coeffs are either 0 or 1
+  /* (So you cannot pass off a linear coeff as a quadratic coeff by multiplying it by a factor of 2 (and considering it as mod 4)!)
+  *****************************************/
+  
+  quadraticcoeff %= 2;
+  linearcoeff %= 2;
+    
+  if(quadraticcoeff == 0)
+    if(linearcoeff == 0)
+      return 2.0+0.0*I;
+    else
+      return 0.0*I;
+  else
+    if(linearcoeff == 0)
+      return 1.0+1.0*I;
+    else
+      return 1.0-1.0*I;
+
+}
+
+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;
+
+}
+
+
+