measurepauli.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <complex.h>
   5
   6 #include "matrix.h"
   7 #include "shrink.h"
   8 #include "extend.h"
   9 #include "measurepauli.h"
  10
  11
  12 /****************************************************************************
  13  * Takes a stabilizer state (n, k, h, G, GBar, Q, D, J) and a Pauli operator
  14  * (m, X, Z) and returns the norm (abs. value) of the projected state.
  15  * If it is non-zero,then it transforms the stabilizer state to the
  16  * projected state.
  17  ****************************************************************************/
  18 // This assumes the order Z-first X-second order when operator is read from left-to-right!
  19 double measurepauli(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int m, int *X, int *Z) {
  20   // Note that X is 'xi' and Z is 'zeta' in Bravyi et al. 'xi' and 'zeta' are instead used for the vectors with elements xi_a and zeta_a.
  21   // Note that D and J must be dynamically allocated arrays as they can change size in this function. This is also the reason that we need their base pointer and are greater in pointer depth by one.
  22
  23   int a, b;
  24
  25   int omega;
  26
  27   int *xi, *zeta, *Xprime;
  28   xi = calloc(*k, sizeof(int));
  29   zeta = calloc(*k, sizeof(int));
  30   Xprime = calloc(n, sizeof(int));
  31
  32   for(a=0; a<*k; a++) {
  33     xi[a] = dotProductMod(GBar[a], X, n, 2);
  34     zeta[a] = dotProductMod(G[a], Z, n, 2);
  35   }
  36
  37   for(a=0; a<n; a++) {
  38     Xprime[a] = 0;
  39     for(b=0; b<*k; b++)
  40       Xprime[a] = Xprime[a] + xi[b]*G[b][a];
  41     Xprime[a] = Xprime[a]%2;
  42   }
  43
  44   omega = 2*m + 4*dotProductMod(Z, h, n, 2);
  45   for(a=0; a<*k; a++) {
  46     omega += (*D)[a]*xi[a];
  47     for(b=a+1; b<*k; b++)
  48       omega += (*J)[a][b]*xi[a]*xi[b];
  49   }
  50   omega = omega%8;
  51
  52   unsigned int xiEqual = 1; // set to 'true' to begin with
  53   for(a=0; a<n; a++)
  54     if(Xprime[a] != X[a])
  55       xiEqual = 0; // xiEqual is 'false'
  56
  57   if(xiEqual && (omega%4 == 0)) { // if xiEqual is true and omega in {0, 4}
  58     int* eta;
  59     eta = calloc(*k, sizeof(int));
  60     for(a=0; a<*k; a++) {
  61       eta[a] = 4*zeta[a];
  62       for(b=0; b<*k; b++)
  63         eta[a] = eta[a] + (*J)[a][b]*xi[b];
  64       eta[a] = (eta[a]/4)%2;
  65     }
  66     int* gamma;
  67     gamma = calloc(n, sizeof(int));
  68     for(a=0; a<n; a++) {
  69       gamma[a] = 0;
  70       for(b=0; b<*k; b++)
  71         gamma[a] = gamma[a] + eta[b]*GBar[b][a];
  72       gamma[a] = gamma[a]%2;
  73     }
  74     int alpha = (omega/4 + dotProductMod(gamma, h, n, 2))%2;
  75     char shrinkResult = shrink(n, k, h, G, GBar, Q, D, J, gamma, alpha);
  76     if(shrinkResult == 'E') {// empty
  77       free(gamma);
  78       free(eta);
  79       free(xi); free(zeta); free(Xprime);
  80       return(0.0); // Gamma = 0.0
  81     } else if(shrinkResult == 'A') {//same
  82       free(gamma);
  83       free(eta);
  84       free(xi); free(zeta); free(Xprime);
  85       return(1.0); // Gamma = 1.0
  86     } else {// success
  87       free(gamma);
  88       free(eta);
  89       free(xi); free(zeta); free(Xprime);
  90       return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
  91     }
  92   } else if(xiEqual && (omega+2)%4 == 0) { // if xiEqual is true and omega in {2, 6}
  93     int* eta;
  94     eta = calloc(*k, sizeof(int));
  95     for(a=0; a<*k; a++) {
  96       eta[a] = 4*zeta[a];
  97       for(b=0; b<*k; b++)
  98         eta[a] = eta[a] + (*J)[a][b]*xi[b];
  99       eta[a] = (eta[a]/4)%2;
 100     }
 101     int sigma = 2 - omega/2;
 102     *Q = (*Q + sigma)%8 < 0 ? (*Q + sigma)%8 + 8 : (*Q + sigma)%8;
 103     for(a=0; a<*k; a++) {
 104       (*D)[a] = ((*D)[a] - 2*sigma*eta[a])%8 < 0 ? ((*D)[a] - 2*sigma*eta[a])%8 + 8 : ((*D)[a] - 2*sigma*eta[a])%8;
 105       for(b=0; b<*k; b++) {
 106         if(a != b)
 107           (*J)[a][b] = ((*J)[a][b] + 4*eta[a]*eta[b])%8;
 108         else
 109           (*J)[a][a] = (2*((*D)[a]))%8;
 110       }
 111     }
 112     free(eta); free(xi); free(zeta); free(Xprime);
 113     return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
 114   } else { // if xiEqual is false
 115     extend(n, k, h, G, GBar, X);
 116     int *newD = calloc(*k, sizeof(int));
 117     for(a=0; a<(*k-1); a++) {
 118       newD[a] = (*D)[a];
 119     }
 120     int *tmpVec = calloc(n, sizeof(int));
 121     addVectorMod(h,X,tmpVec,n,2);
 122     newD[*k-1] = (2*m + 4*dotProductMod(Z,tmpVec,n,2))%8;
 123     free(tmpVec);
 124     if(*k-1 > 0)
 125       free(*D);
 126     *D = calloc(*k, sizeof(int));
 127     for(a=0; a<*k; a++) {
 128       (*D)[a] = newD[a];
 129     }
 130     free(newD);
 131
 132     int **tmpMatrix = calloc(*k, sizeof(int*));
 133     for(a=0; a<(*k-1); a++) {
 134       tmpMatrix[a] = calloc(*k, sizeof(int));
 135       for(b=0; b<(*k-1); b++) {
 136         tmpMatrix[a][b] = (*J)[a][b];
 137       }
 138     }
 139     tmpMatrix[*k-1] = calloc(*k, sizeof(int));
 140     for(a=0; a<(*k-1); a++) {
 141       tmpMatrix[*k-1][a] = 4*zeta[a];
 142       tmpMatrix[a][*k-1] = 4*zeta[a];
 143     }
 144     tmpMatrix[*k-1][*k-1] = (4*m)%8;
 145     if(*k-1 > 0)
 146       deallocate_mem(J, *k-1);
 147     //for(a=0; a<(*k-1); a++)
 148     //  free((*J)[a]);
 149     //free(*J);
 150     *J = calloc(*k, sizeof(int*));
 151     for(a=0; a<*k; a++) {
 152       (*J)[a] = calloc(*k, sizeof(int));
 153       for(b=0; b<*k; b++) {
 154         (*J)[a][b] = tmpMatrix[a][b];
 155       }
 156     }
 157     deallocate_mem(&tmpMatrix, *k);
 158
 159     free(xi); free(zeta); free(Xprime);
 160     return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
 161   }
 162
 163   return(0.0);
 164
 165 }