9 #include "measurepauli.h"
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
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.
27 int *xi, *zeta, *Xprime;
28 xi = calloc(*k, sizeof(int));
29 zeta = calloc(*k, sizeof(int));
30 Xprime = calloc(n, sizeof(int));
33 xi[a] = dotProductMod(GBar[a], X, n, 2);
34 zeta[a] = dotProductMod(G[a], Z, n, 2);
40 Xprime[a] = Xprime[a] + xi[b]*G[b][a];
41 Xprime[a] = Xprime[a]%2;
44 omega = 2*m + 4*dotProductMod(Z, h, n, 2);
46 omega += (*D)[a]*xi[a];
48 omega += (*J)[a][b]*xi[a]*xi[b];
52 unsigned int xiEqual = 1; // set to 'true' to begin with
55 xiEqual = 0; // xiEqual is 'false'
57 if(xiEqual && (omega%4 == 0)) { // if xiEqual is true and omega in {0, 4}
59 eta = calloc(*k, sizeof(int));
63 eta[a] = eta[a] + (*J)[a][b]*xi[b];
64 eta[a] = (eta[a]/4)%2;
67 gamma = calloc(n, sizeof(int));
71 gamma[a] = gamma[a] + eta[b]*GBar[b][a];
72 gamma[a] = gamma[a]%2;
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
79 free(xi); free(zeta); free(Xprime);
80 return(0.0); // Gamma = 0.0
81 } else if(shrinkResult == 'A') {//same
84 free(xi); free(zeta); free(Xprime);
85 return(1.0); // Gamma = 1.0
89 free(xi); free(zeta); free(Xprime);
90 return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
92 } else if(xiEqual && (omega+2)%4 == 0) { // if xiEqual is true and omega in {2, 6}
94 eta = calloc(*k, sizeof(int));
98 eta[a] = eta[a] + (*J)[a][b]*xi[b];
99 eta[a] = (eta[a]/4)%2;
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++) {
107 (*J)[a][b] = ((*J)[a][b] + 4*eta[a]*eta[b])%8;
109 (*J)[a][a] = (2*((*D)[a]))%8;
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++) {
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;
126 *D = calloc(*k, sizeof(int));
127 for(a=0; a<*k; a++) {
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];
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];
144 tmpMatrix[*k-1][*k-1] = (4*m)%8;
146 deallocate_mem(J, *k-1);
147 //for(a=0; a<(*k-1); a++)
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];
157 deallocate_mem(&tmpMatrix, *k);
159 free(xi); free(zeta); free(Xprime);
160 return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)