From d914e7c183347b1340988e0869a6e053dd573991 Mon Sep 17 00:00:00 2001 From: Lucas K Date: Wed, 11 Aug 2021 23:26:51 -0700 Subject: [PATCH] initial commit with license --- LICENSE.txt | 11 ++ Makefile | 83 ++++++++ exponentialsum.c | 416 +++++++++++++++++++++++++++++++++++++++ extend.c | 177 +++++++++++++++++ innerproduct.c | 190 ++++++++++++++++++ matrix.c | 243 +++++++++++++++++++++++ measurepauli.c | 165 ++++++++++++++++ randomstabilizerstate.c | 131 +++++++++++++ shrink.c | 380 ++++++++++++++++++++++++++++++++++++ sparsify.c | 143 ++++++++++++++ supplement.c | 109 +++++++++++ weaksim.c | 327 +++++++++++++++++++++++++++++++ weaksim_relerr.c | 423 ++++++++++++++++++++++++++++++++++++++++ 13 files changed, 2798 insertions(+) create mode 100644 LICENSE.txt create mode 100644 Makefile create mode 100644 exponentialsum.c create mode 100644 extend.c create mode 100644 innerproduct.c create mode 100644 matrix.c create mode 100644 measurepauli.c create mode 100644 randomstabilizerstate.c create mode 100644 shrink.c create mode 100644 sparsify.c create mode 100644 supplement.c create mode 100644 weaksim.c create mode 100644 weaksim_relerr.c diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..a5165de --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,11 @@ +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. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..8c938ba --- /dev/null +++ b/Makefile @@ -0,0 +1,83 @@ +#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 + + diff --git a/exponentialsum.c b/exponentialsum.c new file mode 100644 index 0000000..cbf58d9 --- /dev/null +++ b/exponentialsum.c @@ -0,0 +1,416 @@ +#include +#include +#include +#include +#include + +#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->datadata][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->datadata][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->datadata][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; adata])); + setWalker = setWalker->next; + } + setWalker = setDa; + setWalker2 = setDb; + for(a=0; adata])) + 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; adata])); + setWalker = setWalker->next; + } + setWalker = setDa; + setWalker2 = setDb; + for(a=0; adata])) + 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; adata] + 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; adata][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"); +} + diff --git a/extend.c b/extend.c new file mode 100644 index 0000000..18b2227 --- /dev/null +++ b/extend.c @@ -0,0 +1,177 @@ +#include +#include +#include +#include + +#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; adata][a] = (GBar[setWalker->data][a] + GBar[i][a])%2; + // update G rows, g + for(a=0; adata][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"); +} diff --git a/innerproduct.c b/innerproduct.c new file mode 100644 index 0000000..5d82cdb --- /dev/null +++ b/innerproduct.c @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include + +#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 +#include +#include + +#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; i0) + printf("\n"); + +} + +void printMatrix(int** a, int rows, int cols) +{ + int i, j; + printf("Matrix[%d][%d]\n", rows, cols); + for(i=0; i +#include +#include +#include + +#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 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); + +} diff --git a/randomstabilizerstate.c b/randomstabilizerstate.c new file mode 100644 index 0000000..605affd --- /dev/null +++ b/randomstabilizerstate.c @@ -0,0 +1,131 @@ +#include +#include +#include +#include +#include + +#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 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= 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 +#include +#include +#include +#include + +#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; adata][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(idata != i) { + R[setWalker->data][i] = 0; // reset R + RT[i][setWalker->data] = 0; + } + + // update GBar rows, gbar + for(a=0; adata][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(ic) */ + /* tmp += ((*J)[c][*k-1])*R[c][c]*R[c][*k-1]; */ + /* else if((*k-1)i) */ + /* tmp += ((*J)[i][*k-1])*R[c][i]*R[c][*k-1]; */ + /* else if((*k-1)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"); +} diff --git a/sparsify.c b/sparsify.c new file mode 100644 index 0000000..8e164fc --- /dev/null +++ b/sparsify.c @@ -0,0 +1,143 @@ +#include +#include +#include +#include + +#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 +#include +#include + + +void printBinary(unsigned long number, int T) { + + int i; + + for(i=0; i +#include +#include +#include +#include +#include +#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 0) { + J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int)); + for(k=0; k 0) { + origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int)); + for(k=0; k 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; i0) + 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 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; + +} diff --git a/weaksim_relerr.c b/weaksim_relerr.c new file mode 100644 index 0000000..2cc45fb --- /dev/null +++ b/weaksim_relerr.c @@ -0,0 +1,423 @@ +#include +#include +#include +#include +#include +#include +#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+|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 0) { + J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int)); + for(k=0; k 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; i0) + 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 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; + +} -- 2.30.2