From f36d60f8c85e0879a7d2b52e816638d2b4fdf86a Mon Sep 17 00:00:00 2001 From: Lucas K Date: Tue, 22 Dec 2020 14:54:55 -0800 Subject: [PATCH] added core code, README, and testing scripts --- Makefile | 69 + Pd.txt | 51 + README.txt | 56 + exponentialsum.c | 415 +++++ exponentialsum.h | 3 + extend.c | 177 +++ extend.h | 1 + innerproduct.c | 190 +++ innerproduct.h | 1 + matrix.c | 243 +++ matrix.h | 18 + measurepauli.c | 165 ++ measurepauli.h | 2 + multipauli.c | 301 ++++ randominputPauli.c | 37 + randominputcommutingHermitianPauli.c | 58 + randominputcommutingHermitianPauli2.c | 109 ++ randomstabilizerstate.c | 131 ++ randomstabilizerstate.h | 1 + shrink.c | 380 +++++ shrink.h | 1 + shrinkstar.c | 380 +++++ shrinkstar.h | 1 + strongsim.c | 286 ++++ strongsim1.c | 186 +++ strongsim12.c | 2013 ++++++++++++++++++++++++ strongsim12_relerr.c | 2059 +++++++++++++++++++++++++ strongsim1_relerr.c | 409 +++++ strongsim2.c | 210 +++ strongsim2_relerr.c | 396 +++++ strongsim3.c | 286 ++++ strongsim3_relerr.c | 393 +++++ strongsim6.c | 378 +++++ strongsim6_relerr.c | 432 ++++++ strongsim_relerr.c | 393 +++++ test.bash | 30 + test2.bash | 42 + timing.bash | 11 + 38 files changed, 10314 insertions(+) create mode 100644 Makefile create mode 100644 Pd.txt create mode 100644 README.txt create mode 100644 exponentialsum.c create mode 100644 exponentialsum.h create mode 100644 extend.c create mode 100644 extend.h create mode 100644 innerproduct.c create mode 100644 innerproduct.h create mode 100644 matrix.c create mode 100644 matrix.h create mode 100644 measurepauli.c create mode 100644 measurepauli.h create mode 100644 multipauli.c create mode 100644 randominputPauli.c create mode 100644 randominputcommutingHermitianPauli.c create mode 100644 randominputcommutingHermitianPauli2.c create mode 100644 randomstabilizerstate.c create mode 100644 randomstabilizerstate.h create mode 100644 shrink.c create mode 100644 shrink.h create mode 100644 shrinkstar.c create mode 100644 shrinkstar.h create mode 100644 strongsim.c create mode 100644 strongsim1.c create mode 100644 strongsim12.c create mode 100644 strongsim12_relerr.c create mode 100644 strongsim1_relerr.c create mode 100644 strongsim2.c create mode 100644 strongsim2_relerr.c create mode 100644 strongsim3.c create mode 100644 strongsim3_relerr.c create mode 100644 strongsim6.c create mode 100644 strongsim6_relerr.c create mode 100644 strongsim_relerr.c create mode 100644 test.bash create mode 100644 test2.bash create mode 100644 timing.bash diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..19865ad --- /dev/null +++ b/Makefile @@ -0,0 +1,69 @@ +#IDIR =../include +CC=gcc -std=c99 +CFLAGS=-Wall +LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.so libinnerproduct.so + +strongsim_relerr: strongsim_relerr.c matrix exponentialsum shrink shrinkstar extend measurepauli innerproduct randomstabilizerstate + $(CC) -o $@ strongsim_relerr.c $(CFLAGS) $(LIBS) libshrink.so libshrinkstar.so librandomstabilizerstate.so + +strongsim: strongsim.c matrix exponentialsum shrink extend measurepauli innerproduct + $(CC) -o $@ strongsim.c $(CFLAGS) $(LIBS) libshrink.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 + +randominputcommutingHermitianPauli: randominputcommutingHermitianPauli.c + $(CC) -o randominputcommutingHermitianPauli randominputcommutingHermitianPauli.c + +randominputPauli: randominputPauli.c + $(CC) -o randominputPauli randominputPauli.c + +multipauli: multipauli.c + $(CC) -o multipauli multipauli.c -lm + +.PHONY: clean + +clean: + rm ./strongsim ./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/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..6a53e3a --- /dev/null +++ b/README.txt @@ -0,0 +1,56 @@ +Strong simulation of multi-Pauli measurements using minimal stabilizer state decompositions +------------------------------------------------------------------------------------------- + +BUILD: + +1. We start with the code that scales quadratically with stabilizer rank because it calculates expectation values to zero relative error. Build code for stabilizer rank tensor multiple you'd like to use. Here we choose k=6. +$ cp strongsim6.c strongsim.c +$ make strongsim + +2. Now we build the code that scales linearly with stabilizer rank because it calculates expectation values to non-zero relative error. Build code for stabilizer rank tensor multiple you'd like to use. Here we choose k=6. +$ cp strongsim6_relerr.c strongsim_relerr.c +$ make strongsim_relerr + +3. 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. 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 + +5. Build code for generating random (non-commuting) Pauli inputs (good enough for the strongsim[1-12].c code). +$ make randominputFauli + +6. Build Hilbert vector space code to check our results. +$ make multipauli + + +TEST CODE: +1. Run test diagonositic for strongsim.c (zero relative error). +$ bash ./test.bash + +2. Run test diagonositic for strongsim_relerr.c (non-zero relative error). +$ bash ./test2.bash + + + +Details of source code files: + +strongsim[1-12]_relerror.c takes one argument of the number of stabilizer states to sample and then takes stdin arguments for the number of qubits, t-gate magic states and then the Paulis + +strongsim[1-12]_relerror.c scale linearly with the stabilizer rank and calculate up to finite relative error while strongsim[1-12].c scale quadratically with the stabilizer rank and calculate up to zero relative error. + +strongsim*c need Pauli measurement inputs (with the associated phases) which can be randomly generated with randominputPauli.c and the output from strongsim*c can be verified with multipauli.c + +strongsim[1-12]_relerror.c need commuting Pauli measurement inputs (with the associated phasese) which can be randomly generated with randominputcommutingHermitianPauli2.c and can be verified with multipauli.c + +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 new file mode 100644 index 0000000..0e0bd6e --- /dev/null +++ b/exponentialsum.c @@ -0,0 +1,415 @@ +#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) { + free(setWalker->next); + setWalker->next = setWalker->next->next; // delete 'b' + 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/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.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/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.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/measurepauli.h b/measurepauli.h new file mode 100644 index 0000000..5405627 --- /dev/null +++ b/measurepauli.h @@ -0,0 +1,2 @@ + +double measurepauli(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int m, int *X, int *Z); diff --git a/multipauli.c b/multipauli.c new file mode 100644 index 0000000..2161b72 --- /dev/null +++ b/multipauli.c @@ -0,0 +1,301 @@ +#include +#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/randominputPauli.c b/randominputPauli.c new file mode 100644 index 0000000..da17abe --- /dev/null +++ b/randominputPauli.c @@ -0,0 +1,37 @@ +#include +#include +#include + +// order of matrix elements is [row][column]!!! + +int main( int argc, char *argv[]) +{ + + if(argc != 3) { + printf("randominputPauli arguments: \"number of qubits\" \"number of T gates\"\n"); + exit(0); + } + + int N = atoi(argv[1]); // number of qubits + int T = atoi(argv[2]); // 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) + + printf("%d\n", N); + printf("%d\n", T); + + int r; + srand((unsigned)time(NULL)); + + + int i, j; + + for(i=0; i +#include +#include +#include + +// order of matrix elements is [row][column]!!! + +/*********************************************************/ +/* Generates random Hermitian Paulis that are commuting! */ +/*********************************************************/ +int main( int argc, char *argv[]) +{ + + if(argc != 3) { + printf("randominputPauli arguments: \"number of qubits\" \"number of T gates\"\n"); + exit(0); + } + + int N = atoi(argv[1]); // number of qubits + int T = atoi(argv[2]); // 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) + + printf("%d\n", N); + printf("%d\n", T); + + int Pauli[N][N]; + srand((unsigned)time(NULL)); + + + int i, j, k; + + int sign = 1; + + for(i=0; i +#include +#include +#include + +// order of matrix elements is [row][column]!!! + +/*********************************************************/ +/* Generates random Hermitian Paulis that are commuting! */ +/*********************************************************/ +int main( int argc, char *argv[]) +{ + + if(argc != 4) { + printf("randominputPauli arguments: \"number of qubits\" \"number of T gates\" \"number of random Clifford steps\"\n"); + exit(0); + } + + int N = atoi(argv[1]); // number of qubits + int T = atoi(argv[2]); // 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) + int randsteps = atoi(argv[3]); // number of random Clifford steps + + printf("%d\n", N); + printf("%d\n", T); + + int omega[N]; + int Pauli[N][N]; + srand((unsigned)time(NULL)); + + + int i, j, k; + + int sign = 1; + + // don't initially generate Y's (which is represented by '3'; I=0, Z=1, X=2, Y=3) + for(i=0; i2 and 2->1 (i.e. Z->X and X->Z) + else if(Pauli[j][qubita] ==2) + Pauli[j][qubita] = 1; // 1->2 and 2->1 (i.e. Z->X and X->Z) + } + } else if(gate==1) { // phaseshift + qubita = rand()%N; + //printf("qubita=%d", qubita); + for(j=0; jY + else if(Pauli[j][qubita] == 3) { + Pauli[j][qubita] = 2; // Y->-X + omega[j] = (2*omega[j])%2; + } + } + } else { // CNOT + qubita = rand()%N; + qubitb = qubita; + while(qubitb == qubita) + qubitb = rand()%N; + for(j=0; jZI + Pauli[j][qubitb] == 0; + else if((Pauli[j][qubita] == 1) && (Pauli[j][qubitb] == 0)) // ZI ->ZZ + Pauli[j][qubitb] == 1; + else if((Pauli[j][qubita] == 2) && (Pauli[j][qubitb] == 2)) // XX ->IX + Pauli[j][qubita] == 1; + else if((Pauli[j][qubita] == 0) && (Pauli[j][qubitb] == 2)) // IX ->XX + Pauli[j][qubita] == 2; + else if((Pauli[j][qubita] == 1) && (Pauli[j][qubitb] == 2)) {// ZX ->YY + Pauli[j][qubita] == 3; + Pauli[j][qubitb] == 3; + } else if((Pauli[j][qubita] == 0) && (Pauli[j][qubitb] == 2)) {// YY ->ZX + Pauli[j][qubita] == 1; + Pauli[j][qubitb] == 2; + } else if((Pauli[j][qubita] == 3) && (Pauli[j][qubitb] == 2)) {// YX ->-ZY + Pauli[j][qubita] == 1; + Pauli[j][qubitb] == 3; + omega[j] = -omega[j]; + } else if((Pauli[j][qubita] == 1) && (Pauli[j][qubitb] == 3)) {// ZY ->-YX + Pauli[j][qubita] == 3; + Pauli[j][qubitb] == 2; + omega[j] = -omega[j]; + } + } + } + } + + for(i=0; i +#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/shrink.h b/shrink.h new file mode 100644 index 0000000..543b20c --- /dev/null +++ b/shrink.h @@ -0,0 +1 @@ +char shrink(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int *xi, int alpha); diff --git a/shrinkstar.c b/shrinkstar.c new file mode 100644 index 0000000..ea202f8 --- /dev/null +++ b/shrinkstar.c @@ -0,0 +1,380 @@ +#include +#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/strongsim.c b/strongsim.c new file mode 100644 index 0000000..389a737 --- /dev/null +++ b/strongsim.c @@ -0,0 +1,286 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.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 N; // number of qubits + scanf("%d", &N); + + if(N%3 != 0) { + printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n"); + return 1; + } + + 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l, m; + + double complex coeffa = -0.25*(1.0-I)*(-1.0-I+sqrt(2.0))*csqrt(-I); + double complex coeffb = 0.25*(-1.0-I)*(1.0-I+sqrt(2.0))*csqrt(I); + double complex coeffc = 0.25*(-1.0-I)*(-1.0+I+sqrt(2.0))*csqrt(I); + + int n1 = 3; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1, 1}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0}, (int[]) {1, 1, 0}, (int[]) {1, 0, 1}}; int h1[] = {1, 1, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} }; + int n2 = 3; int k2 = 3; int (*(G2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h2[] = {0, 0, 0}; int Q2 = 2; int D2[] = {2, 2, 0}; int (*(J2[])) = { (int[]) {4, 0, 0}, (int[]) {0, 4, 0}, (int[]) {0, 0, 0} }; + int n3 = 3; int k3 = 3; int (*(G3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h3[] = {0, 0, 0}; int Q3 = 2; int D3[] = {6, 6, 0}; int (*(J3[])) = { (int[]) {4, 4, 4}, (int[]) {4, 4, 4}, (int[]) {4, 4, 0} }; + + int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J; + double complex Gamma[(int)pow(3,N/3)]; // prefactor in front of resultant state + G = calloc(pow(3,N/3),sizeof(int*)); GBar = calloc(pow(3,N/3),sizeof(int*)); + h = calloc(pow(3,N/3),sizeof(int*)); + + J = calloc(pow(3,N/3),sizeof(int*)); D = calloc(pow(3,N/3),sizeof(int*)); Q = calloc(pow(3,N/3),sizeof(int)); + + K = calloc(pow(3,N/3), sizeof(int)); + + double complex origGamma[(int)pow(3,N/3)]; + int *origK, *origQ, **origD, ***origJ; + int ***origG, ***origGBar, **origh; + + origG = calloc(pow(3,N/3),sizeof(int*)); origGBar = calloc(pow(3,N/3),sizeof(int*)); + origh = calloc(pow(3,N/3),sizeof(int*)); + + origJ = calloc(pow(3,N/3),sizeof(int*)); origD = calloc(pow(3,N/3),sizeof(int*)); origQ = calloc(pow(3,N/3),sizeof(int)); + + origK = calloc(pow(3,N/3), sizeof(int)); + + int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together + + + for(j=0; j 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))); + + 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/strongsim1.c b/strongsim1.c new file mode 100644 index 0000000..2277799 --- /dev/null +++ b/strongsim1.c @@ -0,0 +1,186 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.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 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l; + + + int n1 = 1; int k1 = 0; int (*(G1[])) = { (int[]) {1} }; int (*(GBar1[])) = { (int[]) {1} }; int h1[] = {0}; int Q1 = 0; int D1[] = {0}; int (*(J1[])) = { (int[]) {0} }; + int n2 = 1; int k2 = 0; int (*(G2[])) = { (int[]) {1} }; int (*(GBar2[])) = { (int[]) {1} }; int h2[] = {1}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} }; + + int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J; + double complex Gamma[(int)pow(2,N)]; // prefactor in front of resultant state + G = calloc(pow(2,N),sizeof(int*)); GBar = calloc(pow(2,N),sizeof(int*)); + h = calloc(pow(2,N),sizeof(int*)); + + J = calloc(pow(2,N),sizeof(int*)); D = calloc(pow(2,N),sizeof(int*)); Q = calloc(pow(2,N),sizeof(int)); + + K = calloc(pow(2,N), sizeof(int)); + + double complex origGamma[(int)pow(2,N)]; + int *origK, *origQ, **origD, ***origJ; + int ***origG, ***origGBar, **origh; + + origG = calloc(pow(2,N),sizeof(int*)); origGBar = calloc(pow(2,N),sizeof(int*)); + origh = calloc(pow(2,N),sizeof(int*)); + + origJ = calloc(pow(2,N),sizeof(int*)); origD = calloc(pow(2,N),sizeof(int*)); origQ = calloc(pow(2,N),sizeof(int)); + + origK = calloc(pow(2,N), sizeof(int)); + + int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together + + + for(j=0; j 0) { + + J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int)); + } + + origG[j] = calloc(N, sizeof(int*)); origGBar[j] = calloc(N, sizeof(int*)); + origh[j] = calloc(N, sizeof(int)); + + if(K[j] > 0) { + origJ[j] = calloc(K[j], sizeof(int*)); origD[j] = calloc(K[j], sizeof(int)); + } + + for(k=0; k 0) + J[j] = calloc(K[j], sizeof(int*)); + + origG[j][k] = calloc(N, sizeof(int*)); origGBar[j][k] = calloc(N, sizeof(int*)); + + if(K[j] > 0) + origJ[j] = calloc(K[j], sizeof(int*)); + + h[j][k] = (combination & 0x1)*h2[0] + (~combination & 0x1)*h1[0]; + for(l=0; l 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?'+':'-' , cabs(cimag(amplitude))); + else + printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(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/strongsim12.c b/strongsim12.c new file mode 100644 index 0000000..1db1250 --- /dev/null +++ b/strongsim12.c @@ -0,0 +1,2013 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.h" + +#define ZEROTHRESHOLD (0.00000001) + +int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits); + +void outerProductMod(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2, int mod); + +// order of matrix elements is [row][column]!!! + +int main() +{ + + int N; // number of qubits + scanf("%d", &N); + + if(N%12 != 0) { + printf("'N' needs to be a multiple of 12 for a k=12 tensor factor decomposition!\n"); + return 1; + } + + 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l, m; + + double complex coeffb60 = (-16.0+12.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0); + double complex coeffb66 = (96.0-68.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0); + double complex coeffe6 = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffo6 = (-14.0+10.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-14.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffk6 = (7.0-5.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-8.0*PI*I/8.0)*4.0*csqrt(2.0)*cpow(2.0,0.5); + double complex coeffphiprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffphidprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5); + + // b60 + int nn1 = 6; int kk1 = 6; int (*(GG1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GGBar1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int hh1[] = {0, 0, 0, 0, 0, 0}; int QQ1 = 0; int DD1[] = {0, 0, 0, 0, 0, 0}; int (*(JJ1[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} }; + // b66 + int nn2 = 6; int kk2 = 6; int (*(GG2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GGBar2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int hh2[] = {0, 0, 0, 0, 0, 0}; int QQ2 = 4; int DD2[] = {4, 4, 4, 4, 4, 4}; int (*(JJ2[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} }; + // e6 + int nn3 = 6; int kk3 = 5; int (*(GG3[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar3[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh3[] = {1, 0, 0, 0, 0, 0}; int QQ3 = 4; int DD3[] = {0, 0, 0, 0, 0}; int (*(JJ3[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} }; + // o6 + int nn4 = 6; int kk4 = 5; int (*(GG4[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar4[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh4[] = {0, 0, 0, 0, 0, 0}; int QQ4 = 4; int DD4[] = {4, 4, 4, 4, 4}; int (*(JJ4[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} }; + // k6 + int nn5 = 6; int kk5 = 1; int (*(GG5[])) = { (int[]) {1, 1, 1, 1, 1, 1}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1} }; int (*(GGBar5[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1} }; int hh5[] = {1, 1, 1, 1, 1, 1}; int QQ5 = 6; int DD5[] = {2}; int (*(JJ5[])) = { (int[]) {4} }; + // phiprime + int nn6 = 6; int kk6 = 5; int (*(GG6[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar6[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh6[] = {1, 0, 0, 0, 0, 0}; int QQ6 = 0; int DD6[] = {0, 0, 0, 0, 0}; int (*(JJ6[])) = { (int[]) {0, 4, 0, 0, 4}, (int[]) {4, 0, 4, 0, 0}, (int[]) {0, 4, 0, 4, 0}, (int[]) {0, 0, 4, 0, 4}, (int[]) {4, 0, 0, 4, 0} }; + // phidoubleprime + int nn7 = 6; int kk7 = 5; int (*(GG7[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar7[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh7[] = {1, 0, 0, 0, 0, 0}; int QQ7 = 0; int DD7[] = {0, 0, 0, 0, 0}; int (*(JJ7[])) = { (int[]) {0, 0, 4, 4, 0}, (int[]) {0, 0, 0, 4, 4}, (int[]) {4, 0, 0, 0, 4}, (int[]) {4, 4, 0, 0, 0}, (int[]) {0, 4, 4, 0, 0} }; + + // b60 * b60 + int n1 = nn1+nn1; int k1 = kk1+kk1; int **G1; int **GBar1; int *h1; int Q1; int *D1; int **J1; + G1 = calloc(nn1+nn1,sizeof(int*)); + GBar1 = calloc(nn1+nn1,sizeof(int*)); + h1 = calloc(nn1+nn1,sizeof(int)); + D1 = calloc(kk1+kk1,sizeof(int)); + J1 = calloc(kk1+kk1,sizeof(int*)); + 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("%.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))); + + 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/strongsim12_relerr.c b/strongsim12_relerr.c new file mode 100644 index 0000000..8ec8d76 --- /dev/null +++ b/strongsim12_relerr.c @@ -0,0 +1,2059 @@ +#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" + +#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 != 2) { + printf("strongsim12_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + if(N%12 != 0) { + printf("'N' needs to be a multiple of 12 for a k=12 tensor factor decomposition!\n"); + return 1; + } + + 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); + + 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 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))); + + + + 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/strongsim1_relerr.c b/strongsim1_relerr.c new file mode 100644 index 0000000..33210a1 --- /dev/null +++ b/strongsim1_relerr.c @@ -0,0 +1,409 @@ +#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" + +#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 != 2) { + printf("strongsim2_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + //if(N%2 != 0) { + // printf("'N' needs to be a multiple of 2 for a k=2 tensor factor decomposition!\n"); + // return 1; + //} + + 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); + + 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 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; i=0; k--) { // go backwards through the measurements since acting on the bra (doesn't matter if they commute though) + origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]); + //printf("k=%d\n", k); + } + /*printf("origK=%d\n", origK); + printf("origG:\n"); + printMatrix(origG, N, N); + printf("origGBar:\n"); + printMatrix(origGBar, N, N); + printf("origh:\n"); + printVector(origh, N);*/ + + double complex stabstateaverage = 0.0 + 0.0*I; + + for(j=0; j0) + 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))); + + + + 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/strongsim2.c b/strongsim2.c new file mode 100644 index 0000000..3b9b36b --- /dev/null +++ b/strongsim2.c @@ -0,0 +1,210 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.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 N; // number of qubits + scanf("%d", &N); + + if(N%2 != 0) { + printf("'N' needs to be an even number for a k=2 tensor factor decomposition!\n"); + return 1; + } + + 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l, m; + + + int n1 = 2; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1}, (int[]) {0, 1} }; int (*(GBar1[])) = { (int[]) {1, 0}, (int[]) {1, 1} }; int h1[] = {0, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} }; + int n2 = 2; int k2 = 1; int (*(G2[])) = { (int[]) {1, 1}, (int[]) {0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0}, (int[]) {1, 1} }; int h2[] = {0, 1}; int Q2 = 0; int D2[] = {0}; int (*(J2[])) = { (int[]) {0} }; + + int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J; + double complex Gamma[(int)pow(2.0,N/2)]; // prefactor in front of resultant state + G = calloc(pow(2,N/2),sizeof(int*)); GBar = calloc(pow(2,N/2),sizeof(int*)); + h = calloc(pow(2,N/2),sizeof(int*)); + + J = calloc(pow(2,N/2),sizeof(int*)); D = calloc(pow(2,N/2),sizeof(int*)); Q = calloc(pow(2,N/2),sizeof(int)); + + K = calloc(pow(2,N/2), sizeof(int)); + + double complex origGamma[(int)pow(2,N/2)]; + int *origK, *origQ, **origD, ***origJ; + int ***origG, ***origGBar, **origh; + + origG = calloc(pow(2,N/2),sizeof(int*)); origGBar = calloc(pow(2,N/2),sizeof(int*)); + origh = calloc(pow(2,N/2),sizeof(int*)); + + origJ = calloc(pow(2,N/2),sizeof(int*)); origD = calloc(pow(2,N/2),sizeof(int*)); origQ = calloc(pow(2,N/2),sizeof(int)); + + origK = calloc(pow(2,N/2), sizeof(int)); + + int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together + + + for(j=0; j 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?'+':'-' , cabs(cimag(amplitude))); + else + printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(amplitude))); + //printf("%lf %c %lf I\n", creal(amplitude), cimag(amplitude)>0?'+':'-' , cabs(cimag(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/strongsim2_relerr.c b/strongsim2_relerr.c new file mode 100644 index 0000000..ca1f0b8 --- /dev/null +++ b/strongsim2_relerr.c @@ -0,0 +1,396 @@ +#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" + +#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 != 2) { + printf("strongsim2_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + if(N%2 != 0) { + printf("'N' needs to be a multiple of 2 for a k=2 tensor factor decomposition!\n"); + return 1; + } + + 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); + + 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 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; i=0; k--) { // go backwards through the measurements since acting on the bra (doesn't matter if they commute though) + origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]); + //printf("k=%d\n", k); + } + /*printf("origK=%d\n", origK); + printf("origG:\n"); + printMatrix(origG, N, N); + printf("origGBar:\n"); + printMatrix(origGBar, N, N); + printf("origh:\n"); + printVector(origh, N);*/ + + double complex stabstateaverage = 0.0 + 0.0*I; + + for(j=0; j0) + 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))); + + + + 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/strongsim3.c b/strongsim3.c new file mode 100644 index 0000000..389a737 --- /dev/null +++ b/strongsim3.c @@ -0,0 +1,286 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.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 N; // number of qubits + scanf("%d", &N); + + if(N%3 != 0) { + printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n"); + return 1; + } + + 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l, m; + + double complex coeffa = -0.25*(1.0-I)*(-1.0-I+sqrt(2.0))*csqrt(-I); + double complex coeffb = 0.25*(-1.0-I)*(1.0-I+sqrt(2.0))*csqrt(I); + double complex coeffc = 0.25*(-1.0-I)*(-1.0+I+sqrt(2.0))*csqrt(I); + + int n1 = 3; int k1 = 1; int (*(G1[])) = { (int[]) {1, 1, 1}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0}, (int[]) {1, 1, 0}, (int[]) {1, 0, 1}}; int h1[] = {1, 1, 0}; int Q1 = 0; int D1[] = {2}; int (*(J1[])) = { (int[]) {4} }; + int n2 = 3; int k2 = 3; int (*(G2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar2[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h2[] = {0, 0, 0}; int Q2 = 2; int D2[] = {2, 2, 0}; int (*(J2[])) = { (int[]) {4, 0, 0}, (int[]) {0, 4, 0}, (int[]) {0, 0, 0} }; + int n3 = 3; int k3 = 3; int (*(G3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int (*(GBar3[])) = { (int[]) {1, 0, 0}, (int[]) {0, 1, 0}, (int[]) {0, 0, 1} }; int h3[] = {0, 0, 0}; int Q3 = 2; int D3[] = {6, 6, 0}; int (*(J3[])) = { (int[]) {4, 4, 4}, (int[]) {4, 4, 4}, (int[]) {4, 4, 0} }; + + int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J; + double complex Gamma[(int)pow(3,N/3)]; // prefactor in front of resultant state + G = calloc(pow(3,N/3),sizeof(int*)); GBar = calloc(pow(3,N/3),sizeof(int*)); + h = calloc(pow(3,N/3),sizeof(int*)); + + J = calloc(pow(3,N/3),sizeof(int*)); D = calloc(pow(3,N/3),sizeof(int*)); Q = calloc(pow(3,N/3),sizeof(int)); + + K = calloc(pow(3,N/3), sizeof(int)); + + double complex origGamma[(int)pow(3,N/3)]; + int *origK, *origQ, **origD, ***origJ; + int ***origG, ***origGBar, **origh; + + origG = calloc(pow(3,N/3),sizeof(int*)); origGBar = calloc(pow(3,N/3),sizeof(int*)); + origh = calloc(pow(3,N/3),sizeof(int*)); + + origJ = calloc(pow(3,N/3),sizeof(int*)); origD = calloc(pow(3,N/3),sizeof(int*)); origQ = calloc(pow(3,N/3),sizeof(int)); + + origK = calloc(pow(3,N/3), sizeof(int)); + + int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together + + + for(j=0; j 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))); + + 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/strongsim3_relerr.c b/strongsim3_relerr.c new file mode 100644 index 0000000..24ec5b3 --- /dev/null +++ b/strongsim3_relerr.c @@ -0,0 +1,393 @@ +#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" + +#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 != 2) { + printf("strongsim3_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + if(N%3 != 0) { + printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n"); + return 1; + } + + 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); + + 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 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))); + + + + 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/strongsim6.c b/strongsim6.c new file mode 100644 index 0000000..65d73fe --- /dev/null +++ b/strongsim6.c @@ -0,0 +1,378 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "exponentialsum.h" +#include "shrink.h" +#include "extend.h" +#include "measurepauli.h" +#include "innerproduct.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 N; // number of qubits + scanf("%d", &N); + + if(N%6 != 0) { + printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n"); + return 1; + } + + 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); + + int omega; + int alpha[N], beta[N], gamma[N], delta[N]; + int Paulicounter = 0; + + int i, j, k, l, m; + + double complex coeffb60 = (-16.0+12.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0); + double complex coeffb66 = (96.0-68.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0); + double complex coeffe6 = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffo6 = (-14.0+10.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-14.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffk6 = (7.0-5.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-8.0*PI*I/8.0)*4.0*csqrt(2.0)*cpow(2.0,0.5); + double complex coeffphiprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5); + double complex coeffphidprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5); + + // b60 + int n1 = 6; int k1 = 6; int (*(G1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h1[] = {0, 0, 0, 0, 0, 0}; int Q1 = 0; int D1[] = {0, 0, 0, 0, 0, 0}; int (*(J1[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} }; + // b66 + int n2 = 6; int k2 = 6; int (*(G2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GBar2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int h2[] = {0, 0, 0, 0, 0, 0}; int Q2 = 4; int D2[] = {4, 4, 4, 4, 4, 4}; int (*(J2[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} }; + // e6 + int n3 = 6; int k3 = 5; int (*(G3[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar3[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h3[] = {1, 0, 0, 0, 0, 0}; int Q3 = 4; int D3[] = {0, 0, 0, 0, 0}; int (*(J3[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} }; + // o6 + int n4 = 6; int k4 = 5; int (*(G4[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar4[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h4[] = {0, 0, 0, 0, 0, 0}; int Q4 = 4; int D4[] = {4, 4, 4, 4, 4}; int (*(J4[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} }; + // k6 + int n5 = 6; int k5 = 1; int (*(G5[])) = { (int[]) {1, 1, 1, 1, 1, 1}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1} }; int (*(GBar5[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1} }; int h5[] = {1, 1, 1, 1, 1, 1}; int Q5 = 6; int D5[] = {2}; int (*(J5[])) = { (int[]) {4} }; + // phiprime + int n6 = 6; int k6 = 5; int (*(G6[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar6[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h6[] = {1, 0, 0, 0, 0, 0}; int Q6 = 0; int D6[] = {0, 0, 0, 0, 0}; int (*(J6[])) = { (int[]) {0, 4, 0, 0, 4}, (int[]) {4, 0, 4, 0, 0}, (int[]) {0, 4, 0, 4, 0}, (int[]) {0, 0, 4, 0, 4}, (int[]) {4, 0, 0, 4, 0} }; + // phidoubleprime + int n7 = 6; int k7 = 5; int (*(G7[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GBar7[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int h7[] = {1, 0, 0, 0, 0, 0}; int Q7 = 0; int D7[] = {0, 0, 0, 0, 0}; int (*(J7[])) = { (int[]) {0, 0, 4, 4, 0}, (int[]) {0, 0, 0, 4, 4}, (int[]) {4, 0, 0, 0, 4}, (int[]) {4, 4, 0, 0, 0}, (int[]) {0, 4, 4, 0, 0} }; + + int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J; + double complex Gamma[(int)pow(7,N/6)]; // prefactor in front of resultant state + G = calloc(pow(7,N/6),sizeof(int*)); GBar = calloc(pow(7,N/6),sizeof(int*)); + h = calloc(pow(7,N/6),sizeof(int*)); + + J = calloc(pow(7,N/6),sizeof(int*)); D = calloc(pow(7,N/6),sizeof(int*)); Q = calloc(pow(7,N/6),sizeof(int)); + + K = calloc(pow(7,N/6), sizeof(int)); + + double complex origGamma[(int)pow(7,N/6)]; + int *origK, *origQ, **origD, ***origJ; + int ***origG, ***origGBar, **origh; + + origG = calloc(pow(7,N/6),sizeof(int*)); origGBar = calloc(pow(7,N/6),sizeof(int*)); + origh = calloc(pow(7,N/6),sizeof(int*)); + + origJ = calloc(pow(7,N/6),sizeof(int*)); origD = calloc(pow(7,N/6),sizeof(int*)); origQ = calloc(pow(7,N/6),sizeof(int)); + + origK = calloc(pow(7,N/6), sizeof(int)); + + int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together + + + for(j=0; j 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("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(amplitude))); + else + printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+0.00000001)>0?'+':'-' , cabs(cimag(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/strongsim6_relerr.c b/strongsim6_relerr.c new file mode 100644 index 0000000..31bfb1c --- /dev/null +++ b/strongsim6_relerr.c @@ -0,0 +1,432 @@ +#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" + +#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 != 2) { + printf("strongsim6_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + if(N%6 != 0) { + printf("'N' needs to be a multiple of 6 for a k=6 tensor factor decomposition!\n"); + return 1; + } + + 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); + + 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 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))); + + + + 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/strongsim_relerr.c b/strongsim_relerr.c new file mode 100644 index 0000000..24ec5b3 --- /dev/null +++ b/strongsim_relerr.c @@ -0,0 +1,393 @@ +#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" + +#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 != 2) { + printf("strongsim3_rellerr argument: \"number of stabilizer state samples\"\n"); + exit(0); + } + + int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples + + int N; // number of qubits + scanf("%d", &N); + + if(N%3 != 0) { + printf("'N' needs to be a multiple of 3 for a k=3 tensor factor decomposition!\n"); + return 1; + } + + 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); + + 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 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))); + + + + 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/test.bash b/test.bash new file mode 100644 index 0000000..cc380ed --- /dev/null +++ b/test.bash @@ -0,0 +1,30 @@ +#!/bin/bash +# simple Bash script to check if stabilizer rank code works + +# choose the number of qubits and T gates on those qubits +# NOTE: numqubits must be a multiple of your gauss sum tensor multiple! +# e.g. if you test gausssums_multipleof6 then numqubits=6*n for some integer n +numqubits=6 +numTgates=6 + +numruns=20 + +echo "Starting test of $numruns random $numqubits-Pauli expectation values..." +for i in $(seq 1 $numruns) +do +# sleep for 1 second for the pseudorandom number generator +sleep 1;a=$(stdbuf -oL ./randominputPauli $numqubits $numTgates > inputfullPauli.txt && ./strongsim < inputfullPauli.txt | tail -1) +b=$(stdbuf -oL ./multipauli < inputfullPauli.txt | tail -n1) +are=$(echo "$a" | cut -f 1 -d " " | cut -c 1-5); +aim=$(echo "$a" | cut -f 3 -d " " | cut -c 1-5); aimsign=$(echo $a | cut -f 2 -d " "); bimsign=$(echo $a | cut -f 2 -d " "); +bre=$(echo "$b" | cut -f 1 -d " " | cut -c 1-5); +bim=$(echo "$b" | cut -f 3 -d " " | cut -c 1-5); echo "$i: $are $aimsign $aim and $bre $bimsign $bim" +if [ "$are" == "$bre" ] && [ "$aim" == "$bim" ] && [ "$aimsign" == "$bimsign" ] +then +continue +else +echo "NOT EQUAL!" +break +fi +done +echo "Test passed!" diff --git a/test2.bash b/test2.bash new file mode 100644 index 0000000..d5ed868 --- /dev/null +++ b/test2.bash @@ -0,0 +1,42 @@ +#!/bin/bash +# simple Bash diagnostic script to check if non-zero relative error stabilizer rank code works + +# choose the number of qubits and T gates on those qubits +# NOTE: numqubits must be a multiple of your gauss sum tensor multiple! +# e.g. if you test gausssums_multipleof6 then numqubits=6*n for some integer n +numqubits=6 +numTgates=6 +# choose the number of commuting Pauli measurements you want to generate in each run +# (note: can't be greater than $numqubits) +numPaulis=6 + +# number of random stabilizer states to use for stochastic sampling +numsamples=1000 + +# since this is a stochastic calculation, we need a threshold to determine when the calculation is withing close enough relative error +threshold=0.005 # absolute value threshold difference between computed and exact value + +numruns=20 + +echo "Starting test of $numruns random $numqubits-Pauli expectation values calculated by averaging $numsamples random stabilizer states and comparing to threshold $threshold..." + +for i in $(seq 1 $numruns) +do +# sleep for 1 second for the pseudorandom number generator +sleep 1; a=$(stdbuf -oL ./randominputcommutingHermitianPauli2 $numqubits $numTgates $numPaulis > inputPauli.txt && ./strongsim_relerr $numsamples < inputPauli.txt | tail -1) +b=$(stdbuf -oL ./multipauli < inputPauli.txt | tail -n1) +are=$(echo "$a" | cut -f 1 -d " " | cut -c 1-5); +aim=$(echo "$a" | cut -f 3 -d " " | cut -c 1-5); aimsign=$(echo $a | cut -f 2 -d " "); bimsign=$(echo $a | cut -f 2 -d " "); +bre=$(echo "$b" | cut -f 1 -d " " | cut -c 1-5); +bim=$(echo "$b" | cut -f 3 -d " " | cut -c 1-5); echo "$i: $are $aimsign $aim and $bre $bimsign $bim" +rediff=$( printf 'sqrt((%f - %f)^2)\n' "$are" "$bre" | bc -l ) +imdiff=$( printf 'sqrt((%f - %f)^2)\n' "$aim" "$bim" | bc -l ) +if (( $(echo "$rediff < $threshold" |bc -l) )) && (( $(echo "$imdiff < $threshold" |bc -l) )) && [ "$aimsign" == "$bimsign" ] +then + continue + else + echo "NOT EQUAL!" + break +fi +done +echo "Test passed!" diff --git a/timing.bash b/timing.bash new file mode 100644 index 0000000..4ee6c55 --- /dev/null +++ b/timing.bash @@ -0,0 +1,11 @@ +#!/bin/bash +# Bash script to get some timing analysis: +for j in 10 12 14 16 18 +do +for i in {1..1000} +do +./randominputcommutingHermitianPauli $j $j > inputfull.txt && /usr/bin/time -o tmp.txt -a --format=%e ./strongsim < inputfull.txt > /dev/null +sleep 1 +done +echo "" >> tmp.txt +done -- 2.30.2