From d914e7c183347b1340988e0869a6e053dd573991 Mon Sep 17 00:00:00 2001 From: Lucas K Date: Wed, 11 Aug 2021 23:26:51 -0700 Subject: [PATCH 1/3] 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 From 7bc25656311c391af6a4d5814cc82080981c63c4 Mon Sep 17 00:00:00 2001 From: Lucas Date: Tue, 1 Feb 2022 16:35:48 -0800 Subject: [PATCH 2/3] update from v0.1 to v1.0 --- LICENSE.txt | 2 +- Makefile | 26 +- Pd.txt | 51 ++++ README.txt | 59 +++++ exponentialsum.c | 2 +- exponentialsum.h | 3 + extend.h | 1 + innerproduct.h | 1 + innerproduct_equatorial.c | 190 +++++++++++++++ innerproduct_equatorial.h | 1 + matrix.h | 18 ++ measurepauli.h | 2 + multipauli.c | 301 +++++++++++++++++++++++ randomstabilizerstate.c | 7 +- randomstabilizerstate.h | 1 + randomstabilizerstate_equatorial.c | 66 +++++ randomstabilizerstate_equatorial.h | 1 + shrink.h | 1 + shrinkstar.c | 380 +++++++++++++++++++++++++++++ shrinkstar.h | 1 + sparsify.c | 78 +++++- sparsify.h | 4 + supplement.c | 10 +- supplement.h | 2 + supplement2.c | 114 +++++++++ supplement2.h | 2 + timing_t_enforceddelta.bash | 127 ++++++++++ 27 files changed, 1425 insertions(+), 26 deletions(-) create mode 100644 Pd.txt create mode 100644 README.txt create mode 100644 exponentialsum.h create mode 100644 extend.h create mode 100644 innerproduct.h create mode 100644 innerproduct_equatorial.c create mode 100644 innerproduct_equatorial.h create mode 100644 matrix.h create mode 100644 measurepauli.h create mode 100644 multipauli.c create mode 100644 randomstabilizerstate.h create mode 100644 randomstabilizerstate_equatorial.c create mode 100644 randomstabilizerstate_equatorial.h create mode 100644 shrink.h create mode 100644 shrinkstar.c create mode 100644 shrinkstar.h create mode 100644 sparsify.h create mode 100644 supplement.h create mode 100644 supplement2.c create mode 100644 supplement2.h create mode 100644 timing_t_enforceddelta.bash diff --git a/LICENSE.txt b/LICENSE.txt index a5165de..d275d40 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,6 +1,6 @@ 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. +Copyright 2022 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: diff --git a/Makefile b/Makefile index 8c938ba..839e215 100644 --- a/Makefile +++ b/Makefile @@ -1,13 +1,13 @@ #IDIR =../include CC=gcc -std=c99 CFLAGS=-Wall -LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.so libinnerproduct.so +LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.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_relerr: weaksim_relerr.c matrix exponentialsum shrink shrinkstar extend measurepauli innerproduct randomstabilizerstate supplement supplement2 sparsify + $(CC) -o $@ weaksim_relerr.c $(CFLAGS) $(LIBS) libshrink.so libshrinkstar.so librandomstabilizerstate.so libsparsify.so libsupplement.so libsupplement2.so libinnerproduct.so -weaksim: weaksim.c matrix exponentialsum shrink extend measurepauli innerproduct sparsify - $(CC) -o $@ weaksim.c $(CFLAGS) $(LIBS) libshrink.so libsparsify.so +weaksim: weaksim.c matrix exponentialsum shrink extend measurepauli innerproduct supplement supplement2 sparsify + $(CC) -o $@ weaksim.c $(CFLAGS) $(LIBS) libshrink.so libsparsify.so libsupplement.so libsupplement2.so libinnerproduct.so -fopenmp module_sparsify_test: module_sparsify_test matrix sparsify $(CC) -o $@ module_sparsify_test.c $(CFLAGS) libmatrix.so libsparsify.so @@ -32,6 +32,10 @@ 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_equatorial: innerproduct_equatorial.h innerproduct_equatorial.c + $(CC) -c -Wall -fpic innerproduct_equatorial.c + $(CC) -shared -o libinnerproduct_equatorial.so innerproduct_equatorial.o -lm libextend.so libshrink.so libexponentialsum.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 @@ -44,6 +48,10 @@ shrinkstar: shrinkstar.h shrinkstar.c $(CC) -c -Wall -fpic shrinkstar.c $(CC) -shared -o libshrinkstar.so shrinkstar.o -lm libmatrix.so +randomstabilizerstate_equatorial: randomstabilizerstate_equatorial.h randomstabilizerstate_equatorial.c + $(CC) -c -Wall -fpic randomstabilizerstate_equatorial.c + $(CC) -shared -o librandomstabilizerstate_equatorial.so randomstabilizerstate_equatorial.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 @@ -52,9 +60,13 @@ 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 +supplement2: supplement2.h supplement2.c + $(CC) -c -Wall -fpic supplement2.c + $(CC) -shared -o libsupplement2.so supplement2.o -lm + +sparsify: sparsify.h sparsify.c supplement2 $(CC) -c -Wall -fpic sparsify.c - $(CC) -shared -o libsparsify.so sparsify.o -lm libmatrix.so libsupplement.so + $(CC) -shared -o libsparsify.so sparsify.o -lm libmatrix.so libsupplement2.so randominputcommutingHermitianPauli: randominputcommutingHermitianPauli.c $(CC) -o randominputcommutingHermitianPauli randominputcommutingHermitianPauli.c diff --git a/Pd.txt b/Pd.txt new file mode 100644 index 0000000..6245c22 --- /dev/null +++ b/Pd.txt @@ -0,0 +1,51 @@ +50 +0.666667 0.333333 +0.533333 0.4 0.0666667 +0.474074 0.414815 0.103704 0.00740741 +0.446187 0.418301 0.122004 0.0130719 0.00043573 +0.432667 0.419146 0.130983 0.0163729 0.000818644 0.0000132039 +0.42601 0.419354 0.135416 0.0181361 0.00105794 0.0000255953 2.03137E-7 +0.422708 0.419405 0.137617 0.0190453 0.00119033 0.000033598 3.99976E-7 1.57471E-9 +0.421063 0.419418 0.138714 0.0195066 0.0012598 0.0000380989 5.29151E-7 3.12491E-9 6.12727E-12 +0.420242 0.419421 0.139261 0.019739 0.00129537 0.0000404804 6.02387E-7 4.1503E-9 1.22068E-11 1.1944E-14 +0.419832 0.419422 0.139534 0.0198556 0.00131337 0.0000417047 6.41292E-7 4.73395E-9 1.62439E-11 2.38414E-14 1.16527E-17 +0.419627 0.419422 0.139671 0.019914 0.00132242 0.0000423253 6.61333E-7 5.04461E-9 1.85464E-11 3.17575E-14 2.32826E-17 5.68701E-21 +0.419525 0.419422 0.139739 0.0199432 0.00132695 0.0000426377 6.71502E-7 5.2048E-9 1.97731E-11 3.62766E-14 3.10283E-17 1.13685E-20 1.38809E-24 +0.419474 0.419422 0.139773 0.0199579 0.00132923 0.0000427945 6.76624E-7 5.28613E-9 2.0406E-11 3.86856E-14 3.54523E-17 1.51543E-20 2.7755E-24 1.69424E-28 +0.419448 0.419422 0.13979 0.0199652 0.00133036 0.000042873 6.79195E-7 5.3271E-9 2.07274E-11 3.99286E-14 3.78112E-17 1.7317E-20 3.70022E-24 3.38807E-28 1.03402E-32 +0.419435 0.419422 0.139799 0.0199688 0.00133093 0.0000429123 6.80482E-7 5.34766E-9 2.08893E-11 4.05599E-14 3.90285E-17 1.84704E-20 4.22857E-24 4.51715E-28 2.06791E-32 3.15548E-37 +0.419429 0.419422 0.139803 0.0199707 0.00133122 0.0000429319 6.81127E-7 5.35797E-9 2.09706E-11 4.0878E-14 3.96468E-17 1.90656E-20 4.51033E-24 5.1623E-28 2.75713E-32 6.31077E-37 4.81481E-42 +0.419426 0.419422 0.139805 0.0199716 0.00133136 0.0000429418 6.81449E-7 5.36312E-9 2.10113E-11 4.10377E-14 3.99584E-17 1.93679E-20 4.65576E-24 5.50637E-28 3.15096E-32 8.41423E-37 9.62947E-42 3.67338E-47 +0.419424 0.419422 0.139806 0.019972 0.00133143 0.0000429467 6.8161E-7 5.3657E-9 2.10317E-11 4.11177E-14 4.01148E-17 1.95203E-20 4.72962E-24 5.68395E-28 3.361E-32 9.61619E-37 1.28392E-41 7.3467E-47 1.40128E-52 +0.419423 0.419422 0.139807 0.0199723 0.00133146 0.0000429491 6.81691E-7 5.36699E-9 2.10419E-11 4.11577E-14 4.01931E-17 1.95968E-20 4.76684E-24 5.77415E-28 3.4694E-32 1.02572E-36 1.46733E-41 9.79556E-47 2.80254E-52 2.67272E-58 +0.419423 0.419422 0.139807 0.0199724 0.00133148 0.0000429504 6.81731E-7 5.36763E-9 2.1047E-11 4.11778E-14 4.02323E-17 1.96351E-20 4.78553E-24 5.8196E-28 3.52447E-32 1.05881E-36 1.56515E-41 1.11949E-46 3.73672E-52 5.34543E-58 2.5489E-64 +0.419423 0.419422 0.139807 0.0199724 0.00133149 0.000042951 6.81751E-7 5.36796E-9 2.10495E-11 4.11878E-14 4.02519E-17 1.96543E-20 4.79489E-24 5.84242E-28 3.55222E-32 1.07561E-36 1.61564E-41 1.19412E-46 4.27053E-52 7.12723E-58 5.0978E-64 1.21541E-70 +0.419423 0.419422 0.139807 0.0199725 0.0013315 0.0000429513 6.81761E-7 5.36812E-9 2.10508E-11 4.11928E-14 4.02617E-17 1.96639E-20 4.79957E-24 5.85385E-28 3.56614E-32 1.08408E-36 1.64128E-41 1.23264E-46 4.55523E-52 8.1454E-58 6.79706E-64 2.43082E-70 2.89776E-77 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429515 6.81766E-7 5.3682E-9 2.10514E-11 4.11953E-14 4.02667E-17 1.96687E-20 4.80192E-24 5.85957E-28 3.57312E-32 1.08833E-36 1.6542E-41 1.25221E-46 4.70217E-52 8.68843E-58 7.76807E-64 3.24109E-70 5.79552E-77 3.4544E-84 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429515 6.81769E-7 5.36824E-9 2.10518E-11 4.11966E-14 4.02691E-17 1.96711E-20 4.80309E-24 5.86243E-28 3.57661E-32 1.09046E-36 1.66069E-41 1.26207E-46 4.77681E-52 8.9687E-58 8.28594E-64 3.7041E-70 7.72737E-77 6.9088E-84 2.05898E-91 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.8177E-7 5.36826E-9 2.10519E-11 4.11972E-14 4.02703E-17 1.96723E-20 4.80368E-24 5.86386E-28 3.57836E-32 1.09153E-36 1.66394E-41 1.26702E-46 4.81442E-52 9.11106E-58 8.55322E-64 3.95104E-70 8.83127E-77 9.21174E-84 4.11797E-91 6.13625E-99 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81771E-7 5.36827E-9 2.1052E-11 4.11975E-14 4.02709E-17 1.96729E-20 4.80397E-24 5.86458E-28 3.57924E-32 1.09206E-36 1.66557E-41 1.2695E-46 4.8333E-52 9.1828E-58 8.68899E-64 4.0785E-70 9.42003E-77 1.05277E-83 5.49062E-91 1.22725E-98 9.14373E-107 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81771E-7 5.36827E-9 2.1052E-11 4.11976E-14 4.02713E-17 1.96732E-20 4.80412E-24 5.86494E-28 3.57967E-32 1.09233E-36 1.66638E-41 1.27074E-46 4.84276E-52 9.21881E-58 8.75741E-64 4.14324E-70 9.7239E-77 1.12295E-83 6.275E-91 1.63633E-98 1.82875E-106 6.81261E-115 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81771E-7 5.36828E-9 2.10521E-11 4.11977E-14 4.02714E-17 1.96733E-20 4.80419E-24 5.86511E-28 3.57989E-32 1.09246E-36 1.66679E-41 1.27136E-46 4.84749E-52 9.23685E-58 8.79175E-64 4.17586E-70 9.87824E-77 1.15918E-83 6.69333E-91 1.8701E-98 2.43833E-106 1.36252E-114 2.53789E-123 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81771E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02715E-17 1.96734E-20 4.80423E-24 5.8652E-28 3.58E-32 1.09253E-36 1.66699E-41 1.27167E-46 4.84986E-52 9.24588E-58 8.80895E-64 4.19223E-70 9.95603E-77 1.17758E-83 6.90925E-91 1.99477E-98 2.78666E-106 1.8167E-114 5.07579E-123 4.7272E-132 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81771E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02715E-17 1.96734E-20 4.80424E-24 5.86525E-28 3.58006E-32 1.09256E-36 1.66709E-41 1.27182E-46 4.85105E-52 9.2504E-58 8.81757E-64 4.20044E-70 9.99507E-77 1.18685E-83 7.01892E-91 2.05912E-98 2.97244E-106 2.07622E-114 6.76772E-123 9.45439E-132 4.40254E-141 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02715E-17 1.96734E-20 4.80425E-24 5.86527E-28 3.58008E-32 1.09258E-36 1.66714E-41 1.2719E-46 4.85164E-52 9.25265E-58 8.82187E-64 4.20454E-70 1.00146E-76 1.19151E-83 7.07418E-91 2.0918E-98 3.06832E-106 2.21464E-114 7.73453E-123 1.26059E-131 8.80509E-141 2.05009E-150 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02715E-17 1.96734E-20 4.80426E-24 5.86528E-28 3.5801E-32 1.09259E-36 1.66717E-41 1.27194E-46 4.85193E-52 9.25378E-58 8.82403E-64 4.2066E-70 1.00244E-76 1.19384E-83 7.10193E-91 2.10827E-98 3.11703E-106 2.28608E-114 8.25017E-123 1.44067E-131 1.17401E-140 4.10019E-150 4.77325E-160 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96734E-20 4.80426E-24 5.86529E-28 3.5801E-32 1.09259E-36 1.66718E-41 1.27196E-46 4.85208E-52 9.25435E-58 8.8251E-64 4.20763E-70 1.00293E-76 1.195E-83 7.11582E-91 2.11654E-98 3.14157E-106 2.32236E-114 8.5163E-123 1.53671E-131 1.34173E-140 5.46692E-150 9.54649E-160 5.55679E-170 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.09259E-36 1.66719E-41 1.27197E-46 4.85216E-52 9.25463E-58 8.82564E-64 4.20814E-70 1.00318E-76 1.19559E-83 7.12278E-91 2.12068E-98 3.15389E-106 2.34065E-114 8.65148E-123 1.58629E-131 1.43118E-140 6.24791E-150 1.27287E-159 1.11136E-169 3.23448E-180 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27197E-46 4.85219E-52 9.25477E-58 8.82591E-64 4.2084E-70 1.0033E-76 1.19588E-83 7.12626E-91 2.12275E-98 3.16006E-106 2.34983E-114 8.7196E-123 1.61146E-131 1.47734E-140 6.66443E-150 1.4547E-159 1.48181E-169 6.46896E-180 9.41357E-191 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85221E-52 9.25484E-58 8.82605E-64 4.20853E-70 1.00336E-76 1.19603E-83 7.128E-91 2.12379E-98 3.16315E-106 2.35443E-114 8.7538E-123 1.62415E-131 1.50079E-140 6.87941E-150 1.55168E-159 1.6935E-169 8.62528E-180 1.88271E-190 1.36985E-201 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85222E-52 9.25488E-58 8.82611E-64 4.20859E-70 1.00339E-76 1.1961E-83 7.12887E-91 2.12431E-98 3.16469E-106 2.35673E-114 8.77093E-123 1.63052E-131 1.51261E-140 6.98861E-150 1.60174E-159 1.8064E-169 9.85746E-180 2.51029E-190 2.73971E-201 9.96701E-213 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.2549E-58 8.82615E-64 4.20862E-70 1.00341E-76 1.19614E-83 7.1293E-91 2.12457E-98 3.16547E-106 2.35788E-114 8.7795E-123 1.63371E-131 1.51854E-140 7.04364E-150 1.62716E-159 1.86467E-169 1.05146E-179 2.8689E-190 3.65295E-201 1.9934E-212 3.62598E-224 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.2549E-58 8.82617E-64 4.20864E-70 1.00341E-76 1.19615E-83 7.12952E-91 2.1247E-98 3.16585E-106 2.35846E-114 8.78379E-123 1.63531E-131 1.52151E-140 7.07126E-150 1.63998E-159 1.89427E-169 1.08538E-179 3.06016E-190 4.1748E-201 2.65787E-212 7.25195E-224 6.59561E-236 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82617E-64 4.20865E-70 1.00342E-76 1.19616E-83 7.12963E-91 2.12476E-98 3.16605E-106 2.35875E-114 8.78594E-123 1.63611E-131 1.523E-140 7.0851E-150 1.64641E-159 1.90918E-169 1.10261E-179 3.15887E-190 4.45312E-201 3.03756E-212 9.66927E-224 1.31912E-235 5.99867E-248 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12968E-91 2.1248E-98 3.16614E-106 2.35889E-114 8.78701E-123 1.63651E-131 1.52375E-140 7.09203E-150 1.64963E-159 1.91667E-169 1.11129E-179 3.20901E-190 4.59676E-201 3.24007E-212 1.10506E-223 1.75883E-235 1.19973E-247 2.72788E-260 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12971E-91 2.12481E-98 3.16619E-106 2.35896E-114 8.78755E-123 1.63671E-131 1.52412E-140 7.09549E-150 1.65124E-159 1.92042E-169 1.11565E-179 3.23428E-190 4.66973E-201 3.34459E-212 1.17873E-223 2.01009E-235 1.59965E-247 5.45576E-260 6.20248E-273 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12973E-91 2.12482E-98 3.16622E-106 2.359E-114 8.78781E-123 1.63681E-131 1.5243E-140 7.09722E-150 1.65205E-159 1.9223E-169 1.11783E-179 3.24696E-190 4.7065E-201 3.39767E-212 1.21675E-223 2.1441E-235 1.82817E-247 7.27435E-260 1.2405E-272 7.05141E-286 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12973E-91 2.12482E-98 3.16623E-106 2.35901E-114 8.78795E-123 1.63686E-131 1.5244E-140 7.09809E-150 1.65245E-159 1.92324E-169 1.11892E-179 3.25332E-190 4.72496E-201 3.42443E-212 1.23607E-223 2.21326E-235 1.95004E-247 8.31354E-260 1.654E-272 1.41028E-285 4.00826E-299 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16623E-106 2.35902E-114 8.78802E-123 1.63688E-131 1.52444E-140 7.09852E-150 1.65265E-159 1.92371E-169 1.11947E-179 3.2565E-190 4.7342E-201 3.43786E-212 1.2458E-223 2.24839E-235 2.01295E-247 8.86778E-260 1.89028E-272 1.88038E-285 8.01652E-299 1.13922E-312 +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16624E-106 2.35903E-114 8.78805E-123 1.6369E-131 1.52447E-140 7.09874E-150 1.65275E-159 1.92394E-169 1.11974E-179 3.25809E-190 4.73883E-201 3.44459E-212 1.25069E-223 2.2661E-235 2.0449E-247 9.15384E-260 2.0163E-272 2.149E-285 1.06887E-298 2.27843E-312 0. +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16624E-106 2.35903E-114 8.78807E-123 1.6369E-131 1.52448E-140 7.09885E-150 1.6528E-159 1.92406E-169 1.11988E-179 3.25889E-190 4.74114E-201 3.44795E-212 1.25313E-223 2.27498E-235 2.061E-247 9.29913E-260 2.08134E-272 2.29227E-285 1.22157E-298 3.03791E-312 0. 0. +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16624E-106 2.35903E-114 8.78807E-123 1.6369E-131 1.52448E-140 7.0989E-150 1.65283E-159 1.92412E-169 1.11995E-179 3.25928E-190 4.7423E-201 3.44964E-212 1.25436E-223 2.27944E-235 2.06909E-247 9.37236E-260 2.11438E-272 2.36621E-285 1.303E-298 3.4719E-312 0. 0. 0. +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16624E-106 2.35903E-114 8.78808E-123 1.63691E-131 1.52449E-140 7.09893E-150 1.65284E-159 1.92415E-169 1.11998E-179 3.25948E-190 4.74288E-201 3.45048E-212 1.25497E-223 2.28166E-235 2.07313E-247 9.40911E-260 2.13103E-272 2.40377E-285 1.34504E-298 3.70336E-312 0. 0. 0. 0. +0.419422 0.419422 0.139807 0.0199725 0.0013315 0.0000429516 6.81772E-7 5.36828E-9 2.10521E-11 4.11978E-14 4.02716E-17 1.96735E-20 4.80426E-24 5.86529E-28 3.58011E-32 1.0926E-36 1.66719E-41 1.27198E-46 4.85223E-52 9.25491E-58 8.82618E-64 4.20865E-70 1.00342E-76 1.19617E-83 7.12974E-91 2.12483E-98 3.16624E-106 2.35903E-114 8.78808E-123 1.63691E-131 1.52449E-140 7.09894E-150 1.65285E-159 1.92416E-169 1.12E-179 3.25958E-190 4.74317E-201 3.4509E-212 1.25528E-223 2.28278E-235 2.07516E-247 9.42752E-260 2.13938E-272 2.4227E-285 1.36638E-298 3.82282E-312 0. 0. 0. 0. 0. diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..4d4c2fe --- /dev/null +++ b/README.txt @@ -0,0 +1,59 @@ +Weak simulation of multi-Pauli measurements using approximate correlated stabilizer state decompositions +------------------------------------------------------------------------------------------------------- + +See the draft at [https://arxiv.org/abs/XXXXXX] for an introduction and description of this code. The uncorrelated version is described in [PRX QUANTUM 2, 010345 (2021)] by Seddon, Regula, Pashayan, Ouyang, Campbell (SRPOC). + +Note: the correlated sampling is currently only implemented for a maximum of 2t-1 supplemental correlated states! + +Note: The stabilizer state representation used in this does not use the CH-form, which ostensibly will speed up simulations here by a factor of approximately 50. We use the implementation described in the Supporting Information of arXiv:1601.07601. A 50-fold improvement is described in arXiv:1808.00128v2. + +PREREQUISITES: + +You will need the openMP library to run weaksim.c (in particular, the omp.h header). +For instance, in Debian, you can obtain this by running: +$ sudo apt install libomp-dev + +BUILD: + +1. We start with the code that scales quadratically with approximate stabilizer rank. +$ make weaksim + +2. This produces libraries (*.so) of core functions like SHRINK, INNERPRODUCT, EXTEND, etc. since you might want to use them in other code. This means we need to put them somewhere the linker can find them. Follow 3a) for a temporary and easy method that will only be valid for the current shell session and follow 3b) for a persistent method that requires root access. + +3a) Update LD_LIBRARY_PATH to see the library: +$ export LD_LIBRARY_PATH=$PWD:$LD_LIBRARY_PATH +or +3b) If you have root privileges put the *so library /usr/local/lib or whatever library directory in your path. Then, use ldconfig to write the path in the config file: +$ sudo cp ./*so /usr/local/lib +$ sudo echo "/usr/local/lib" >> /etc/ld.so.conf +$ sudo ldconfig + +4. To run weaksim you need to specify how many parallel processors for openMP to use +a) if you have no hyperthreading +export OMP_NUM_THREADS=8 +b) if you have hyperthreading +export OMP_PROC_BIND=true +export OMP_NUM_THREADS=8 + + +5. Build code for generating random commuting Pauli inputs (necessary for the strongsim[1-12]_relerr.c code since it requires a Hermitian observable to get away with using the square root number of stabilizer states). +$ make randominputcommutingHermitianPauli2 + +6. Build Hilbert vector space code to check our results. +$ make multipauli + +7. Run timing analysis for fixed t and delta +$ bash timing_t_enforceddelta.bash +Output will be in tmp_$[t]_t.txt for the correlated case and in tmp_$[t]_iid.txt for the i.i.d. case and provide the second-longest and longest runtime, and the largest error (closest to the target delta). + +Details of source code files: + +weaksim.c takes four arguments corresponding to the target additive error, the angle of the magic state, a "coherent sampling" flag to indicate whether sampling should be correlated and if so how many maximal supplemental states to use (0=no; 1=t; 2=2t-1; 3=xi^3t/2), a seed to use for the i.i.d. random stabilizer state sampling. There is also an optional argument for the "sample number prefactor" which only takes the corresponding fraction of supplemental states indicated by the "coherent sampling" flag, i.e. if you want 0.5 t supplemental states then this argument should be 0.5 and the "coherent sampling" should be 1. + +weaksim_relerr.c is a prototype function that also uses random stabilizer state sampling to calculate outcome probabilities from the Haar measure. This is still in a beta stage. + +randominputcommutingHermitianPauli.c uses a brute-force random sampling of the commuting Paulis whereas *2.c uses random Clifford operations on initially separable Paulis to generate the commuting Paulis. *2.c is the only practically viable method for more than 30 Paulis on a desktop computer. + +NOTE: current strongsim*c only support the same number of T gates as number of qubits! + +Pd.txt is a tabulated list of P(d) values for strongsim[1-12]_relerror.c. diff --git a/exponentialsum.c b/exponentialsum.c index cbf58d9..53360e1 100644 --- a/exponentialsum.c +++ b/exponentialsum.c @@ -331,7 +331,7 @@ complex double exponentialsum(int *k, int *Q, int *D, int **J) { void append(struct Node** head_ref, int new_data) -{ +{ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node)); struct Node *last = *head_ref; diff --git a/exponentialsum.h b/exponentialsum.h new file mode 100644 index 0000000..0be116f --- /dev/null +++ b/exponentialsum.h @@ -0,0 +1,3 @@ +#define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 + +complex double exponentialsum(int *k, int *Q, int *D, int **J); diff --git a/extend.h b/extend.h new file mode 100644 index 0000000..24221c4 --- /dev/null +++ b/extend.h @@ -0,0 +1 @@ +void extend(int n, int *k, int *h, int **G, int **GBar, int *xi); diff --git a/innerproduct.h b/innerproduct.h new file mode 100644 index 0000000..09dc5aa --- /dev/null +++ b/innerproduct.h @@ -0,0 +1 @@ +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); diff --git a/innerproduct_equatorial.c b/innerproduct_equatorial.c new file mode 100644 index 0000000..37f1e2c --- /dev/null +++ b/innerproduct_equatorial.c @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include + +#include "matrix.h" +#include "shrink.h" +#include "extend.h" +#include "exponentialsum.h" +#include "innerproduct_equatorial.h" + + +/**************************************************************************** + * Note: Arguments are copied and not modified! + ****************************************************************************/ + +double complex innerproduct_equatorial(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 + +void deallocate_mem(complex double ***arr, int rows); +void printMatrix(complex double** a, int rows, int cols); +complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2); +complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2); +complex double** transp(complex double **a, int rows, int cols); +complex double** conjtransp(complex double **a, int rows, int cols); +complex double trace(complex double **a, int rows, int cols); +complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols); +complex double** addMatrix(complex double **A, complex double **B, int rows, int cols); +int readPaulicoeffs(int* omega, int* alpha, int* beta, int* gamma, int* delta, int numqubits); + +// order of matrix elements is [row][column]!!! + +static complex double (*(I2[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, 1.0+0.0*I} }; +static complex double (*(X[])) = { (double complex[]) {0.0*I, 1.0+0.0*I}, (double complex[]) {1.0+0.0*I, 0.0*I} }; +static complex double (*(Y[])) = { (double complex[]) {0.0*I, 0.0-1.0*I}, (double complex[]) {0.0+1.0*I, 0.0*I} }; +static complex double (*(Z[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, -1.0+0.0*I} }; + +int main() +{ + + int i, j; + + int N; // number of qubits + scanf("%d", &N); + + int K; // 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", &K); + + complex double **IN; // N-qubit identity matrix + IN = I2; + for(i=1; i 0) { + psiN = psiT; + for(i=1; i N) { + printf("Error: Number of Paulis is greater than N!\n"); + return 1; + } + + + printf("%d\n", omega); + for(i=0; i0) + printf("%.10lf %c %.10lf I\n", cabs(creal(tr)), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr))); + else + printf("%.10lf %c %.10lf I\n", creal(tr), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr))); + //printf("%lf %c %lf I\n", creal(tr), cimag(tr)>0?'+':'-' , cabs(cimag(tr))); + //printf("%lf\n", cabs(creal(tr))); + /* cabs the creal part because it's always positive, but sometimes the 0.0 gets a minus sign which is annoying to see when comparing outputs */ + + deallocate_mem(&psiNfinal, pow(2,N)); + + + // deallocate_mem(&psiN, pow(2,N)); + + return 0; +} + + +complex double** addMatrix(complex double **A, complex double **B, int rows, int cols) +{ + int i, j; + + complex double** C; + + C = calloc(cols, sizeof(complex double*)); + for(i=0; i 1) { + printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i); + exit(0); + } + } + return 1; + } else + return 0; + +} + + + diff --git a/randomstabilizerstate.c b/randomstabilizerstate.c index 605affd..60f5cd8 100644 --- a/randomstabilizerstate.c +++ b/randomstabilizerstate.c @@ -41,8 +41,12 @@ void randomstabilizerstate(int n, int* k, int** h, int*** G, int*** GBar, int* Q //d = rand()%(n+1); //printf("d=%d total=%f randPd=%f\n", d, total, randPd); //printf("!!!d=%d\n", d); + //d = 0; - *k = n; // we will get it to be k=n-d by caling shrinkstar d times below + printf("!!\n"); + fflush(stdout); + + *k = n; // we will get it to be k=n-d by calling shrinkstar d times below randomX = calloc(n*d,sizeof(float)); randomXcopy = calloc(n*d,sizeof(float)); @@ -77,7 +81,6 @@ void randomstabilizerstate(int n, int* k, int** h, int*** G, int*** GBar, int* Q } *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 "randomstabilizerstate_equatorial.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_equatorial(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J) +{ + // vector and matrix pointers should all be passed unallocated! (no freeing of memory performed) + + int d; + + int i, j; + + d = 0; + + *k = n; // we will get it to be k=n-d by calling shrinkstar d times below + + int info; + int numsingvals = -1; + + *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 "shrinkstar.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 shrinkstar(int n, int *k, int *h, int **G, int **GBar, 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/shrinkstar.h b/shrinkstar.h new file mode 100644 index 0000000..88aba74 --- /dev/null +++ b/shrinkstar.h @@ -0,0 +1 @@ +char shrinkstar(int n, int *k, int *h, int **G, int **GBar, int *xi, int alpha); diff --git a/sparsify.c b/sparsify.c index 8e164fc..268dd98 100644 --- a/sparsify.c +++ b/sparsify.c @@ -5,6 +5,7 @@ #include "matrix.h" #include "supplement.h" +#include "supplement2.h" #include "sparsify.h" /**************************************************************************** @@ -24,32 +25,79 @@ // 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 sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor) { 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)); + // fraction of optimal T supplemental states + //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/((double)T); + //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/(4.0+pow(delta,2.0))/((double)T); + double alpha; + + // fraction of T supplemental states that we actually need for optimal scaling + // if first term is kept from Seddon et al. + //double fractionSupplement = alpha*0.25*delta*pow(L1Norm,2.0)/((double)T); + // if first term isn't kept + double fractionSupplement; + printf("delta = %lf\n", delta); - //printf("L1Norm = %lf\n", L1Norm); + printf("L1Norm = %lf\n", L1Norm); double gamma; int increment; - if(coherentSampling) { + if(coherentSampling!=0) { + if(coherentSampling==2) { + alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(2*T-1)); + fractionSupplement = alpha*1.0; + *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(2*T-1))+1)))*(floor(fractionSupplement*(2*T-1))+1); + increment = floor(fractionSupplement*(2*T-1))+1; + } else if(coherentSampling==1) { + alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(T)); + fractionSupplement = alpha*1.0; + *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(T))+1)))*(floor(fractionSupplement*(T))+1); + increment = floor(fractionSupplement*(T))+1; + } + 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 { + printf("Using an optimal fraction of T supplemental states!\n"); + printf("Optimal fraction = %lf\n", fractionSupplement); + // since we will be taking a multiple of (fractionSupplement*T+1) samples, we round numStabStates up to the nearest multiple of (fractionSupplement*T+1) + // if first term is kept from Seddon et al. + //*numStabStates = ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + // if first term isn't kept + //*numStabStates = ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + //*numStabStates = ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1); + + printf("increment= %d\n", increment); + //printf("prebinned numStabStates=%d\n", (int)ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0))); + //printf("prebinned numStabStates=%d\n", (int)ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0)))); + //printf("prebinned numStabStates=%d\n", (int)ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0)))); + printf("prebinned numStabStates=%d\n", (int)ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta)); + fflush(stdout); + if(fractionSupplement > 1.0) { + printf("Too many supplemental states requested!\n"); + printf("Note: the correlated sampling is currently only implemented for a maximum of (2t-1) supplemental correlated states!\n"); + return 1; + } + } else if(coherentSampling==0) { //gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T; + // SPARSIFY //*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)); + // if first term is kept from Seddon et al. + //*numStabStates = ceil((2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0)); + // if first term isn't kept + *numStabStates = ceil(samplePrefactor*(1.0*2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0)); + // multiply the factor of 2.0 by 1.0 to get the prior state-of-the-art + // multiply the factor of 2.0 by 0.0 to get the tighter bound on the prior state-of-the-art increment = 1; + } else { + printf("Error: Unknown value for \"coherent sampling\" flag"); + return 1; } @@ -108,10 +156,14 @@ int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, dou } (*stabStateIndices)[i] = index; //printf("first index = %d\n", index); - if(coherentSampling) - for(k=1; k<=T; k++) + if(coherentSampling==2) + //for(k=1; k<=T; k++) + //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index + supplement2(&(*stabStateIndices)[i], T, fractionSupplement); + else if(coherentSampling==1) + //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); + supplement(&(*stabStateIndices)[i], T, fractionSupplement); break; } //printf("%lf\n", sum); diff --git a/sparsify.h b/sparsify.h new file mode 100644 index 0000000..c13cff7 --- /dev/null +++ b/sparsify.h @@ -0,0 +1,4 @@ +#define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 + +double binomcoeff(int m, int n); +int sparsify(long** stabStateIndices, int* numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor); diff --git a/supplement.c b/supplement.c index f74574c..b3a5cbe 100644 --- a/supplement.c +++ b/supplement.c @@ -13,7 +13,7 @@ void printBinary(unsigned long number, int T) { } -void supplement(long stabStateIndices[], int T) +void supplement(long stabStateIndices[], int T, double fractionSupplement) { int K = round(log(T)/log(2)); // T = 2^K @@ -31,9 +31,12 @@ void supplement(long stabStateIndices[], int T) unsigned long bitstring; int bitstringCounter = 1; // bitstringCounter starts at 1 since 0 is the given bitstring + int bitstringCounterMax = floor(fractionSupplement*T); int i; + if(bitstringCounter > bitstringCounterMax) return; + // first added bitstring is alpha ... alpha for alpha = 01 bitstring = 1; for(i=1; i bitstringCounterMax) return; // second added bitstring is beta ... beta for beta = 10 bitstring = 2; @@ -53,6 +57,7 @@ void supplement(long stabStateIndices[], int T) //printf("%lu\n", bitstring); stabStateIndices[bitstringCounter] = given^bitstring; bitstringCounter++; + if(bitstringCounter > bitstringCounterMax) return; maskList[0] = round(pow(2,pow(2,K)))-1; @@ -78,6 +83,7 @@ void supplement(long stabStateIndices[], int T) //printf("new bitstring="); //printBinary(stabStateIndices[bitstringCounter], T); bitstringCounter++; + if(bitstringCounter > bitstringCounterMax) return; } maskList[treeDepth] = levelMask; @@ -97,7 +103,7 @@ int main(int argc, char *argv[]) stabStateIndices[0] = (int)(pow(2,T))-1; - supplement(&stabStateIndices[0], T); + supplement(&stabStateIndices[0], T, 1.0); int i; for(i=0; i +#include +#include +#include + +void printBinary(unsigned long number, int T) { + + int i; + + for(i=0; i1 will lead to more pessimistic runtime analysis, due to parallelization overhead, though the runs should complete faster! + +## For the most accurate time analysis, run with +## export OMP_NUM_THREADS=1 +## export OMP_PROC_BIND=false +#export OMP_NUM_THREADS=6 +#export OMP_PROC_BIND=true +export OMP_NUM_THREADS=1 +export OMP_PROC_BIND=false + +delta=0.1 +deltasquared=`bc -l <<< $delta^2` + +numruns=10 #1000 + +#for t in 4 8 +for t in 8 +#for t in 16 +#for t in 32 +#for t in 2 4 8 #16 #32 +#for t in 8 16 32 +do + #delta=`bc -l <<< 10^-$deltapower` + echo "t=$t" + # alpha is the fraction of t supplemental states that we will take to keep the cost the same + # alpha = delta*xi^t/t + #echo "0.25*$delta*1.17^$t/$t" + #alpha=`bc -l <<< 0.25*$delta*1.17^$t/$t` + #%echo "alpha=$alpha" + TIMEFORMAT='%3R'; + #number of runs + for (( i=1; i<=$numruns; i++ )) + do + # Generate two t-Pauli measurements + ./randominputcommutingHermitianPauli2 $t $t 1000 | head -n $[2*$t+4] > inputPauli_$[i].txt + sleep 1 + done + for (( i=1; i<=$numruns; i++ )) + do + # Exact measurements + ./multipauli < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 1 > tmp_$[t]_exact_$[i].txt & + #seed=`od -A n -t d -N 4 /dev/urandom |tr -d ' '|tr -d '-'` + #./weaksim 0.001 0.25 2 $seed < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 3 > tmp_$[t]_exact_$[i].txt & + #./weaksim 0.01 0.25 2 $seed < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 3 > tmp_$[t]_exact_$[i].txt & + # minus one below so as not to count grep search job + numthreads=$[`ps -aux | grep multipauli | wc -l`-1] + #numthreads=$[`ps -aux | grep weaksim | wc -l`-1] + while [ $numthreads -gt $[$maxthreads-1] ] + do + sleep 10 + numthreads=$[`ps -aux | grep multipauli | wc -l`-1] + #numthreads=$[`ps -aux | grep weaksim | wc -l`-1] + done + done + + for (( i=1; i<=$numruns; i++ )) + do + echo "run: $i" + + approxerror=1.0 + approxcoherror=1.0 + + samplePrefactor=0.001 + seed=`od -A n -t d -N 4 /dev/urandom |tr -d ' '|tr -d '-'` + while (( $(echo "$approxerror > $deltasquared" |bc -l) )); do + echo "samplePrefactor=$samplePrefactor" + ( time ( ./weaksim $delta 0.25 0 $seed $samplePrefactor < inputPauli_$[i].txt >> tmp_$[t]_iid_$[i].txt ) ) 2> tmp_$[t]_iid_$[i]_timing.txt + + exactanswer=`tail -1 tmp_$[t]_exact_$[i].txt` + approxanswer=`tail -1 tmp_$[t]_iid_$[i].txt | cut -d ' ' -f 3` + approxerror=`bc -l <<< "sqrt(($exactanswer-$approxanswer)^2)"` + echo $approxerror + samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor*1.1)/1"` + done + samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor/1.1)/1"` + echo "samplePrefactor used=$samplePrefactor" + echo $approxerror >> tmp_$[t]_iid_$[i]_error.txt + + samplePrefactor=0.001 + while (( $(echo "$approxcoherror > $deltasquared" |bc -l) )); do + echo "samplePrefactor=$samplePrefactor" + ( time ( ./weaksim $delta 0.25 2 $seed $samplePrefactor < inputPauli_$[i].txt >> tmp_$[t]_t_$[i].txt ) ) 2> tmp_$[t]_t_$[i]_timing.txt + + exactanswer=`tail -1 tmp_$[t]_exact_$[i].txt` + approxcoherentanswer=`tail -1 tmp_$[t]_t_$[i].txt | cut -d ' ' -f 3` + approxcoherror=`bc -l <<< "sqrt(($exactanswer-$approxcoherentanswer)^2)"` + echo $approxcoherror + samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor*1.1)/1"` + done + samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor/1.1)/1"` + echo "samplePrefactor used=$samplePrefactor" + echo $approxcoherror >> tmp_$[t]_t_$[i]_error.txt + + done + + for (( i=1; i<=$numruns; i++ )) + do + rm tmp_$[t]_exact_$[i].txt; + rm tmp_$[t]_iid_$[i].txt; + rm tmp_$[t]_t_$[i].txt; + rm inputPauli_$[i].txt; + #echo "*******************" >> tmp_$[t]_iid.txt + #echo "*******************" >> tmp_$[t]_t.txt + # sleep 1 + done + echo "Longest runtime:" >> tmp_$[t]_iid.txt + cat tmp_$[t]_iid_[1-9]*_timing.txt | sort -n | tail -2 >> tmp_$[t]_iid.txt + rm tmp_$[t]_iid_[1-9]*_timing.txt + echo "Longest runtime:" >> tmp_$[t]_t.txt + cat tmp_$[t]_t_[1-9]*_timing.txt | sort -n | tail -2 >> tmp_$[t]_t.txt + rm tmp_$[t]_t_[1-9]*_timing.txt + echo "Largest error:" >> tmp_$[t]_iid.txt + cat tmp_$[t]_iid_[1-9]*_error.txt | sort -n | tail -1 >> tmp_$[t]_iid.txt + rm tmp_$[t]_iid_[1-9]*_error.txt + echo "Largest error:" >> tmp_$[t]_t.txt + cat tmp_$[t]_t_[1-9]*_error.txt | sort -n | tail -1 >> tmp_$[t]_t.txt + rm tmp_$[t]_t_[1-9]*_error.txt + echo "delta=$delta done!" +done -- 2.30.2 From d0e2a9dbedede81c404357d1e381a1a5564ea516 Mon Sep 17 00:00:00 2001 From: Lucas Date: Tue, 1 Feb 2022 16:36:46 -0800 Subject: [PATCH 3/3] final v1.0 files --- weaksim.c | 82 ++++++++++--- weaksim_relerr.c | 298 +++++++++++++++++++++++------------------------ 2 files changed, 214 insertions(+), 166 deletions(-) diff --git a/weaksim.c b/weaksim.c index 8ec6222..b2c93f4 100644 --- a/weaksim.c +++ b/weaksim.c @@ -11,6 +11,7 @@ #include "measurepauli.h" #include "innerproduct.h" #include "sparsify.h" +#include "omp.h" int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits); @@ -19,14 +20,18 @@ int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, i 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"); + if(argc != 5 && argc != 6) { + printf("weaksim argument: \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\" \"seed\" \"(optional) sample number prefactor\"\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 coherentSampling = atoi(argv[3]); // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3) + int seed = atoi(argv[4]); // random number seed + double samplePrefactor = 1.0; + if(argc == 6) + samplePrefactor = atof(argv[5]); // factor that multiplies k, the number of samples used int N; // number of qubits scanf("%d", &N); @@ -57,16 +62,17 @@ int main(int argc, char *argv[]) long* stabStateIndices; int numStabStates; - srand((unsigned)time(NULL)); // seeding the random number generator for sparsify() + //srand((unsigned)time(NULL)); // seeding the random number generator for sparsify() + srand((unsigned)seed); // seeding the random number generator for sparsify() - if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling)) + if(sparsify(&stabStateIndices, &numStabStates, T, phi, additiveErrorDelta, coherentSampling, samplePrefactor)) return 1; - //printf("checking: numStabStateIndices:\n"); - //for(i=0; i0?'+':'-' , cabs(cimag(lastamplitude))); + printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude))); + + probability *= amplitude/lastamplitude; + //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement + + printf("probability = %lf\n", cabs(probability)); + printf("!\n"); } //printf("numStabStates=%d\n", numStabStates); @@ -295,6 +341,8 @@ int main(int argc, char *argv[]) printf("\nabs(amplitude):\n"); printf("%lf\n", cabs(amplitude)); + printf("probability = %lf\n", cabs(probability)); + return 0; } diff --git a/weaksim_relerr.c b/weaksim_relerr.c index 2cc45fb..09bf002 100644 --- a/weaksim_relerr.c +++ b/weaksim_relerr.c @@ -6,7 +6,7 @@ #include #include "matrix.h" #include "exponentialsum.h" -#include "shrinkstar.h" +#include "shrink.h" #include "extend.h" #include "measurepauli.h" #include "innerproduct.h" @@ -23,25 +23,25 @@ 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"); + printf("weaksim_relerr argument: \"number of stabilizer state samples\" \"additive error delta\" \"phi (times PI)\" \"coherent sampling (0=no; 1=t; 2=2t-1; 3=xi^3t/2)\\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 coherentSampling = atoi(argv[4]); // perform coherent sampling (false=0; linear t=1; linear 2t-1=2; exponential xi^3t/2=3) - int N; // number of qubits + 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); + 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 omega; + int alpha[N], beta[N], gamma[N], delta[N]; int Paulicounter = 0; int i, j, k, l, m; @@ -49,7 +49,7 @@ int main(int argc, char *argv[]) FILE *fp; char buff[255]; - srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate() + srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate _equatorial() fp = fopen("Pd.txt", "r"); @@ -63,7 +63,7 @@ int main(int argc, char *argv[]) int PdN = atoi(buff); Pd = calloc(PdN, sizeof(double*)); - for(i=0; i = 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 @@ -98,15 +100,14 @@ int main(int argc, char *argv[]) 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) { + // origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int)); + // for(k=0; k N) { + while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) { + + Paulicounter++; + if(Paulicounter > N) { printf("Error: Number of Paulis is greater than N!\n"); return 1; } @@ -287,116 +342,61 @@ int main(int argc, char *argv[]) // 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?'+':'-' , cabs(cimag(lastamplitude))); + printf("amplitude: %lf %c %lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude))); + + probability *= amplitude/lastamplitude; + //amplitude = amplitude/lastamplitude; // for NORMALIZE-SPARSIFY you normalize the algorithm after every Pauli measurement - deallocate_mem(&origG, N); - deallocate_mem(&origGBar, N); - free(origh); - deallocate_mem(&origJ, origK); - free(origD); + printf("probability = %lf\n", cabs(probability)); + printf("!\n"); } //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))); + 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("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude))); + 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)); - - for(i=0; i