--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrink.h"
+#include "extend.h"
+#include "measurepauli.h"
+
+
+/****************************************************************************
+ * Takes a stabilizer state (n, k, h, G, GBar, Q, D, J) and a Pauli operator
+ * (m, X, Z) and returns the norm (abs. value) of the projected state.
+ * If it is non-zero,then it transforms the stabilizer state to the
+ * projected state.
+ ****************************************************************************/
+// This assumes the order Z-first X-second order when operator is read from left-to-right!
+double measurepauli(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int m, int *X, int *Z) {
+ // 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.
+ // 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.
+
+ int a, b;
+
+ int omega;
+
+ int *xi, *zeta, *Xprime;
+ xi = calloc(*k, sizeof(int));
+ zeta = calloc(*k, sizeof(int));
+ Xprime = calloc(n, sizeof(int));
+
+ for(a=0; a<*k; a++) {
+ xi[a] = dotProductMod(GBar[a], X, n, 2);
+ zeta[a] = dotProductMod(G[a], Z, n, 2);
+ }
+
+ for(a=0; a<n; a++) {
+ Xprime[a] = 0;
+ for(b=0; b<*k; b++)
+ Xprime[a] = Xprime[a] + xi[b]*G[b][a];
+ Xprime[a] = Xprime[a]%2;
+ }
+
+ omega = 2*m + 4*dotProductMod(Z, h, n, 2);
+ for(a=0; a<*k; a++) {
+ omega += (*D)[a]*xi[a];
+ for(b=a+1; b<*k; b++)
+ omega += (*J)[a][b]*xi[a]*xi[b];
+ }
+ omega = omega%8;
+
+ unsigned int xiEqual = 1; // set to 'true' to begin with
+ for(a=0; a<n; a++)
+ if(Xprime[a] != X[a])
+ xiEqual = 0; // xiEqual is 'false'
+
+ if(xiEqual && (omega%4 == 0)) { // if xiEqual is true and omega in {0, 4}
+ int* eta;
+ eta = calloc(*k, sizeof(int));
+ for(a=0; a<*k; a++) {
+ eta[a] = 4*zeta[a];
+ for(b=0; b<*k; b++)
+ eta[a] = eta[a] + (*J)[a][b]*xi[b];
+ eta[a] = (eta[a]/4)%2;
+ }
+ int* gamma;
+ gamma = calloc(n, sizeof(int));
+ for(a=0; a<n; a++) {
+ gamma[a] = 0;
+ for(b=0; b<*k; b++)
+ gamma[a] = gamma[a] + eta[b]*GBar[b][a];
+ gamma[a] = gamma[a]%2;
+ }
+ int alpha = (omega/4 + dotProductMod(gamma, h, n, 2))%2;
+ char shrinkResult = shrink(n, k, h, G, GBar, Q, D, J, gamma, alpha);
+ if(shrinkResult == 'E') {// empty
+ free(gamma);
+ free(eta);
+ free(xi); free(zeta); free(Xprime);
+ return(0.0); // Gamma = 0.0
+ } else if(shrinkResult == 'A') {//same
+ free(gamma);
+ free(eta);
+ free(xi); free(zeta); free(Xprime);
+ return(1.0); // Gamma = 1.0
+ } else {// success
+ free(gamma);
+ free(eta);
+ free(xi); free(zeta); free(Xprime);
+ return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
+ }
+ } else if(xiEqual && (omega+2)%4 == 0) { // if xiEqual is true and omega in {2, 6}
+ int* eta;
+ eta = calloc(*k, sizeof(int));
+ for(a=0; a<*k; a++) {
+ eta[a] = 4*zeta[a];
+ for(b=0; b<*k; b++)
+ eta[a] = eta[a] + (*J)[a][b]*xi[b];
+ eta[a] = (eta[a]/4)%2;
+ }
+ int sigma = 2 - omega/2;
+ *Q = (*Q + sigma)%8 < 0 ? (*Q + sigma)%8 + 8 : (*Q + sigma)%8;
+ for(a=0; a<*k; a++) {
+ (*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;
+ for(b=0; b<*k; b++) {
+ if(a != b)
+ (*J)[a][b] = ((*J)[a][b] + 4*eta[a]*eta[b])%8;
+ else
+ (*J)[a][a] = (2*((*D)[a]))%8;
+ }
+ }
+ free(eta); free(xi); free(zeta); free(Xprime);
+ return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
+ } else { // if xiEqual is false
+ extend(n, k, h, G, GBar, X);
+ int *newD = calloc(*k, sizeof(int));
+ for(a=0; a<(*k-1); a++) {
+ newD[a] = (*D)[a];
+ }
+ int *tmpVec = calloc(n, sizeof(int));
+ addVectorMod(h,X,tmpVec,n,2);
+ newD[*k-1] = (2*m + 4*dotProductMod(Z,tmpVec,n,2))%8;
+ free(tmpVec);
+ if(*k-1 > 0)
+ free(*D);
+ *D = calloc(*k, sizeof(int));
+ for(a=0; a<*k; a++) {
+ (*D)[a] = newD[a];
+ }
+ free(newD);
+
+ int **tmpMatrix = calloc(*k, sizeof(int*));
+ for(a=0; a<(*k-1); a++) {
+ tmpMatrix[a] = calloc(*k, sizeof(int));
+ for(b=0; b<(*k-1); b++) {
+ tmpMatrix[a][b] = (*J)[a][b];
+ }
+ }
+ tmpMatrix[*k-1] = calloc(*k, sizeof(int));
+ for(a=0; a<(*k-1); a++) {
+ tmpMatrix[*k-1][a] = 4*zeta[a];
+ tmpMatrix[a][*k-1] = 4*zeta[a];
+ }
+ tmpMatrix[*k-1][*k-1] = (4*m)%8;
+ if(*k-1 > 0)
+ deallocate_mem(J, *k-1);
+ //for(a=0; a<(*k-1); a++)
+ // free((*J)[a]);
+ //free(*J);
+ *J = calloc(*k, sizeof(int*));
+ for(a=0; a<*k; a++) {
+ (*J)[a] = calloc(*k, sizeof(int));
+ for(b=0; b<*k; b++) {
+ (*J)[a][b] = tmpMatrix[a][b];
+ }
+ }
+ deallocate_mem(&tmpMatrix, *k);
+
+ free(xi); free(zeta); free(Xprime);
+ return(1.0/sqrt(2.0)); // Gamma = 2^(-1/2)
+ }
+
+ return(0.0);
+
+}
projects / strong_simulation_stabilizer_rank.git / blobdiff