From 7bc25656311c391af6a4d5814cc82080981c63c4 Mon Sep 17 00:00:00 2001 From: Lucas Date: Tue, 1 Feb 2022 16:35:48 -0800 Subject: [PATCH] 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