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);
10 // order of matrix elements is [row][column]!!!
17 int N; // number of qubits
20 int alpha[N], beta[N], gamma[N], delta[N];
22 double complex summand, sum;
24 while(readPaulicoeffs(alpha, beta, gamma, delta, N)) {
27 // printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
33 summand += 0.25*((cpow(I,delta[i+1]-gamma[i])*Gausssum1d(0,beta[i]+beta[i+1]+gamma[i]+delta[i+1])
34 +cpow(-1.0,beta[i]+delta[i])*cpow(I,delta[i]+delta[i+1])*Gausssum1d(0,beta[i]+beta[i+1]+delta[i]+delta[i+1]))
35 *Kroneck(gamma[i]+delta[i]+gamma[i+1]+delta[i+1])
36 +(csqrt(-I)*cpow(I,delta[i]+delta[i+1])*Gausssum1d(1,beta[i]+beta[i+1]+delta[i]+delta[i+1])
37 +csqrt(I)*cpow(-1.0,beta[i+1]+gamma[i+1])*cpow(I,delta[i]-gamma[i+1])*(-I)*Gausssum1d(1,beta[i]+beta[i+1]+delta[i]+gamma[i+1]))
38 *Kroneck(gamma[i]+delta[i]+gamma[i+1]+delta[i+1]+1));
41 //printf("%lf+%lfI\n", creal(sum), cimag(sum));
42 printf("%lf\n", cabs(creal(sum)));
50 complex double Kroneck(int arg)
52 arg = (arg+1)%2; // output 1 if argument is 0 mod 2 and 0 otherwise
53 return ((complex double)arg);
56 complex double Gausssum1d(int quadraticcoeff, int linearcoeff)
58 // we assume coeffs are either 0 or 1
62 if(quadraticcoeff == 0)
75 int readPaulicoeffs(int *alpha, int *beta, int *gamma, int *delta, int numqubits)
80 if(scanf("%d %d %d %d", &alpha[0], &beta[0], &gamma[0], &delta[0]) != EOF) {
81 for(i=1; i<numqubits; i++) {
82 scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]);