--- /dev/null
+BSD three-clause license
+
+Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+ * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--- /dev/null
+#IDIR =../include
+CC=gcc -std=c99
+CFLAGS=-Wall
+LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.so libinnerproduct.so
+
+weaksim_relerr: weaksim_relerr.c matrix exponentialsum shrink shrinkstar extend measurepauli innerproduct randomstabilizerstate supplement sparsify
+ $(CC) -o $@ weaksim_relerr.c $(CFLAGS) $(LIBS) libshrink.so libshrinkstar.so librandomstabilizerstate.so libsparsify.so
+
+weaksim: weaksim.c matrix exponentialsum shrink extend measurepauli innerproduct sparsify
+ $(CC) -o $@ weaksim.c $(CFLAGS) $(LIBS) libshrink.so libsparsify.so
+
+module_sparsify_test: module_sparsify_test matrix sparsify
+ $(CC) -o $@ module_sparsify_test.c $(CFLAGS) libmatrix.so libsparsify.so
+
+matrix: matrix.h matrix.c
+ $(CC) -c -Werror -Wall -fpic matrix.c
+ $(CC) -shared -o libmatrix.so matrix.o
+
+exponentialsum: exponentialsum.h exponentialsum.c
+ $(CC) -c -Wall -fpic exponentialsum.c
+ $(CC) -shared -o libexponentialsum.so exponentialsum.o -lm libmatrix.so
+
+shrink: shrink.h shrink.c
+ $(CC) -c -Wall -fpic shrink.c
+ $(CC) -shared -o libshrink.so shrink.o -lm libmatrix.so
+
+extend: extend.h extend.c
+ $(CC) -c -Werror -Wall -fpic extend.c
+ $(CC) -shared -o libextend.so extend.o -lm libmatrix.so
+
+measurepauli: measurepauli.h measurepauli.c
+ $(CC) -c -Wall -fpic measurepauli.c
+ $(CC) -shared -o libmeasurepauli.so measurepauli.o -lm libextend.so libshrink.so libmatrix.so
+
+innerproduct: innerproduct.h innerproduct.c
+ $(CC) -c -Wall -fpic innerproduct.c
+ $(CC) -shared -o libinnerproduct.so innerproduct.o -lm libextend.so libshrink.so libexponentialsum.so libmatrix.so
+
+innerproductintersection: innerproductintersection.h innerproductintersection.c
+ $(CC) -c -Wall -fpic innerproductintersection.c
+ $(CC) -shared -o libinnerproductintersection.so innerproductintersection.o -lm libextend.so libshrink.so libexponentialsum.so libmatrix.so
+
+shrinkstar: shrinkstar.h shrinkstar.c
+ $(CC) -c -Wall -fpic shrinkstar.c
+ $(CC) -shared -o libshrinkstar.so shrinkstar.o -lm libmatrix.so
+
+randomstabilizerstate: randomstabilizerstate.h randomstabilizerstate.c
+ $(CC) -c -Wall -fpic randomstabilizerstate.c
+ $(CC) -shared -o librandomstabilizerstate.so randomstabilizerstate.o -lm libmatrix.so libshrinkstar.so -llapacke
+
+supplement: supplement.h supplement.c
+ $(CC) -c -Wall -fpic supplement.c
+ $(CC) -shared -o libsupplement.so supplement.o -lm
+
+sparsify: sparsify.h sparsify.c supplement
+ $(CC) -c -Wall -fpic sparsify.c
+ $(CC) -shared -o libsparsify.so sparsify.o -lm libmatrix.so libsupplement.so
+
+randominputcommutingHermitianPauli: randominputcommutingHermitianPauli.c
+ $(CC) -o randominputcommutingHermitianPauli randominputcommutingHermitianPauli.c
+
+randominputcommutingHermitianPauli2: randominputcommutingHermitianPauli2.c
+ $(CC) -o randominputcommutingHermitianPauli2 randominputcommutingHermitianPauli2.c
+
+randominputPauli: randominputPauli.c
+ $(CC) -o randominputPauli randominputPauli.c
+
+multipauli: multipauli.c
+ $(CC) -o multipauli multipauli.c -lm
+
+.PHONY: clean
+
+clean:
+ rm ./weaksim_relerr ./weaksim ./matrix.o ./libmatrix.so ./exponentialsum.o ./libexponentialsum.so ./shrink.o ./libshrink.so ./extend.o ./libextend.so ./measurepauli.o ./libmeasurepauli.so ./libinnerproduct.so
+
+# you might want to update LD_LIBRARY_PATH to see the library:
+# export LD_LIBRARY_PATH=$PWD:$LD_LIBRARY_PATH
+# or if you have root privileges put the library /usr/local/lib or whatever library directory in your path
+# Then, use ldconfig to write the path in the config file:
+# sudo echo "/usr/local/lib" >> /etc/ld.so.conf
+# sudo ldconfig
+
+
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "exponentialsum.h"
+
+// linked list
+struct Node {
+ int data;
+ struct Node* next;
+};
+
+void printLinkedList(struct Node *node);
+void append(struct Node** head_ref, int new_data);
+void append2(struct Node** head_ref, struct Node** last_ref, int new_data);
+void freeList(struct Node** head_ref);
+
+complex double exponentialsum(int *k, int *Q, int *D, int **J) {
+
+ int a, b, c, d;
+
+ int tmp;
+
+ struct Node* setS = NULL;
+ struct Node* setWalker = NULL;
+ struct Node* setWalker2 = NULL;
+
+ struct Node* setDa = NULL; // first element in \mathcal D
+ struct Node* setDb = NULL; // second element in \mathcal D (sets should have the same size)
+
+ int setSleftover = -1; // -1 means nothing left over
+
+ int** R; // basis change matrix R and R^T
+ int** RT;
+ R = calloc(*k, sizeof(int*));
+ RT = calloc(*k, sizeof(int*));
+ // initially set it to the k-dimensional identity matrix
+ for(a=0; a<*k; a++) {
+ R[a] = calloc(*k, sizeof(int));
+ R[a][a] = 1;
+ RT[a] = calloc(*k, sizeof(int));
+ RT[a][a] = 1;
+ }
+
+ int* tmpVec; // temporary vector used in Eq. 49 to update D
+ int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
+ tmpMatrix = calloc(*k, sizeof(int*));
+ for(a=0; a<*k; a++)
+ tmpMatrix[a] = calloc(*k, sizeof(int));
+
+ for(a=0; a<*k; a++) {
+ if(D[a] == 2 || D[a] == 6) {
+ append(&setS, a);
+ }
+ }
+
+ // Choose 'a' to be the first element in the linked list setS; a=setS->data
+
+ if(setS != NULL) {
+ setWalker = setS;
+ while(setWalker->next != NULL) {
+ setWalker = setWalker->next;
+ R[setWalker->data][setS->data] += 1;
+ }
+
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;
+ /* //tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) // don't confuse this for-loop a with 'a'=setS->data */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*D[a]; */
+ /* for(a=0;a<*k;a++) */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + J[a][b]*R[c][a]*R[c][b]; */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ if(c!=setS->data)
+ tmp += R[c][c]*D[c];
+ tmp += R[c][setS->data]*D[setS->data];
+ if(setS->data>c)
+ tmp += J[c][setS->data]*R[c][c]*R[c][setS->data];
+ else if(setS->data<c)
+ tmp += J[setS->data][c]*R[c][setS->data]*R[c][c];
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=setS->data) {
+ if(d!=setS->data) {
+ tmp += R[c][c]*J[c][d]*R[d][d];
+ tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
+ tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
+ }else {
+ tmp += R[c][c]*J[c][setS->data]*R[d][setS->data];
+ }
+ } else {
+ if(d!=setS->data) {
+ tmp += R[c][setS->data]*J[setS->data][d]*R[d][d];
+ }
+ }
+ tmp += R[c][setS->data]*J[setS->data][setS->data]*R[d][setS->data];
+ tmp %= 8;
+ tmpMatrix[c][d] = tmp;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(J[c], tmpMatrix[c], sizeof(int)*(*k));
+ /* transp(R,RT,*k,*k); */
+ /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
+
+ // remove everything but first element from setS
+ //tmp = setS->data;
+ //freeList(&setS);
+ //append(&setS, tmp);
+ //setS->next = NULL;
+ setSleftover = setS->data;
+ }
+
+ struct Node* setE = NULL;
+
+ for(a=0;a<*k;a++) {
+ //if((setS != NULL && a != setS->data) || setS == NULL)
+ if((setSleftover != -1 && a != setSleftover) || setSleftover == -1)
+ append(&setE, a);
+ }
+
+ struct Node* setM = NULL;
+ int setMCount = 0; // number of elements in M
+
+ int r = 0;
+
+ while(setE != NULL) {
+ // let 'a' be setE->data (the first element in setE)
+ struct Node* setK = NULL;
+ setWalker = setE;
+ while(setWalker->next != NULL) {
+ setWalker = setWalker->next;
+ if(J[setE->data][setWalker->data] == 4)
+ append(&setK, setWalker->data);
+ }
+
+ if(setK == NULL) { // Found new monomer {'a'}
+ append(&setM, setE->data);
+ setMCount++;
+
+ if(setE->next != NULL) {
+ setE->data = setE->next->data;
+ setE->next = setE->next->next;
+ }
+ else setE = NULL;
+ } else {
+ // let 'b' be setK->data (the first element in setK)
+ setWalker = setE->next;
+ // reset 'R' to the k-dimensional identity matrix
+ for(a=0; a<*k; a++) {
+ memset(R[a], 0, sizeof(int)*(*k));
+ //for(b=0; b<*k; b++)
+ // R[a][b] = 0;
+ R[a][a] = 1;
+ }
+
+ while(setWalker != NULL) {
+ if(setWalker->data == setK->data) {
+ setWalker = setWalker->next;
+ continue; // 'c' \in E minus {'a', 'b'} where 'a'=setE->data, 'b'=setK->data and 'c'=setWalker->data
+ }
+ R[setWalker->data][setK->data] = J[setE->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
+ R[setWalker->data][setE->data] = J[setK->data][setWalker->data]/4; // divide by 4 to get binary matrix \bs J from J
+ setWalker = setWalker->next;
+ }
+ // update D and J
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;
+ //tmpVec[c] = 0;
+ /* for(a=0;a<*k;a++) {// don't confuse this for-loop a with 'a'=setS->data */
+ /* tmp+= R[c][a]*D[a]; */
+ /* //for(a=0;a<*k;a++) */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmp += J[a][b]*R[c][a]*R[c][b]; */
+ /* } */
+ if(c!=setK->data && c!=setE->data)
+ tmp += R[c][c]*D[c];
+ tmp += R[c][setK->data]*D[setK->data];
+ tmp += R[c][setE->data]*D[setE->data];
+ if(setK->data>setE->data)
+ tmp += J[setE->data][setK->data]*R[c][setE->data]*R[c][setK->data];
+ else if(setE->data>setK->data)
+ tmp += J[setK->data][setE->data]*R[c][setK->data]*R[c][setE->data];
+ if(setK->data>c)
+ tmp += J[c][setK->data]*R[c][c]*R[c][setK->data];
+ else if(setK->data<c)
+ tmp += J[setK->data][c]*R[c][setK->data]*R[c][c];
+ if(setE->data>c)
+ tmp += J[c][setE->data]*R[c][c]*R[c][setE->data];
+ else if(setE->data<c)
+ tmp += J[setE->data][c]*R[c][setE->data]*R[c][c];
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=setK->data && c!=setE->data) {
+ if(d!=setK->data && d!=setE->data) {
+ tmp += R[c][c]*J[c][d]*R[d][d];
+ tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
+ tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
+ tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
+ tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
+ }else {
+ tmp += R[c][c]*J[c][setK->data]*R[d][setK->data];
+ tmp += R[c][c]*J[c][setE->data]*R[d][setE->data];
+ }
+ } else {
+ if(d!=setK->data && d!=setE->data) {
+ tmp += R[c][setK->data]*J[setK->data][d]*R[d][d];
+ tmp += R[c][setE->data]*J[setE->data][d]*R[d][d];
+ }
+ }
+ tmp += R[c][setK->data]*J[setK->data][setK->data]*R[d][setK->data];
+ tmp += R[c][setE->data]*J[setE->data][setK->data]*R[d][setK->data];
+ tmp += R[c][setK->data]*J[setK->data][setE->data]*R[d][setE->data];
+ tmp += R[c][setE->data]*J[setE->data][setE->data]*R[d][setE->data];
+ tmp %= 8;
+ tmpMatrix[c][d] = tmp;
+ }
+ }
+ J = tmpMatrix;
+ /* transp(R,RT,*k,*k); */
+ /* multMatrixMod(J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,J,*k,*k,*k,*k,8); */
+
+ // Now {'a','b'} form a dimer
+ r++;
+ append(&setDa,setE->data);
+ append(&setDb,setK->data);
+
+ // E = E minus {'a', 'b'}
+ setWalker = setE; // 'a'=setE->data
+ while(setWalker != NULL) {
+ if(setWalker->next->data == setK->data) {
+ setWalker2 = setWalker->next;
+ setWalker->next = setWalker->next->next; // delete 'b'
+ free(setWalker2);
+ break;
+ }
+ setWalker = setWalker->next;
+ }
+ setWalker = setE;
+ setE = setE->next; // delete 'a'
+ free(setWalker);
+ }
+
+ freeList(&setK);
+ }
+
+ complex double W = 0.0 + 0.0*I;
+
+//if(setS == NULL) {
+ if(setSleftover == -1) {
+ W = cexp(I*PI*0.25*(*Q));
+ setWalker = setM;
+ for(a=0; a<setMCount; a++) {
+ W *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
+ setWalker = setWalker->next;
+ }
+ setWalker = setDa;
+ setWalker2 = setDb;
+ for(a=0; a<r; a++) {
+ W *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data])) + cexp(I*PI*0.25*(D[setWalker2->data])) - cexp(I*PI*0.25*(D[setWalker->data] + D[setWalker2->data])));
+ setWalker = setWalker->next;
+ setWalker2 = setWalker2->next;
+ }
+ } else {
+ complex double W0 = 0.0 + 0.0*I;
+ complex double W1 = 0.0 + 0.0*I;
+
+ W0 = cexp(I*PI*0.25*(*Q));
+ setWalker = setM;
+ for(a=0; a<setMCount; a++) {
+ W0 *= (1.0 + cexp(I*PI*0.25*D[setWalker->data]));
+ setWalker = setWalker->next;
+ }
+ setWalker = setDa;
+ setWalker2 = setDb;
+ for(a=0; a<r; a++) {
+ W0 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data])) + cexp(I*PI*0.25*(D[setWalker2->data])) - cexp(I*PI*0.25*(D[setWalker->data] + D[setWalker2->data])));
+ setWalker = setWalker->next;
+ setWalker2 = setWalker2->next;
+ }
+ // W1 = cexp(I*PI*0.25*(*Q+D[setS->data]));
+ W1 = cexp(I*PI*0.25*(*Q+D[setSleftover]));
+ setWalker = setM;
+ for(a=0; a<setMCount; a++) {
+ // W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setS->data])));
+ W1 *= (1.0 + cexp(I*PI*0.25*(D[setWalker->data] + J[setWalker->data][setSleftover])));
+ setWalker = setWalker->next;
+ }
+ setWalker = setDa;
+ setWalker2 = setDb;
+ for(a=0; a<r; a++) {
+ // W1 *= (1.0 + cexp(I*PI*0.25*(J[setWalker->data][setS->data] + D[setWalker->data])) + cexp(I*PI*0.25*(J[setWalker2->data][setS->data] + D[setWalker2->data])) - cexp(I*PI*0.25*(J[setWalker->data][setS->data] + J[setWalker2->data][setS->data] + D[setWalker->data] + D[setWalker2->data])));
+ W1 *= (1.0 + cexp(I*PI*0.25*(J[setWalker->data][setSleftover] + D[setWalker->data])) + cexp(I*PI*0.25*(J[setWalker2->data][setSleftover] + D[setWalker2->data])) - cexp(I*PI*0.25*(J[setWalker->data][setSleftover] + J[setWalker2->data][setSleftover] + D[setWalker->data] + D[setWalker2->data])));
+ setWalker = setWalker->next;
+ setWalker2 = setWalker2->next;
+ }
+ W = W0 + W1;
+ }
+
+ deallocate_mem(&tmpMatrix, *k);
+ deallocate_mem(&R, *k);
+ deallocate_mem(&RT, *k);
+
+ freeList(&setDa); freeList(&setDb); freeList(&setM); freeList(&setE); freeList(&setS);
+
+ return(W);
+
+}
+
+
+void append(struct Node** head_ref, int new_data)
+{
+ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
+ struct Node *last = *head_ref;
+
+ new_node->data = new_data;
+ new_node->next = NULL;
+
+ // if the linked list is empty, then make the new node as head
+ if (*head_ref == NULL) {
+ *head_ref = new_node;
+ return;
+ }
+
+ // else traverse till the last node
+ while (last->next != NULL) {
+ last = last->next;
+ }
+ last->next = new_node;
+
+ return;
+}
+
+/*void freeList(struct Node** head_ref)
+{
+
+ struct Node *second_last_walker = *head_ref;
+ struct Node *walker = *head_ref;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ while(walker->next != NULL) {
+ while(walker->next != NULL) {
+ //printf("!%d\n", walker->data);
+ second_last_walker = walker;
+ walker = walker->next;
+ //printf("%d\n", walker->data);
+ //printf("!%d\n", second_last_walker->data);
+ }
+ free(walker);
+ second_last_walker->next = NULL;
+ walker = *head_ref;
+ //printf("!!%d\n", second_last_walker->data);
+ }
+ }
+ free(walker);
+
+ return;
+ }*/
+
+void freeList(struct Node** head_ref)
+{
+
+ struct Node *walker = *head_ref;
+ struct Node *walker2;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ walker2 = walker->next;
+ while(walker2 != NULL) {
+ free(walker);
+ walker = walker2;
+ walker2 = walker->next;
+ }
+ }
+ free(walker);
+
+ return;
+}
+
+void printLinkedList(struct Node *node)
+{
+ if(node == NULL)
+ printf("NULL\n");
+ else {
+ while (node != NULL)
+ {
+ printf(" %d ", node->data);
+ node = node->next;
+ }
+ }
+ printf("\n");
+}
+
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "extend.h"
+
+// linked list
+struct Node {
+ int data;
+ struct Node* next;
+};
+
+void printLinkedList(struct Node *node);
+void append(struct Node** head_ref, int new_data);
+void freeList(struct Node** head_ref);
+
+/****************************************************************************
+ * Attempts to extend the dimension of the affine space by one dimension.
+ * The affine space is specified by
+ * the first 'k' columns of the 'n'x'n' matrix 'G', translated by 'h',
+ * that will be extended to include the vector 'xi'. 'GBar' is 'G^(-1)^T'.
+ ****************************************************************************/
+
+void extend(int n, int *k, int *h, int **G, int **GBar, int *xi) {
+
+ int a;
+
+ struct Node* setS = NULL;
+ struct Node* setT = NULL;
+ struct Node* setWalker = NULL;
+
+ int i;
+
+ for(a=0; a<*k; a++) {
+ if(dotProductMod(xi, GBar[a], n, 2) == 1) {
+ append(&setS, a);
+ }
+ }
+ unsigned int firstT = 1; // set initially to true
+ for(a=*k; a<n; a++) {
+ if(dotProductMod(xi, GBar[a], n, 2) == 1) {
+ if(firstT) {
+ append(&setT, a);
+ i = a;
+ firstT = 0;
+ }
+ else {
+ append(&setS, a);
+ append(&setT, a);
+ }
+ }
+ }
+
+ if(setT == NULL) {
+ freeList(&setS);
+ freeList(&setT);
+ return; // xi is in the linear space so do nothing
+ }
+
+ setWalker = setS;
+ while(setWalker != NULL) {
+ // update GBar rows, gbar
+ for(a=0; a<n; a++)
+ GBar[setWalker->data][a] = (GBar[setWalker->data][a] + GBar[i][a])%2;
+ // update G rows, g
+ for(a=0; a<n; a++)
+ G[i][a] = (G[i][a] + G[setWalker->data][a])%2;
+
+ setWalker = setWalker->next;
+ }
+
+ // Swap 'i' and 'k+1'
+ int* tmpVec;
+ tmpVec = G[i];
+ G[i] = G[*k];
+ G[*k] = tmpVec;
+ tmpVec = GBar[i];
+ GBar[i] = GBar[*k];
+ GBar[*k] = tmpVec;
+
+ freeList(&setS);
+ freeList(&setT);
+
+ // increment 'k'
+ *k = *k + 1;
+
+ return;
+
+}
+
+
+void append(struct Node** head_ref, int new_data)
+{
+ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
+ struct Node *last = *head_ref;
+
+ new_node->data = new_data;
+ new_node->next = NULL;
+
+ // if the linked list is empty, then make the new node as head
+ if (*head_ref == NULL) {
+ *head_ref = new_node;
+ return;
+ }
+
+ // else traverse till the last node
+ while (last->next != NULL) {
+ last = last->next;
+ }
+ last->next = new_node;
+
+ return;
+}
+
+/*void freeList(struct Node** head_ref)
+{
+
+ struct Node *second_last_walker = *head_ref;
+ struct Node *walker = *head_ref;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ while(walker->next != NULL) {
+ while(walker->next != NULL) {
+ //printf("!%d\n", walker->data);
+ second_last_walker = walker;
+ walker = walker->next;
+ //printf("%d\n", walker->data);
+ //printf("!%d\n", second_last_walker->data);
+ }
+ free(walker);
+ second_last_walker->next = NULL;
+ walker = *head_ref;
+ //printf("!!%d\n", second_last_walker->data);
+ }
+ }
+ free(walker);
+
+ return;
+ }*/
+
+void freeList(struct Node** head_ref)
+{
+
+ struct Node *walker = *head_ref;
+ struct Node *walker2;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ walker2 = walker->next;
+ while(walker2 != NULL) {
+ free(walker);
+ walker = walker2;
+ walker2 = walker->next;
+ }
+ }
+ free(walker);
+
+ return;
+}
+
+
+void printLinkedList(struct Node *node)
+{
+ if(node == NULL)
+ printf("NULL\n");
+ else {
+ while (node != NULL)
+ {
+ printf(" %d ", node->data);
+ node = node->next;
+ }
+ }
+ printf("\n");
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrink.h"
+#include "extend.h"
+#include "exponentialsum.h"
+#include "innerproduct.h"
+
+
+/****************************************************************************
+ * Note: Arguments are copied and not modified!
+ ****************************************************************************/
+
+double complex innerproduct(int n1, int k1, int *h1, int **G1, int **GBar1, int Q1, int *D1, int **J1, int n2, int k2, int *h2, int **G2, int **GBar2, int Q2, int *D2, int **J2) {
+ // n1 should equal n2 (we assume this)
+
+ int a, b, c;
+
+ int n = n1;
+ int k = k1;
+ int *h = calloc(n, sizeof(int));
+ int **G = calloc(n, sizeof(int*));
+ int **GBar = calloc(n, sizeof(int*));
+ for(a=0; a<n; a++) {
+ h[a] = h1[a];
+ G[a] = calloc(n, sizeof(int));
+ GBar[a] = calloc(n, sizeof(int));
+ for(b=0; b<n; b++) {
+ G[a][b] = G1[a][b];
+ GBar[a][b] = GBar1[a][b];
+ }
+ }
+
+ int Q1copy = Q1;
+ int *D1copy = calloc(k1, sizeof(int));
+ int **J1copy = calloc(k1, sizeof(int*));
+ for(a=0; a<k1; a++) {
+ D1copy[a] = D1[a];
+ J1copy[a] = calloc(k1, sizeof(int));
+ for(b=0; b<k1; b++) {
+ J1copy[a][b] = J1[a][b];
+ }
+ }
+
+ int Q2copy = Q2;
+ int *D2copy = calloc(k2, sizeof(int));
+ int **J2copy = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++) {
+ D2copy[a] = D2[a];
+ J2copy[a] = calloc(k2, sizeof(int));
+ for(b=0; b<k2; b++) {
+ J2copy[a][b] = J2[a][b];
+ }
+ }
+
+ int alpha;
+ int **R = calloc(k, sizeof(int*));
+ for(a=0; a<k; a++)
+ R[a] = calloc(k2, sizeof(int));
+ int **RT = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ RT[a] = calloc(k, sizeof(int));
+ char shrinkResult;
+
+ for(a=k2; a<n2; a++) {
+ alpha = dotProductMod(h2, GBar2[a], n2, 2);
+
+ shrinkResult = shrink(n, &k, h, G, GBar, &Q1copy, &D1copy, &J1copy, GBar2[a], alpha);
+ if(shrinkResult == 'E') {// if EMPTY
+ free(h);
+ deallocate_mem(&G, n); deallocate_mem(&GBar, n);
+ free(D1copy);
+ deallocate_mem(&J1copy, k);
+ free(D2copy);
+ deallocate_mem(&J2copy, k2);
+ deallocate_mem(&R, k);
+ deallocate_mem(&RT, k2);
+ return(0.0+0.0*I);
+ }
+ }
+
+ // Now K = K1 intersect K2 = (n1, k1, h1, G1, GBar1)
+ int *y = calloc(k2, sizeof(int));
+ int *tmpVec = calloc(n2, sizeof(int));
+ addVectorMod(h,h2,tmpVec,n2,2);
+ for(a=0; a<k2; a++) {
+ y[a] = dotProductMod(tmpVec,GBar2[a],n2,2);
+ for(b=0; b<k; b++) {
+ R[b][a] = dotProductMod(G[b],GBar2[a],n2,2);
+ }
+ }
+
+ // IF YOU UNCOMMENT THIS THEN YOU WILL BE CHANGING h2!!!
+ /* for(a=0; a<k2; a++) { */
+ /* for(b=0; b<n2; b++) { */
+ /* h2[b] = (h2[b] + y[a]*G2[a][b])%2; */
+ /* } */
+ /* } */
+ /* printf("h2 (should now match h):\n"); */
+ /* printVector(h2, n2); */
+ /* printf("h:\n"); */
+ /* printVector(h, n2); */
+
+ free(tmpVec);
+
+ // update Q2 and D2 using Eqs. 52 and 53 with y
+ for(a=0; a<k2; a++) {
+ Q2copy = Q2copy + D2copy[a]*y[a];
+ for(b=a+1; b<k2; b++) {
+ Q2copy = Q2copy + J2copy[a][b]*y[a]*y[b];
+ }
+ }
+ Q2copy = Q2copy%8;
+
+ for(a=0; a<k2; a++) {
+ for(b=0; b<k2; b++)
+ D2copy[a] = D2copy[a] + J2copy[a][b]*y[b];
+ D2copy[a] = D2copy[a]%8;
+ }
+
+ free(y);
+
+ // update D2copy and J2copy using Eqs. 49 and 50 with R
+ // (note that R is a kxk2 matrix so this may reduce the dimension of D2copy and J2copy)
+ y = calloc(k, sizeof(int));
+ for(c=0;c<k;c++) {
+ y[c] = 0;
+ for(a=0;a<k2;a++) {
+ y[c] = y[c] + R[c][a]*D2copy[a];
+ for(b=a+1;b<k2;b++)
+ y[c] = y[c] + J2copy[a][b]*R[c][a]*R[c][b];
+ }
+ y[c] = y[c]%8;
+ }
+ free(D2copy);
+ D2copy = calloc(k, sizeof(int));
+ memcpy(D2copy, y, sizeof(int)*k);
+ free(y);
+
+ transp(R,RT,k,k2);
+ int **tmpMatrix = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ tmpMatrix[a] = calloc(k, sizeof(int));
+ multMatrixMod(J2copy,RT,tmpMatrix,k2,k2,k2,k,8);
+ deallocate_mem(&J2copy, k2);
+ J2copy = calloc(k, sizeof(int*));
+ for(a=0;a<k;a++)
+ J2copy[a] = calloc(k, sizeof(int));
+ multMatrixMod(R,tmpMatrix,J2copy,k,k2,k2,k,8);
+
+ deallocate_mem(&tmpMatrix, k2);
+
+ deallocate_mem(&R, k);
+ deallocate_mem(&RT, k2);
+
+ // Now (Q1, D1, J1) and (Q2, D2, J2) are defined in the same basis
+
+ int Q = (Q1copy - Q2copy)%8 < 0 ? (Q1copy - Q2copy)%8 + 8 : (Q1copy - Q2copy)%8;
+ int *D = calloc(k, sizeof(int));
+ int **J = calloc(k, sizeof(int*));
+
+ for(a=0; a<k; a++) { // again, if k=0 then no need to fill these
+ J[a] = calloc(k, sizeof(int));
+ D[a] = (D1copy[a] - D2copy[a])%8 < 0 ? (D1copy[a] - D2copy[a])%8 + 8 : (D1copy[a] - D2copy[a])%8;
+ for(b=0; b<k; b++)
+ J[a][b] = (J1copy[a][b] - J2copy[a][b])%8 < 0 ? (J1copy[a][b] - J2copy[a][b])%8 + 8 : (J1copy[a][b] - J2copy[a][b])%8;
+
+ }
+
+ free(h);
+ deallocate_mem(&G, n); deallocate_mem(&GBar, n);
+
+ double complex amplitude = pow(2,-0.5*(k1+k2))*exponentialsum(&k, &Q, D, J);
+
+ free(D);
+ deallocate_mem(&J, k);
+
+ free(D1copy);
+ deallocate_mem(&J1copy, k);
+
+ free(D2copy);
+ deallocate_mem(&J2copy, k);
+
+ return(amplitude);
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "matrix.h"
+
+
+int** addMatrix(int **A, int **B, int rows, int cols)
+{
+ int i, j;
+
+ int** C;
+
+ C = calloc(cols, sizeof(int*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(int));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = A[i][j] + B[i][j];
+
+ return C;
+}
+
+void addMatrixMod(int **A, int **B, int **C, int rows, int cols, int mod)
+{
+ int i, j;
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = (A[i][j] + B[i][j])%mod;
+
+ return;
+}
+
+void addVectorMod(int *A, int *B, int *C, int length, int mod)
+{
+ int i;
+
+ for(i=0; i<length; i++)
+ C[i] = (A[i] + B[i])%mod;
+
+ return;
+}
+
+void scalarmultMatrix(int scalar, int **a, int **b, int rows, int cols)
+{
+ int i, j;
+
+ //b = calloc(cols, sizeof(int*));
+ //for(i=0; i<cols; i++)
+ // b[i] = calloc(rows, sizeof(int));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ b[i][j] = scalar*a[i][j];
+
+}
+
+int trace(int **a, int rows, int cols)
+{
+ int i;
+ int tr = 0.0;
+
+ for(i=0; i<rows; i++)
+ tr += a[i][i];
+
+ return tr;
+}
+
+// b = a^T
+void transp(int **a, int **b, int rows, int cols)
+{
+ int i, j;
+
+ for(i=0; i<cols; i++) {
+ for(j=0; j<rows; j++)
+ b[i][j] = a[j][i];
+ }
+
+}
+
+void printVector(int* a, int length)
+{
+ int i;
+ printf("Vector[%d]\n", length);
+ for(i=0; i<length; i++)
+ printf("%d ", a[i]);
+ if(length>0)
+ printf("\n");
+
+}
+
+void printMatrix(int** a, int rows, int cols)
+{
+ int i, j;
+ printf("Matrix[%d][%d]\n", rows, cols);
+ for(i=0; i<rows; i++) {
+ for(j=0; j<cols; j++) {
+ printf("%d ", a[i][j]);
+ }
+ printf("\n");
+ }
+}
+
+int** multMatrix(int **A, int **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k;
+ int **C;
+ C = calloc(ro1, sizeof(int*));
+ for(i=0; i<ro1; i++)
+ C[i] = calloc(co2, sizeof(int));
+
+ for(i=0; i<ro1; i++) {
+ for(j=0; j<co2; j++) {
+ C[i][j] = 0;
+ for(k=0; k<co1; k++)
+ C[i][j] += A[i][k] * B[k][j];
+ }
+ }
+
+ return C;
+}
+
+// A times B = C
+void multMatrixMod(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2, int mod)
+{
+ int i, j, k, tmp;
+
+ for(i=0; i<ro1; i++) {
+ for(j=0; j<co2; j++) {
+ tmp = 0;
+ for(k=0; k<co1; k++)
+ tmp += (A[i][k] * B[k][j]);
+ C[i][j] = tmp%mod;
+ }
+ }
+
+}
+
+int** outerMatrix(int **A, int **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k, l;
+ int **C;
+ C = calloc(ro1*ro2, sizeof(int*));
+ for(i=0; i<ro1*ro2; i++)
+ C[i] = calloc(co1*co2, sizeof(int));
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ for(k=0; k<co1; k++)
+ for(l=0; l<co2; l++) {
+ C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
+ }
+
+ return C;
+}
+
+void outerMatrixMod(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2, int mod)
+{
+ int i, j, k, l;
+ //int **C;
+ //C = calloc(ro1*ro2, sizeof(int*));
+ //for(i=0; i<ro1*ro2; i++)
+ // C[i] = calloc(co1*co2, sizeof(int));
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ for(k=0; k<co1; k++)
+ for(l=0; l<co2; l++) {
+ C[j+ro2*i][l+co2*k] = (A[i][k]* B[j][l])%mod;
+ }
+
+}
+
+void outerVectorMod(int *A, int *B, int *C, int ro1, int ro2, int mod)
+{
+ int i, j;
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ C[j+ro2*i]= (A[i]* B[j])%mod;
+
+}
+
+void addSubMatrix(int **A, int **B, int ro1, int co1, int rooff1, int cooff1, int rooff2, int cooff2)
+{
+ // rooff1 is row offset to start range ro1 in first matrix etc.
+ int i, j;
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<co1; j++)
+ B[rooff2+i][cooff2+j] = A[rooff1+i][cooff1+j];
+
+}
+
+void appendBlockMatrix(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2)
+{
+ // rooff1 is row offset to start range ro1 in first matrix etc.
+ int i, j;
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<co1; j++)
+ C[i][j] = A[i][j];
+
+ for(i=0; i<ro2; i++)
+ for(j=0; j<co2; j++)
+ C[ro1+i][co1+j] = B[i][j];
+
+}
+
+
+void appendVector(int *A, int *B, int *C, int ro1, int ro2)
+{
+ int i;
+
+ for(i=0; i<ro1; i++)
+ C[i] = A[i];
+ for(i=0; i<ro2; i++)
+ C[ro1+i] = B[i];
+
+}
+
+
+int dotProductMod(int *a, int *b, int length, int mod)
+{
+ int i;
+ int dotproduct = 0;
+
+ for(i=0; i<length; i++)
+ dotproduct += a[i]*b[i];
+
+ return(dotproduct%mod);
+}
+
+
+void deallocate_mem(int ***arr, int rows)
+{
+ int i;
+ for(i=0; i<rows; i++)
+ free((*arr)[i]);
+ free(*arr);
+}
--- /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);
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <string.h>
+#include <math.h>
+
+#include "lapacke.h"
+#include "matrix.h"
+#include "shrinkstar.h"
+
+#include "randomstabilizerstate.h"
+
+#define ZEROTHRESHOLD (0.00000001)
+
+// Note: Make sure you seed the random number generator before calling this!
+// i.e. srand((unsigned)time(NULL));
+void randomstabilizerstate(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J, double** Pd)
+{
+ // vector and matrix pointers should all be passed unallocated! (no freeing of memory performed)
+
+ float *randomX, *randomXcopy; // d*n matrix that is in array-form for LAPACKE
+ float *randomXsingvals; // singular values of randomX
+ float *U, *VT, *superb;
+ int d;
+
+ int i, j;
+
+ //srand((unsigned)time(NULL));
+
+ double randPd = (double)rand()/(double)(RAND_MAX);
+
+ double cumprob = 0.0;
+ for(i=0;i<n+1;i++) {
+ cumprob += Pd[n-1][i]; // cumprob goes up to 1.0 once you count all Pd[n][0:n+1]
+ if(cumprob > randPd) {
+ d = i;
+ break;
+ }
+ d=i;
+ }
+ //d = rand()%(n+1);
+ //printf("d=%d total=%f randPd=%f\n", d, total, randPd);
+ //printf("!!!d=%d\n", d);
+
+ *k = n; // we will get it to be k=n-d by caling shrinkstar d times below
+
+ randomX = calloc(n*d,sizeof(float));
+ randomXcopy = calloc(n*d,sizeof(float));
+ randomXsingvals = calloc(d,sizeof(float));
+ superb = calloc((d),sizeof(float));
+ U = calloc((d*d),sizeof(float));
+ VT = calloc((n*n),sizeof(float));
+
+ int info;
+ int numsingvals = -1;
+
+ while(numsingvals != d) {
+ for(i=0; i<d; i++) {
+ for(j=0; j<n; j++) {
+ randomX[i*n+j] = (float)(rand()%2);
+ randomXcopy[i*n+j] = randomX[i*n+j];
+ //printf("%1.0f ", randomX[i*n+j]);
+ }
+ //printf("\n");
+ }
+
+ info = LAPACKE_sgesvd( LAPACK_ROW_MAJOR, 'N', 'N', d, n, randomXcopy, n, randomXsingvals, U, n, VT, n, superb );
+
+ numsingvals = 0;
+
+ for(i=0; i<d; i++) {
+ if(fabs(randomXsingvals[i]) >= ZEROTHRESHOLD)
+ numsingvals++;
+ }
+
+ //printf("Number of singular values: %d\n", numsingvals);
+ }
+
+ *G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
+
+ *h = calloc(n, sizeof(int));
+
+ for(i=0; i<n; i++) {
+ (*G)[i] = calloc(n, sizeof(int));
+ (*GBar)[i] = calloc(n, sizeof(int));
+
+ (*G)[i][i] = 1;
+ (*GBar)[i][i] = 1;
+ }
+
+ int* xi = calloc(n, sizeof(int));
+ for(i=0; i<d; i++) {
+ for(j=0; j<n; j++){
+ xi[j] = (int)randomX[i*n+j];
+ //printf("%1.0f ", randomX[i*n+j]);
+ }
+ //printf("\n");
+ //printVector(xi, n);
+ shrinkstar(n, k, *h, *G, *GBar, xi, 0);
+ }
+ free(xi);
+ //printf("n-d=%d\n", n-d);
+ //printf("k=%d\n", *k);
+ //printMatrix(*G, n, n);
+
+ for(i=0; i<n; i++)
+ (*h)[i] = rand()%2;
+
+ *Q = (rand()%8); // Q \in Z/8Z
+
+ *D = calloc(*k, sizeof(int));
+ for(i=0; i<*k; i++)
+ (*D)[i] = (rand()%4)*2; // D_a \in {0, 2, 4, 6}
+
+ *J = calloc(*k, sizeof(int*));
+ for(i=0; i<*k; i++) {
+ (*J)[i] = calloc(*k, sizeof(int));
+ for(j=(i+1); j<*k; j++)
+ (*J)[i][j] = (rand()%2)*4; // J_{a,b} \in {0, 4} for a!=b
+ (*J)[i][i] = (2*(*D)[i])%8;
+ }
+ for(i=0; i<*k; i++) {
+ for(j=(i+1); j<*k; j++) {
+ (*J)[j][i] = (*J)[i][j];
+ }
+ }
+
+ free(randomX); free(randomXcopy); free(randomXsingvals); free(superb); free(U); free(VT);
+ return;
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrink.h"
+
+// linked list
+struct Node {
+ int data;
+ struct Node* next;
+};
+
+void printLinkedList(struct Node *node);
+void append(struct Node** head_ref, int new_data);
+void freeList(struct Node** head_ref);
+
+/****************************************************************************
+ * Attempts to shrinks the dimension of the affine space by one dimension.
+ * The affine space is specified by
+ * the first 'k' columns of the 'n'x'n' matrix 'G', translated by 'h',
+ * that have an inner product with 'xi' of 'alpha'. 'GBar' is 'G^(-1)^T',
+ * and ('Q', 'D', 'J') specify the quadratic form on this affine space.
+ * ---
+ * If result is an EMPTY affine space then return value is 'E'.
+ * If result is the SAME affine space then return value is 'A'.
+ * If result is a SUCCESSfully shrunk affine space then return value is 'U'.
+ * ***
+ * Note that D, and J are assumed to be nx1 and nxn even though they are
+ * actually kx1 and kxk (extra elements are ignored).
+ ****************************************************************************/
+
+char shrink(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int *xi, int alpha) {
+ // 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, c, d;
+ int beta;
+
+ struct Node* setS = NULL;
+ struct Node* setWalker = NULL;
+
+ for(a=0; a<*k; a++) {
+ if(dotProductMod(xi, G[a], n, 2) == 1)
+ append(&setS, a);
+ }
+
+ beta = (alpha + dotProductMod(xi, h, n, 2))%2;
+
+ if(setS == NULL) {
+ if(beta == 1) {
+ freeList(&setS);
+ return 'E'; // EMPTY
+ }
+ else {// beta = 0
+ freeList(&setS);
+ return 'A'; // SAME
+ }
+ }
+
+ // we define tmpVec, tmpMatrix, R and RT now for later use since we are sure that *k!=0 (since otherwise setS==NULL) and so they will be finite-sized
+ int* tmpVec; // temporary vector used in swapping G and GBar rows and in Eq. 49 for updating D
+ int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
+ tmpMatrix = calloc(*k, sizeof(int*));
+ for(a=0; a<*k; a++)
+ tmpMatrix[a] = calloc(*k, sizeof(int));
+
+ int** R; // basis change matrix R and R^T
+ int** RT;
+ R = calloc(*k, sizeof(int*));
+ RT = calloc(*k, sizeof(int*));
+ // initially set it to the k-dimensional identity matrix
+ for(a=0; a<*k; a++) {
+ R[a] = calloc(*k, sizeof(int));
+ R[a][a] = 1;
+ RT[a] = calloc(*k, sizeof(int));
+ RT[a][a] = 1;
+ }
+
+ int i = setS->data; // pick first element in setS to be added to the new space's shift vector
+ // remove from setS
+ if(setS->next != NULL) {
+ setWalker = setS->next;
+ setS->data = setS->next->data;
+ setS->next = setS->next->next;
+ free(setWalker);
+ } else
+ setS = NULL;
+
+ int tmp;
+ setWalker = setS;
+ while(setWalker != NULL) {
+ // update G rows, g
+ for(a=0; a<n; a++)
+ G[setWalker->data][a] = (G[setWalker->data][a] + G[i][a])%2;
+
+ // update D and J
+ R[setWalker->data][i] = 1;
+ RT[i][setWalker->data] = 1;
+
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ if(c!=i)
+ tmp += R[c][c]*((*D)[c]);
+ tmp += R[c][i]*((*D)[i]);
+ if(i>c)
+ tmp += ((*J)[c][i])*R[c][c]*R[c][i];
+ else if(i<c)
+ tmp += ((*J)[i][c])*R[c][i]*R[c][c];
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=i) {
+ if(d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ }
+ } else {
+ if(d!=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }
+ }
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp %= 8;
+ tmpMatrix[c][d] = tmp;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+
+ if(setWalker->data != i) {
+ R[setWalker->data][i] = 0; // reset R
+ RT[i][setWalker->data] = 0;
+ }
+
+ // update GBar rows, gbar
+ for(a=0; a<n; a++)
+ GBar[i][a] = (GBar[i][a] + GBar[setWalker->data][a])%2;
+
+ setWalker = setWalker->next;
+ }
+
+ // Swap 'i' and 'k'
+ R[*k-1][*k-1] = 0; R[i][i] = 0;
+ R[i][*k-1] = 1; R[*k-1][i] = 1;
+ RT[*k-1][*k-1] = 0; RT[i][i] = 0;
+ RT[i][*k-1] = 1; RT[*k-1][i] = 1;
+
+ tmpVec = G[i];
+ G[i] = G[*k-1];
+ G[*k-1] = tmpVec;
+ tmpVec = GBar[i];
+ GBar[i] = GBar[*k-1];
+ GBar[*k-1] = tmpVec;
+
+ // update D and J
+ //(Note: can't quite directly use previous faster algorithm since R[c][c]!=1 forall c. However, this can still be brought down by a factor of N by the same algorithm appropriately tweaked.)
+ if(i!=(*k-1)) {
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) { */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* } */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ tmp += R[c][c]*((*D)[c]); // iff c!=i and c!=(*k-1) this can be non-zero
+ tmp += R[c][*k-1]*((*D)[*k-1]); // iff c=i this can be non-zero
+ tmp += R[c][i]*((*D)[i]); // iff c=*k-1 this can be non-zero
+
+ /* if(i>c) */
+ /* tmp += ((*J)[c][i])*R[c][c]*R[c][i]; */
+ /* else if(i<c) */
+ /* tmp += ((*J)[i][c])*R[c][i]*R[c][c]; */
+ /* if((*k-1)>c) */
+ /* tmp += ((*J)[c][*k-1])*R[c][c]*R[c][*k-1]; */
+ /* else if((*k-1)<c) // pretty sure this can't happen */
+ /* tmp += ((*J)[*k-1][c])*R[c][*k-1]*R[c][c]; */
+ /* if((*k-1)>i) */
+ /* tmp += ((*J)[i][*k-1])*R[c][i]*R[c][*k-1]; */
+ /* else if((*k-1)<i) */
+ /* tmp += ((*J)[*k-1][i])*R[c][*k-1]*R[c][i]; */
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=(*k-1) && c!=i) {
+ if(d!=(*k-1) && d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][c]*((*J)[c][*k-1])*R[d][*k-1];
+ }
+ } else {
+ if(d!=(*k-1) && d !=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ tmp += R[c][*k-1]*((*J)[*k-1][d])*R[d][d];
+ } else {
+ tmp += R[c][i]*((*J)[i][*k-1])*R[d][*k-1];
+ tmp += R[c][*k-1]*((*J)[*k-1][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp += R[c][*k-1]*((*J)[*k-1][*k-1])*R[d][*k-1];
+ }
+ }
+ tmpMatrix[c][d] = tmp%8;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+ }
+ // reset R & RT
+ R[*k-1][*k-1] = 1; R[i][i] = 1;
+ R[i][*k-1] = 0; R[*k-1][i] = 0;
+ RT[*k-1][*k-1] = 1; RT[i][i] = 1;
+ RT[i][*k-1] = 0; RT[*k-1][i] = 0;
+
+ for(a=0; a<n; a++)
+ h[a] = (h[a] + beta*G[*k-1][a])%2;
+
+ // update Q & D using Eqs. 52 & 53 with y = beta*G[*k-1][:]
+ // since h is expressed in F_2^n, in F_2^n, y_i = beta*delta_{i,*k-1}
+ // (and so the second sum of Eq. 52 is always over zero terms)
+ *Q = (*Q + ((*D)[*k-1])*beta)%8;
+
+ for(a=0; a<*k; a++) {
+ (*D)[a] = ((*D)[a] + (*J)[a][*k-1]*beta)%8;
+ }
+
+ // remove the kth row & column from J
+ // remove the ith element from D
+
+ // we only need tmpMatrix to be (k-1)x(k-1) and it is kxk
+ for(a=0; a<(*k-1); a++) {
+ //tmpMatrix[a] = calloc(*k-1, sizeof(int)); // no need!
+ for(b=0; b<(*k-1); b++)
+ tmpMatrix[a][b] = (*J)[a][b];
+ }
+ deallocate_mem(J, *k);
+ *J = calloc(*k-1, sizeof(int*));
+ for(a=0; a<(*k-1); a++) {
+ (*J)[a] = calloc(*k-1, sizeof(int));
+ for(b=0; b<(*k-1); b++)
+ (*J)[a][b] = tmpMatrix[a][b];
+ }
+ deallocate_mem(&tmpMatrix, *k);
+
+ tmpVec = calloc(*k-1, sizeof(int));
+ for(a=0; a<(*k-1); a++)
+ tmpVec[a] = (*D)[a];
+ free(*D);
+ *D = calloc(*k-1, sizeof(int));
+ memcpy(*D, tmpVec, sizeof(int)*(*k-1));
+ free(tmpVec);
+
+ deallocate_mem(&R, *k);
+ deallocate_mem(&RT, *k);
+
+ // decrement 'k'
+ *k = *k - 1;
+
+ freeList(&setS);
+
+ return 'U'; // SUCCESS
+
+}
+
+
+void append(struct Node** head_ref, int new_data)
+{
+ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
+ struct Node *last = *head_ref;
+
+ new_node->data = new_data;
+ new_node->next = NULL;
+
+ // if the linked list is empty, then make the new node as head
+ if (*head_ref == NULL) {
+ *head_ref = new_node;
+ return;
+ }
+
+ // else traverse till the last node
+ while (last->next != NULL) {
+ last = last->next;
+ }
+ last->next = new_node;
+
+ return;
+}
+
+/*void freeList(struct Node** head_ref)
+{
+
+ struct Node *second_last_walker = *head_ref;
+ struct Node *walker = *head_ref;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ while(walker->next != NULL) {
+ while(walker->next != NULL) {
+ //printf("!%d\n", walker->data);
+ second_last_walker = walker;
+ walker = walker->next;
+ //printf("%d\n", walker->data);
+ //printf("!%d\n", second_last_walker->data);
+ }
+ free(walker);
+ second_last_walker->next = NULL;
+ walker = *head_ref;
+ //printf("!!%d\n", second_last_walker->data);
+ }
+ }
+ free(walker);
+
+ return;
+ }*/
+
+void freeList(struct Node** head_ref)
+{
+
+ struct Node *walker = *head_ref;
+ struct Node *walker2;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ walker2 = walker->next;
+ while(walker2 != NULL) {
+ free(walker);
+ walker = walker2;
+ walker2 = walker->next;
+ }
+ }
+ free(walker);
+
+ return;
+}
+
+void printLinkedList(struct Node *node)
+{
+ if(node == NULL)
+ printf("NULL\n");
+ else {
+ while (node != NULL)
+ {
+ printf(" %d ", node->data);
+ node = node->next;
+ }
+ }
+ printf("\n");
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "supplement.h"
+#include "sparsify.h"
+
+/****************************************************************************
+ * SPARSIFY procedure introduced by Bravyi et al.
+ * (S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and
+ * M. Howard, Simulation of Quantum Circuits by Low-Rank Sta-
+ * bilizer Decompositions, Quantum 3, 181 (2019).)
+ * Outputs stabilizer state indices in stabStateIndices
+ * as understood by weaksim.c in the tilde'd "computational" basis.
+ ***Needs address of stabStateIndices since it allocates it.***
+ * 'numStabStates' is the number of stabilizer state indices,
+ * 'T' is the number of T-gate magic states, 'phi' is the angle of the
+ * intermediate diagonal magic state, and 'delta' is the additive
+ * error of the state.
+ * For phi=pi/4 this should give you a uniform distro over [0,2^T-1].
+ ****************************************************************************/
+
+// Note: Make sure you seed the random number generator before calling this!
+// i.e. srand((unsigned)time(NULL));
+int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling) {
+
+ int i, j, k, l;
+
+ double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
+ //printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0));
+
+ printf("delta = %lf\n", delta);
+
+ //printf("L1Norm = %lf\n", L1Norm);
+
+ double gamma;
+ int increment;
+
+ if(coherentSampling) {
+ printf("Performing coherent sampling!\n");
+ gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
+ // since we will be taking a multiple of (T+1) samples, we round numStabStates up to the nearest multiple of (T+1)
+ *numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
+ increment = T+1;
+ } else {
+ //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
+ //*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
+ gamma = 1.0;
+ *numStabStates = ceil((pow(L1Norm,2.0)-gamma)/pow(delta,2.0));
+ increment = 1;
+ }
+
+
+ printf("numStabStates = %d\n", *numStabStates);
+ fflush(stdout);
+
+ double cumL1Norm, sum;
+ long index, tmpIndex;
+ int shiftCounter, digit;
+
+ *stabStateIndices = calloc(*numStabStates,sizeof(long));
+
+ for(i=0; i<*numStabStates; i=i+increment) {
+ // We will pick our stabilizer index by first selecting how many tilde'd 1s it will have
+ // by randomly sampling weighted by the binomial expansion of that in the L1norm,
+ // and then uniformly randomly picking which digits to make that many tilde 1's
+
+ // Get random number and scale by the L_1 norm of the t-qubit full stabilizer decomposition
+ cumL1Norm = ((((double)rand())/((double)(RAND_MAX)))*L1Norm); // this is the "cumulative" L1 norm of the randomly selected state if you index by the tilde'd computational basis states truncated to an integer
+ sum = 0.0;
+ //printf("%d cumL1Norm = %lf\n", i, cumL1Norm);
+ for(j=T; j>=0; j--) {
+ // add cumulative L1 norm to sum until you are greater than cumL1Norm of random state
+ //printf("binomcoeff(%d,%d)=%lf\n", T, j, binomcoeff(T,j));
+ if(j < T)
+ sum += binomcoeff(T,j)*(sqrt(1.0-cos(phi)))*pow(sqrt(1.0 - cos(0.25*PI)),T-1-j)*pow(sqrt(1.0-sin(phi)),j);
+ else
+ sum += pow((sqrt(1.0-sin(phi))),T);
+ //printf("%d sum=%lf\n", j, sum);
+ if(cumL1Norm <= sum) {
+ // the randomly chosen tilde'd computational state will have j tilde'd 1s
+ // Now uniformly randomly choose a state with j tilde'd 1s
+ //printf("%d tilde'd ones\n", (T-j));
+ index = 0;
+ for(k=0; k<(T-j); k++) {
+ digit = (int)((((double)rand())/((double)(RAND_MAX)))*(T-k));
+ //printf("digit=%d\n", digit);
+ tmpIndex = index;
+ shiftCounter = 0;
+ // first shift past leading 1's
+ while(tmpIndex%2 == 1) {
+ tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
+ shiftCounter++;
+ }
+ // now shift past 'digit' 0's
+ for(l=0; l<digit; l++){
+ tmpIndex /= 2; // shift to the left by one (base 2 arithmetic)
+ shiftCounter++;
+ while(tmpIndex%2 == 1) {
+ tmpIndex /= 2; // keep shifting to the left past any 1's that have already been assigned
+ shiftCounter++;
+ }
+ }
+ //printf("shiftCounter = %d\n", shiftCounter);
+ index += pow(2,shiftCounter);
+ }
+ (*stabStateIndices)[i] = index;
+ //printf("first index = %d\n", index);
+ if(coherentSampling)
+ for(k=1; k<=T; k++)
+ //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
+ supplement(&(*stabStateIndices)[i], T);
+ break;
+ }
+ //printf("%lf\n", sum);
+ }
+ //printf("index:%ld\n ", index);
+ //printf("(*stabStateIndices)[%d]:%ld\n ", i, (*stabStateIndices)[i]);
+ }
+ //printf("\n");
+
+ return 0;
+
+}
+
+// binomial coefficient
+double binomcoeff(int m, int n) {
+ //((double)fac(T))/((double)(fac(j)*fac(T-j)))
+ int i;
+ double coeff = 1.0;
+ if (n<=m) {
+ for(i=n+1;i<=m;i++) {
+ coeff *= (double)i/((double)i-n);
+ //printf("%lf ", coeff);
+ }
+ //printf("\n");
+ return coeff;
+ }
+ else
+ return 1;
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+
+void printBinary(unsigned long number, int T) {
+
+ int i;
+
+ for(i=0; i<T; i++)
+ printf("%d", (int)(number/pow(2,i))%2);
+ printf("\n");
+
+}
+
+void supplement(long stabStateIndices[], int T)
+{
+
+ int K = round(log(T)/log(2)); // T = 2^K
+
+ //printf("K=%d\n", K);
+
+ int treeDepth, levelMaskDepth, maskIndex;
+
+ unsigned long levelMask;
+
+ unsigned long long maskList[K];
+
+ unsigned long given = (unsigned long)stabStateIndices[0]; // the given bitstring
+
+ unsigned long bitstring;
+
+ int bitstringCounter = 1; // bitstringCounter starts at 1 since 0 is the given bitstring
+
+ int i;
+
+ // first added bitstring is alpha ... alpha for alpha = 01
+ bitstring = 1;
+ for(i=1; i<T/2; i++) {
+ bitstring *= 4;
+ bitstring += 1;
+ }
+ //printf("%lu\n", bitstring);
+ stabStateIndices[bitstringCounter] = given^(bitstring);
+ bitstringCounter++;
+
+ // second added bitstring is beta ... beta for beta = 10
+ bitstring = 2;
+ for(i=1; i<T/2; i++) {
+ bitstring *= 4;
+ bitstring += 2;
+ }
+ //printf("%lu\n", bitstring);
+ stabStateIndices[bitstringCounter] = given^bitstring;
+ bitstringCounter++;
+
+ maskList[0] = round(pow(2,pow(2,K)))-1;
+
+ for(treeDepth=1; treeDepth<=K-1; treeDepth++) {
+
+ levelMask = round(pow(2,pow(2,K-treeDepth))-1);
+ // repeat this unit-cell mask to the other bits on levelMask
+ for(levelMaskDepth=2; levelMaskDepth<=round(pow(2,treeDepth-1)); levelMaskDepth++) {
+ levelMask += round(pow(2,pow(2,K-treeDepth+1)))*levelMask;
+ }
+
+ //printf("levelMask=");
+ //printBinary(levelMask, T);
+
+ bitstring = stabStateIndices[1] ^ levelMask;
+
+ for(maskIndex=0; maskIndex < (int)(pow(2,treeDepth)); maskIndex++) {
+ stabStateIndices[bitstringCounter] = bitstring;
+ for(i=0; i<treeDepth; i++){
+ // if ith binary digit of maskIndex is 1 then XOR with ith string in maskList
+ stabStateIndices[bitstringCounter] ^= ((int)(maskIndex/(pow(2,i)))%2)*maskList[i];
+ }
+ //printf("new bitstring=");
+ //printBinary(stabStateIndices[bitstringCounter], T);
+ bitstringCounter++;
+ }
+
+ maskList[treeDepth] = levelMask;
+
+ }
+
+}
+
+
+int main(int argc, char *argv[])
+{
+
+ int T = (int)pow(2,4);
+ long* stabStateIndices;
+
+ stabStateIndices = calloc(T+1,sizeof(long));
+
+ stabStateIndices[0] = (int)(pow(2,T))-1;
+
+ supplement(&stabStateIndices[0], T);
+
+ int i;
+ for(i=0; i<T+1; i++)
+ printBinary(stabStateIndices[i], T);
+ printf("\n");
+
+ return 0;
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+#include <time.h>
+#include "matrix.h"
+#include "exponentialsum.h"
+#include "shrink.h"
+#include "extend.h"
+#include "measurepauli.h"
+#include "innerproduct.h"
+#include "sparsify.h"
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+int main(int argc, char *argv[])
+{
+
+ if(argc != 4) {
+ printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+ exit(0);
+ }
+
+ double additiveErrorDelta = atof(argv[1]); // additive error delta
+ double phi = PI*atof(argv[2]);//PI/4.0; // PI/4.0 is the T gate magic state
+ int coherentSampling = atoi(argv[3]); // perform coherent sampling (true=1 or false=0)
+
+ int N; // number of qubits
+ scanf("%d", &N);
+
+ int T; // number of T gate magic states (set to the first 'K' of the 'N' qubits -- the rest are set to the '0' computational basis state)
+ scanf("%d", &T);
+
+ printf("phi = %lf\n", phi);
+
+ int omega;
+ int alpha[N], beta[N], gamma[N], delta[N];
+ int Paulicounter = 0;
+
+ int i, j, k, l, m;
+
+ //|T> = e^(pi*i/8)/2*sqrt(4-2*sqrt(2))* (sqrt(-i)*(|0>+i|1>)/sqrt(2) + (|0>+|1>)/sqrt(2))
+
+ //double additiveErrorDelta = 0.1;
+
+ double complex coeffa = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*cexp(-PI*I*0.25)*I/sqrt(2.0)*(-I+cexp(-0.25*PI*I))*(-I+cexp(I*phi)))); // factor of cexp(PI*I/8.0)*cexp(-PI*I*0.25) comes from converting (|0>+|1>)/sqrt(2) under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ double complex coeffb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*phi)))); // factor of cexp(PI*I/8.0) comes from converting |0> under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ // alternative coefficient to use instead of coeffb to get overall entangled state
+ double complex coeffbb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*0.25*PI))));
+
+ int n1 = 1; int k1 = 1; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
+ int n2 = 1; int k2 = 1; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {0}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
+
+ long* stabStateIndices;
+ int numStabStates;
+
+ srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+
+ if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+ return 1;
+
+ //printf("checking: numStabStateIndices:\n");
+ //for(i=0; i<numStabStates; i++)
+ // printf("%d ", stabStateIndices[i]);
+ //printf("\n");
+ //fflush(stdout);
+
+ int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
+ double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
+ G = calloc(numStabStates,sizeof(int*)); GBar = calloc(numStabStates,sizeof(int*));
+ h = calloc(numStabStates,sizeof(int*));
+
+ J = calloc(numStabStates,sizeof(int*)); D = calloc(numStabStates,sizeof(int*)); Q = calloc(numStabStates,sizeof(int));
+
+ K = calloc(numStabStates, sizeof(int));
+
+ double complex origGamma[(int)numStabStates];
+ int *origK, *origQ, **origD, ***origJ;
+ int ***origG, ***origGBar, **origh;
+
+ origG = calloc(numStabStates,sizeof(int*)); origGBar = calloc(numStabStates,sizeof(int*));
+ origh = calloc(numStabStates,sizeof(int*));
+
+ origJ = calloc(numStabStates,sizeof(int*)); origD = calloc(numStabStates,sizeof(int*)); origQ = calloc(numStabStates,sizeof(int));
+
+ origK = calloc(numStabStates, sizeof(int));
+
+ int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+
+ double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
+
+ for(j=0; j<numStabStates; j++) {
+
+ combination = stabStateIndices[j];
+
+ K[j] = 0.0;
+
+ for(k=0; k<N; k++) {
+ K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ combination /= 2;
+ }
+ combination = stabStateIndices[j];
+ origK[j] = K[j];
+
+ Gamma[j] = 1.0;
+ Gamma[j] *= L1Norm/((double)numStabStates);
+
+ // the coefficients which are a product of 'coeffa', 'coeffb', 'coeffbb' (that are subsequently multiplied into Gamma[j]) is multiplied by 'norm'
+ //Gamma[j] *= norm;
+
+ G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
+ h[j] = calloc(N, sizeof(int));
+
+ if(K[j] > 0) {
+ J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
+ for(k=0; k<K[j]; k++)
+ J[j][k] = calloc(K[j], sizeof(int));
+ }
+
+ origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*));
+ origh[j] = calloc(N, sizeof(int));
+
+ if(K[j] > 0) {
+ origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int));
+ for(k=0; k<K[j]; k++)
+ origJ[j][k] = calloc(K[j], sizeof(int));
+ }
+
+ for(k=0; k<N; k++) {
+ G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
+ origG[j][k] = calloc(N, sizeof(int)); origGBar[j][k] = calloc(N, sizeof(int));
+ }
+
+ int Kcounter = 0; // Kcounter keeps track of the K<=N that we have added already to the G rows etc. for each combination that is indexed by the digits (base 3) of 'j' in that we go through with 'k'
+ int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
+
+ // if combination contains at least one instance of the second state, i.e. contains the 0 digit in binary, then we want to have it have one instance of coeffb instead of coeffbb
+ for(k=0; k<N; k++) {
+ if(combination%2==1) {
+ Gamma[j] *= coeffb/coeffbb;
+ break; // break out of loop
+ }
+ combination /= 2; // shift to the right by one (in base-2 arithmetic)
+ }
+ combination = stabStateIndices[j];
+
+ for(k=0; k<N; k++) {
+
+ Q[j] += ((combination%2)==1)*Q2 + ((combination%2)==0)*Q1;
+
+ Gamma[j] *= (((combination%2)==1)*coeffbb + ((combination%2)==0)*coeffa); // only assign coeffbb instead of coeffb; coeffb replaces one instance of coeffbb before this loop
+
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ for(l=0; l<Kcombo; l++) {
+ // D1 may have a different number of rows 'l' than D2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ D[j][Kcounter+l] = D1[l];
+ break;
+ case 1:
+ D[j][Kcounter+l] = D2[l];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ for(m=0; m<Kcombo; m++) {
+ // J1 may have a different number of rows 'l' than J2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ J[j][Kcounter+l][Kcounter+m] = J1[l][m];
+ break;
+ case 1:
+ J[j][Kcounter+l][Kcounter+m] = J2[l][m];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ }
+ }
+
+ for(l=0; l<n1; l++) { // assuming n1=n2
+ h[j][k*n1+l] = ((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l];
+ }
+ // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
+ for(l=0; l<Kcombo; l++) {
+ for(m=0; m<n1; m++) { // assuming n1=n2
+ G[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
+ GBar[j][Kcounter+l][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
+ }
+ }
+ Kcounter = Kcounter + Kcombo;
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+ //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
+
+ //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
+
+ combination /= 2; // shift to the right by one (in base-2 arithmetic)
+ }
+ //printf("!\n");
+
+ // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
+ combination = stabStateIndices[j];
+ for(k=0; k<(N); k++) {
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ //printf("Kcounter=%d\n", Kcounter);
+ // G and GBar rows that are outside the first 'k' spanning basis states
+ for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
+ //printf("l=%d\n", l);
+ for(m=0; m<n1; m++) { // assuming n1=n2=n3
+ /* printf("m=%d\n", m); */
+ /* printf("Kcounter+l=%d\n", Kcounter+l); */
+ /* printf("k*n1+m=%d\n", k*n1+m); */
+ G[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m];
+ GBar[j][Kcounter+l-Kcombo][k*n1+m] = ((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m];
+ }
+ }
+ Kcounter = Kcounter + (n1-Kcombo);
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+
+ combination /= 2;
+ }
+ for(k=0; k<N; k++) {
+ memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
+ }
+ for(k=0; k<K[j]; k++) {
+ memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
+ }
+
+ memcpy(origh[j], h[j], N*sizeof(int));
+ memcpy(origD[j], D[j], K[j]*sizeof(int));
+
+ }
+ //exit(0);
+ memcpy(origGamma, Gamma, numStabStates*sizeof(double complex));
+
+ memcpy(origQ, Q, numStabStates*sizeof(int));
+
+ while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) {
+
+ Paulicounter++;
+ if(Paulicounter > N) {
+ printf("Error: Number of Paulis is greater than N!\n");
+ return 1;
+ }
+
+ // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
+ // Y_i = -I*Z*X
+ for(i=0; i<N; i++) {
+ if(delta[i]){
+ omega += 3; // -I = I^3
+ beta[i] = delta[i];
+ gamma[i] = delta[i];
+ }
+ }
+
+
+ for(j=0; j<numStabStates; j++) { // the kets
+
+ Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+
+ }
+
+ }
+
+ double complex amplitude = 0.0 + 0.0*I;
+ for(i=0; i<numStabStates; i++) { // the bras
+ for(j=0; j<numStabStates; j++) {
+ double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+ amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+ }
+ }
+
+ //printf("numStabStates=%d\n", numStabStates);
+ printf("L1Norm=%lf\n", L1Norm);
+
+ printf("\namplitude:\n");
+ if(creal(amplitude+0.00000001)>0)
+ printf("%lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+ else
+ printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude)));
+ //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude)));
+ printf("\nabs(amplitude):\n");
+ printf("%lf\n", cabs(amplitude));
+
+ return 0;
+
+}
+
+
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int newomega, newalpha, newbeta, newgamma, newdelta;
+ int i;
+
+ if(scanf("%d", &newomega) != EOF) {
+ *omega = newomega;
+ for(i=0; i<numqubits; i++) {
+ if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
+ printf("Error: Too few input coeffs!\n");
+ exit(0);
+ }
+ if(newalpha+newbeta+newgamma+newdelta > 1) {
+ printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+ exit(0);
+ }
+ alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
+ }
+ return 1;
+ } else
+ return 0;
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+#include <time.h>
+#include "matrix.h"
+#include "exponentialsum.h"
+#include "shrinkstar.h"
+#include "extend.h"
+#include "measurepauli.h"
+#include "innerproduct.h"
+#include "randomstabilizerstate.h"
+#include "sparsify.h"
+
+#define ZEROTHRESHOLD (0.00000001)
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+int main(int argc, char *argv[])
+{
+
+ if(argc != 5) {
+ printf("weaksim_rellerr argument: \"number of stabilizer state samples\" \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=yes)\"\n");
+ exit(0);
+ }
+
+ int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples
+ double additiveErrorDelta = atof(argv[2]); // additive error delta
+ double phi = PI*atof(argv[3]);//PI/4.0; // PI/4.0 is the T gate magic state
+ int coherentSampling = atoi(argv[4]); // perform coherent sampling (true=1 or false=0)
+
+ int N; // number of qubits
+ scanf("%d", &N);
+
+ int T; // number of T gate magic states (set to the first 'K' of the 'N' qubits -- the rest are set to the '0' computational basis state)
+ scanf("%d", &T);
+
+ printf("phi = %lf\n", phi);
+
+ int omega[N]; // max of N measurements
+ int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
+ int Paulicounter = 0;
+
+ int i, j, k, l, m;
+
+ FILE *fp;
+ char buff[255];
+
+ srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
+
+ fp = fopen("Pd.txt", "r");
+
+ if(fscanf(fp, "%s", buff) == EOF) {
+ printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
+ return 1;
+ }
+
+ double** Pd;
+
+ int PdN = atoi(buff);
+
+ Pd = calloc(PdN, sizeof(double*));
+ for(i=0; i<PdN; i++)
+ Pd[i] = calloc(PdN+1, sizeof(double));
+
+ double tmp;
+
+ for(i=1; i<PdN; i++) {
+ tmp = 0.0;
+ for(j=0; j<=i; j++) {
+ if(fscanf(fp, "%s", buff) == EOF) {
+ printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
+ return 1;
+ }
+ Pd[i][j] = atof(buff);
+ //printf("%e ", Pd[i][j]);
+ tmp += Pd[i][j];
+ }
+ //printf("\n");
+ //printf("total=%f\n", tmp);
+ }
+
+ double complex amplitude;
+
+ double complex coeffa = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*cexp(-PI*I*0.25)*I/sqrt(2.0)*(-I+cexp(-0.25*PI*I))*(-I+cexp(I*phi)))); // factor of cexp(PI*I/8.0)*cexp(-PI*I*0.25) comes from converting (|0>+|1>)/sqrt(2) under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ double complex coeffb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*phi)))); // factor of cexp(PI*I/8.0) comes from converting |0> under e^(pi*i/8) H S^\dagger to take it from |H> to |T>
+ // alternative coefficient to use instead of coeffb to get overall entangled state
+ double complex coeffbb = cexp(I*carg(cexp(PI*I/8.0)*0.5*csqrt(4.0-2.0*csqrt(2.0))*I/sqrt(2.0)*(1.0+cexp(-0.25*PI*I))*(1.0-cexp(I*0.25*PI))));
+
+ int n1 = 1; int k1 = 1; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} };
+ int n2 = 1; int k2 = 1; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {0}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} };
+
+ long* stabStateIndices;
+ int numStabStates;
+
+ srand((unsigned)time(NULL)); // seeding the random number generator for sparsify()
+
+
+ if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling))
+ return 1;
+
+ //printf("checking: numStabStateIndices:\n");
+ //for(i=0; i<numStabStates; i++)
+ // printf("%ld ", stabStateIndices[i]);
+ //printf("\n");
+ //fflush(stdout);
+
+ int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
+ double complex Gamma[(int)numStabStates]; // prefactor in front of resultant state
+ G = calloc(numStabStates,sizeof(int*)); GBar = calloc(numStabStates,sizeof(int*));
+ h = calloc(numStabStates,sizeof(int*));
+
+ J = calloc(numStabStates,sizeof(int*)); D = calloc(numStabStates,sizeof(int*)); Q = calloc(numStabStates,sizeof(int));
+
+ K = calloc(numStabStates, sizeof(int));
+
+ int origK, origQ, *origD;
+ int **origJ;
+ int **origG, **origGBar;
+ int *origh;
+ double complex origGamma;
+
+ long combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
+
+ double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
+
+ for(j=0; j<numStabStates; j++) {
+
+ combination = stabStateIndices[j];
+
+ K[j] = 0.0;
+
+ for(k=0; k<N; k++) {
+ K[j] += (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ combination /= 2;
+ }
+ combination = j;
+
+ Gamma[j] = 1.0;
+ Gamma[j] *= L1Norm/((double)numStabStates);
+
+ // the coefficients which are a product of 'coeffa', 'coeffb', 'coeffbb' (that are subsequently multiplied into Gamma[j]) is multiplied by 'norm'
+ //Gamma[j] *= norm;
+
+ G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
+ h[j] = calloc(N, sizeof(int));
+
+ if(K[j] > 0) {
+ J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
+ for(k=0; k<K[j]; k++)
+ J[j][k] = calloc(K[j], sizeof(int));
+ }
+
+ for(k=0; k<N; k++) {
+ G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
+ }
+
+ int Kcounter = 0; // Kcounter keeps track of the K<=N that we have added already to the G rows etc. for each combination that is indexed by the digits (base 3) of 'j' in that we go through with 'k'
+ int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
+
+ // if combination contains at least one instance of the second state, i.e. contains the 0 digit in binary, then we want to have it have one instance of coeffb instead of coeffbb
+ for(k=0; k<N; k++) {
+ if(combination%2==1) {
+ Gamma[j] *= coeffb/coeffbb;
+ break; // break out of loop
+ }
+ combination /= 2; // shift to the right by one (in base-2 arithmetic)
+ }
+ combination = stabStateIndices[j];
+
+ for(k=0; k<N; k++) {
+
+ Q[j] += (((combination%2)==1)*Q2 + ((combination%2)==0)*Q1);
+
+
+ Gamma[j] *= (((combination%2)==1)*coeffbb + ((combination%2)==0)*coeffa); // only assign coeffbb instead of coeffb; coeffb replaces one instance of coeffbb before this loop
+
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ for(l=0; l<Kcombo; l++) {
+ // D1 has a different number of rows 'l' than D2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ D[j][Kcounter+l] = D1[l];
+ break;
+ case 1:
+ D[j][Kcounter+l] = D2[l];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ for(m=0; m<Kcombo; m++) {
+ // J1 has a different number of rows 'l' than J2 so you need to use something like 'switch' to check combination%2 without going out of bound of J1
+ switch(combination%2) {
+ case 0:
+ J[j][Kcounter+l][Kcounter+m] = J1[l][m];
+ break;
+ case 1:
+ J[j][Kcounter+l][Kcounter+m] = J2[l][m];
+ break;
+ default:
+ printf("error");
+ return 1;
+ }
+ }
+ }
+
+ for(l=0; l<n1; l++) { // assuming n1=n2
+ h[j][k*n1+l] = (((combination%2)==1)*h2[l] + ((combination%2)==0)*h1[l]);
+ }
+ // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
+ for(l=0; l<Kcombo; l++) {
+ for(m=0; m<n1; m++) { // assuming n1=n2
+ G[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
+ GBar[j][Kcounter+l][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+ }
+ }
+ Kcounter = Kcounter + Kcombo;
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+ //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
+
+ //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
+
+ combination /= 2; // shift to the right by one (in base-7 arithmetic)
+ }
+ //printf("!\n");
+
+ // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
+ combination = j;
+ for(k=0; k<(N); k++) {
+ Kcombo = (((combination%2)==1)*k2 + ((combination%2)==0)*k1);
+ //printf("Kcounter=%d\n", Kcounter);
+ // G and GBar rows that are outside the first 'k' spanning basis states
+ for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
+ //printf("l=%d\n", l);
+ for(m=0; m<n1; m++) { // assuming n1=n2=n3
+ /* printf("m=%d\n", m); */
+ /* printf("Kcounter+l=%d\n", Kcounter+l); */
+ /* printf("k*n1+m=%d\n", k*n1+m); */
+ G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*G2[l][m] + ((combination%2)==0)*G1[l][m]);
+ GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%2)==1)*GBar2[l][m] + ((combination%2)==0)*GBar1[l][m]);
+ }
+ }
+ Kcounter = Kcounter + (n1-Kcombo);
+
+ /* printf("intermediate G[%d]:\n", j); */
+ /* printMatrix(G[j], N, N); */
+ /* printf("intermediate GBar[%d]:\n", j); */
+ /* printMatrix(GBar[j], N, N); */
+
+ combination /= 2;
+ }
+
+ /*printf("G[%d]:\n", j);
+ printMatrix(G[j], N, N);
+ printf("GBar[%d]:\n", j);
+ printMatrix(GBar[j], N, N);
+
+ printf("h[%d]:\n", j);
+ printVector(h[j], N);
+
+ printf("J[%d]:\n", j);
+ printMatrix(J[j], K[j], K[j]);
+
+ printf("D[%d]:\n", j);
+ printVector(D[j], K[j]);
+
+ printf("Q[%d]=%d\n", j, Q[j]);*/
+
+ }
+ //exit(0);
+
+ while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
+
+ if((Paulicounter+1) > N) {
+ printf("Error: Number of Paulis is greater than N!\n");
+ return 1;
+ }
+
+ // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
+ // Y_i = -I*Z*X
+ for(i=0; i<N; i++) {
+ if(delta[Paulicounter][i]){
+ omega[Paulicounter] += 3; // -I = I^3
+ beta[Paulicounter][i] = delta[Paulicounter][i];
+ gamma[Paulicounter][i] = delta[Paulicounter][i];
+ }
+ }
+
+ /*printf("*******\n");
+ printf("*******\n");
+ printf("omega=%d\n", omega);
+ printf("X:\n");
+ printVector(gamma, N);
+ printf("Z:\n");
+ printVector(beta, N);
+ printf("*******\n");
+ printf("*******\n");*/
+
+ //for(j=0; j<numStabStates; j++) { // the kets
+
+ /*printf("========\n");
+ printf("before:\n");
+ printf("K=%d\n", K[j]);
+ printf("h:\n");
+ printVector(h[j], N);
+ printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
+ printf("G:\n");
+ printMatrix(G[j], N, N);
+ printf("GBar:\n");
+ printMatrix(GBar[j], N, N);
+ printf("Q=%d\n", Q[j]);
+ printf("D:\n");
+ printVector(D[j], K[j]);
+ printf("J:\n");
+ printMatrix(J[j], K[j], K[j]);*/
+ //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+ /*printf("\nafter:\n");
+ printf("K=%d\n", K[j]);
+ printf("h:\n");
+ printVector(h[j], N);
+ printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
+ printf("G:\n");
+ printMatrix(G[j], N, N);
+ printf("GBar:\n");
+ printMatrix(GBar[j], N, N);
+ printf("Q=%d\n", Q[j]);
+ printf("D:\n");
+ printVector(D[j], K[j]);
+ printf("J:\n");
+ printMatrix(J[j], K[j], K[j]);*/
+
+ //}
+
+ Paulicounter++;
+ }
+
+ amplitude = 0.0 + 0.0*I;
+ for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
+ //printf("i=%d\n", i);
+
+ randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
+
+ origGamma = 1.0 + 0.0*I;
+
+ for(k=0; k<Paulicounter; k++) {
+ origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
+ //printf("k=%d\n", k);
+ }
+ /*printf("origK=%d\n", origK);
+ printf("origG:\n");
+ printMatrix(origG, N, N);
+ printf("origGBar:\n");
+ printMatrix(origGBar, N, N);
+ printf("origh:\n");
+ printVector(origh, N);*/
+
+ double complex stabstateaverage = 0.0 + 0.0*I;
+
+ for(j=0; j<numStabStates; j++) {
+ //printf("j=%d\n", j);
+ double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK, origh, origG, origGBar, origQ, origD, origJ);
+ stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
+ }
+ amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
+
+ deallocate_mem(&origG, N);
+ deallocate_mem(&origGBar, N);
+ free(origh);
+ deallocate_mem(&origJ, origK);
+ free(origD);
+ }
+
+ //printf("numStabStates=%d\n", numStabStates);
+ printf("L1Norm=%lf\n", L1Norm);
+
+ printf("\namplitude:\n");
+ if(creal(amplitude+ZEROTHRESHOLD)>0)
+ printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+ else
+ printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
+ printf("\nabs(amplitude):\n");
+ printf("%lf\n", cabs(amplitude));
+
+
+ for(i=0; i<PdN; i++)
+ free(Pd[i]);
+ free(Pd);
+
+ return 0;
+}
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int newomega, newalpha, newbeta, newgamma, newdelta;
+ int i;
+
+ if(scanf("%d", &newomega) != EOF) {
+ *omega = newomega;
+ for(i=0; i<numqubits; i++) {
+ if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
+ printf("Error: Too few input coeffs!\n");
+ exit(0);
+ }
+ if(newalpha+newbeta+newgamma+newdelta > 1) {
+ printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+ exit(0);
+ }
+ alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
+ }
+ return 1;
+ } else
+ return 0;
+
+}
projects / weak_simulation_stab_extent.git / commitdiff