gausssums_multipleof3.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <complex.h>
   4 #include <math.h>
   5
   6 int readPaulicoeffs(int* alpha, int* beta, int* gamma, int* delta, int numqubits);
   7 complex double Gausssum1d(int quadraticcoeff, int linearcoeff);
   8 complex double Kroneck(int arg);
   9
  10 // order of matrix elements is [row][column]!!!
  11
  12 int main()
  13 {
  14
  15   int i, x, y;
  16
  17   int N;              // number of qubits
  18   scanf("%d", &N);
  19   if(N%3!=0) {
  20     printf("Error: N must be a multiple of 3!\n");
  21     return 1;
  22   }
  23
  24   int alpha[N], beta[N], gamma[N], delta[N];
  25
  26   double complex summand, sum;
  27
  28   while(readPaulicoeffs(alpha, beta, gamma, delta, N)) {
  29
  30     /* for(i=0; i<N; i++) */
  31     /*   printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]); */
  32
  33     sum = 1.0+0.0*I;
  34     for(i=0; i<N; i+=3) {
  35       summand = 0.0*I;
  36       for(y=0; y<=1; y++) {
  37         //for(x=0;x<=(y*(gamma[i]+gamma[i+1])+(y+1)*(gamma[i]+delta[i+1]))%2; x++) {
  38         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++) {
  39           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]);
  40           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))*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]);
  41           //summand += 0.125*pow(2.0,gamma[i]+delta[i]+x*y*(beta[i]+alpha[i])*alpha[i+1]+pow(x-1,2)*y*(alpha[i]+beta[i])*beta[i+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]);
  42           //}
  43         //for(x=((y+1)*(alpha[i]+beta[i]))%2;x<=((beta[i]+alpha[i])+y*(gamma[i]+gamma[i+1])+(y+1)*(gamma[i]+delta[i+1]))%2; x++) {
  44           //summand += 0.125*pow(2.0,gamma[i]+delta[i]+x*y*(beta[i]+alpha[i])*alpha[i+1]+pow(x-1,2)*y*(alpha[i]+beta[i])*beta[i+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]);
  45         }
  46       }
  47       for(y=0; y<=(1+alpha[i+2])%2; y++) {
  48         for(x=0;x<=(1+beta[i+2])%2; x++) {
  49           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]);
  50         }
  51       }
  52       for(y=0; y<=1; y++) {
  53         for(x=(y*(delta[i]+delta[i+1]))%2;x<=(1+y*(gamma[i]+gamma[i+1]))%2; x++) {
  54           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]);
  55         }
  56       }
  57       sum *= summand;
  58     }
  59     //printf("%lf+%lfI\n", creal(sum), cimag(sum));
  60     printf("%lf\n", cabs(creal(sum)));
  61
  62
  63   }
  64
  65   return 0;
  66 }
  67
  68 complex double Kroneck(int arg)
  69 {
  70   arg = (arg+1)%2; // output 1 if argument is 0 mod 2 and 0 otherwise
  71   return ((complex double)arg);
  72 }
  73
  74 complex double Gausssum1d(int quadraticcoeff, int linearcoeff)
  75 {
  76   /*****************************************
  77   /* NOTE! we assume coeffs are either 0 or 1
  78   /* (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)!)
  79   *****************************************/
  80
  81   quadraticcoeff %= 2;
  82   linearcoeff %= 2;
  83
  84   if(quadraticcoeff == 0)
  85     if(linearcoeff == 0)
  86       return 2.0+0.0*I;
  87     else
  88       return 0.0*I;
  89   else
  90     if(linearcoeff == 0)
  91       return 1.0+1.0*I;
  92     else
  93       return 1.0-1.0*I;
  94
  95 }
  96
  97 int readPaulicoeffs(int *alpha, int *beta, int *gamma, int *delta, int numqubits)
  98 {
  99
 100   int i;
 101
 102   if(scanf("%d %d %d %d", &alpha[0], &beta[0], &gamma[0], &delta[0]) != EOF) {
 103     for(i=1; i<numqubits; i++) {
 104       scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]);
 105     }
 106     return 1;
 107   } else
 108     return 0;
 109
 110 }
 111
 112
 113