added code, README and Makefile
[strong_simulation_gauss_sum_rank.git] / gausssums_multipleof1.c
diff --git a/gausssums_multipleof1.c b/gausssums_multipleof1.c
new file mode 100644 (file)
index 0000000..3986ee0
--- /dev/null
@@ -0,0 +1,87 @@
+#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 N;              // number of qubits
+  scanf("%d", &N);
+
+  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++) {
+      summand = 0.0*I;
+      summand += 0.5*(Kroneck(gamma[i]+delta[i]) + pow(-1,beta[i])*Kroneck(gamma[i]+delta[i])
+                          +csqrt(I)*cpow(-I,delta[i])*Kroneck(alpha[i]+beta[i]) + csqrt(-I)*cpow(I,delta[i])*Kroneck(alpha[i]+beta[i]));
+      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)
+{
+  // we assume coeffs are either 0 or 1
+  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;
+
+}
+
+
+