BSD three-clause license
-Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software.
+Copyright 2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#IDIR =../include
CC=gcc -std=c99
CFLAGS=-Wall
-LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.so libinnerproduct.so
+LIBS=-lm libmatrix.so libexponentialsum.so libextend.so libmeasurepauli.so
-weaksim_relerr: weaksim_relerr.c matrix exponentialsum shrink shrinkstar extend measurepauli innerproduct randomstabilizerstate supplement sparsify
- $(CC) -o $@ weaksim_relerr.c $(CFLAGS) $(LIBS) libshrink.so libshrinkstar.so librandomstabilizerstate.so libsparsify.so
+weaksim_relerr: weaksim_relerr.c matrix exponentialsum shrink shrinkstar extend measurepauli innerproduct randomstabilizerstate supplement supplement2 sparsify
+ $(CC) -o $@ weaksim_relerr.c $(CFLAGS) $(LIBS) libshrink.so libshrinkstar.so librandomstabilizerstate.so libsparsify.so libsupplement.so libsupplement2.so libinnerproduct.so
-weaksim: weaksim.c matrix exponentialsum shrink extend measurepauli innerproduct sparsify
- $(CC) -o $@ weaksim.c $(CFLAGS) $(LIBS) libshrink.so libsparsify.so
+weaksim: weaksim.c matrix exponentialsum shrink extend measurepauli innerproduct supplement supplement2 sparsify
+ $(CC) -o $@ weaksim.c $(CFLAGS) $(LIBS) libshrink.so libsparsify.so libsupplement.so libsupplement2.so libinnerproduct.so -fopenmp
module_sparsify_test: module_sparsify_test matrix sparsify
$(CC) -o $@ module_sparsify_test.c $(CFLAGS) libmatrix.so libsparsify.so
$(CC) -c -Wall -fpic measurepauli.c
$(CC) -shared -o libmeasurepauli.so measurepauli.o -lm libextend.so libshrink.so libmatrix.so
+innerproduct_equatorial: innerproduct_equatorial.h innerproduct_equatorial.c
+ $(CC) -c -Wall -fpic innerproduct_equatorial.c
+ $(CC) -shared -o libinnerproduct_equatorial.so innerproduct_equatorial.o -lm libextend.so libshrink.so libexponentialsum.so libmatrix.so
+
innerproduct: innerproduct.h innerproduct.c
$(CC) -c -Wall -fpic innerproduct.c
$(CC) -shared -o libinnerproduct.so innerproduct.o -lm libextend.so libshrink.so libexponentialsum.so libmatrix.so
$(CC) -c -Wall -fpic shrinkstar.c
$(CC) -shared -o libshrinkstar.so shrinkstar.o -lm libmatrix.so
+randomstabilizerstate_equatorial: randomstabilizerstate_equatorial.h randomstabilizerstate_equatorial.c
+ $(CC) -c -Wall -fpic randomstabilizerstate_equatorial.c
+ $(CC) -shared -o librandomstabilizerstate_equatorial.so randomstabilizerstate_equatorial.o -lm libmatrix.so
+
randomstabilizerstate: randomstabilizerstate.h randomstabilizerstate.c
$(CC) -c -Wall -fpic randomstabilizerstate.c
$(CC) -shared -o librandomstabilizerstate.so randomstabilizerstate.o -lm libmatrix.so libshrinkstar.so -llapacke
$(CC) -c -Wall -fpic supplement.c
$(CC) -shared -o libsupplement.so supplement.o -lm
-sparsify: sparsify.h sparsify.c supplement
+supplement2: supplement2.h supplement2.c
+ $(CC) -c -Wall -fpic supplement2.c
+ $(CC) -shared -o libsupplement2.so supplement2.o -lm
+
+sparsify: sparsify.h sparsify.c supplement2
$(CC) -c -Wall -fpic sparsify.c
- $(CC) -shared -o libsparsify.so sparsify.o -lm libmatrix.so libsupplement.so
+ $(CC) -shared -o libsparsify.so sparsify.o -lm libmatrix.so libsupplement2.so
randominputcommutingHermitianPauli: randominputcommutingHermitianPauli.c
$(CC) -o randominputcommutingHermitianPauli randominputcommutingHermitianPauli.c
--- /dev/null
+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.
--- /dev/null
+Weak simulation of multi-Pauli measurements using approximate correlated stabilizer state decompositions
+-------------------------------------------------------------------------------------------------------
+
+See the draft at [https://arxiv.org/abs/XXXXXX] for an introduction and description of this code. The uncorrelated version is described in [PRX QUANTUM 2, 010345 (2021)] by Seddon, Regula, Pashayan, Ouyang, Campbell (SRPOC).
+
+Note: the correlated sampling is currently only implemented for a maximum of 2t-1 supplemental correlated states!
+
+Note: The stabilizer state representation used in this does not use the CH-form, which ostensibly will speed up simulations here by a factor of approximately 50. We use the implementation described in the Supporting Information of arXiv:1601.07601. A 50-fold improvement is described in arXiv:1808.00128v2.
+
+PREREQUISITES:
+
+You will need the openMP library to run weaksim.c (in particular, the omp.h header).
+For instance, in Debian, you can obtain this by running:
+$ sudo apt install libomp-dev
+
+BUILD:
+
+1. We start with the code that scales quadratically with approximate stabilizer rank.
+$ make weaksim
+
+2. This produces libraries (*.so) of core functions like SHRINK, INNERPRODUCT, EXTEND, etc. since you might want to use them in other code. This means we need to put them somewhere the linker can find them. Follow 3a) for a temporary and easy method that will only be valid for the current shell session and follow 3b) for a persistent method that requires root access.
+
+3a) Update LD_LIBRARY_PATH to see the library:
+$ export LD_LIBRARY_PATH=$PWD:$LD_LIBRARY_PATH
+or
+3b) If you have root privileges put the *so library /usr/local/lib or whatever library directory in your path. Then, use ldconfig to write the path in the config file:
+$ sudo cp ./*so /usr/local/lib
+$ sudo echo "/usr/local/lib" >> /etc/ld.so.conf
+$ sudo ldconfig
+
+4. To run weaksim you need to specify how many parallel processors for openMP to use
+a) if you have no hyperthreading
+export OMP_NUM_THREADS=8
+b) if you have hyperthreading
+export OMP_PROC_BIND=true
+export OMP_NUM_THREADS=8
+
+
+5. Build code for generating random commuting Pauli inputs (necessary for the strongsim[1-12]_relerr.c code since it requires a Hermitian observable to get away with using the square root number of stabilizer states).
+$ make randominputcommutingHermitianPauli2
+
+6. Build Hilbert vector space code to check our results.
+$ make multipauli
+
+7. Run timing analysis for fixed t and delta
+$ bash timing_t_enforceddelta.bash
+Output will be in tmp_$[t]_t.txt for the correlated case and in tmp_$[t]_iid.txt for the i.i.d. case and provide the second-longest and longest runtime, and the largest error (closest to the target delta).
+
+Details of source code files:
+
+weaksim.c takes four arguments corresponding to the target additive error, the angle of the magic state, a "coherent sampling" flag to indicate whether sampling should be correlated and if so how many maximal supplemental states to use (0=no; 1=t; 2=2t-1; 3=xi^3t/2), a seed to use for the i.i.d. random stabilizer state sampling. There is also an optional argument for the "sample number prefactor" which only takes the corresponding fraction of supplemental states indicated by the "coherent sampling" flag, i.e. if you want 0.5 t supplemental states then this argument should be 0.5 and the "coherent sampling" should be 1.
+
+weaksim_relerr.c is a prototype function that also uses random stabilizer state sampling to calculate outcome probabilities from the Haar measure. This is still in a beta stage.
+
+randominputcommutingHermitianPauli.c uses a brute-force random sampling of the commuting Paulis whereas *2.c uses random Clifford operations on initially separable Paulis to generate the commuting Paulis. *2.c is the only practically viable method for more than 30 Paulis on a desktop computer.
+
+NOTE: current strongsim*c only support the same number of T gates as number of qubits!
+
+Pd.txt is a tabulated list of P(d) values for strongsim[1-12]_relerror.c.
void append(struct Node** head_ref, int new_data)
-{
+{
struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
struct Node *last = *head_ref;
--- /dev/null
+#define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
+
+complex double exponentialsum(int *k, int *Q, int *D, int **J);
--- /dev/null
+void extend(int n, int *k, int *h, int **G, int **GBar, int *xi);
--- /dev/null
+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);
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrink.h"
+#include "extend.h"
+#include "exponentialsum.h"
+#include "innerproduct_equatorial.h"
+
+
+/****************************************************************************
+ * Note: Arguments are copied and not modified!
+ ****************************************************************************/
+
+double complex innerproduct_equatorial(int n1, int k1, int *h1, int **G1, int **GBar1, int Q1, int *D1, int **J1, int n2, int k2, int *h2, int **G2, int **GBar2, int Q2, int *D2, int **J2) {
+ // n1 should equal n2 (we assume this)
+
+ int a, b, c;
+
+ int n = n1;
+ int k = k1;
+ int *h = calloc(n, sizeof(int));
+ int **G = calloc(n, sizeof(int*));
+ int **GBar = calloc(n, sizeof(int*));
+ for(a=0; a<n; a++) {
+ h[a] = h1[a];
+ G[a] = calloc(n, sizeof(int));
+ GBar[a] = calloc(n, sizeof(int));
+ for(b=0; b<n; b++) {
+ G[a][b] = G1[a][b];
+ GBar[a][b] = GBar1[a][b];
+ }
+ }
+
+ int Q1copy = Q1;
+ int *D1copy = calloc(k1, sizeof(int));
+ int **J1copy = calloc(k1, sizeof(int*));
+ for(a=0; a<k1; a++) {
+ D1copy[a] = D1[a];
+ J1copy[a] = calloc(k1, sizeof(int));
+ for(b=0; b<k1; b++) {
+ J1copy[a][b] = J1[a][b];
+ }
+ }
+
+ int Q2copy = Q2;
+ int *D2copy = calloc(k2, sizeof(int));
+ int **J2copy = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++) {
+ D2copy[a] = D2[a];
+ J2copy[a] = calloc(k2, sizeof(int));
+ for(b=0; b<k2; b++) {
+ J2copy[a][b] = J2[a][b];
+ }
+ }
+
+ int alpha;
+ int **R = calloc(k, sizeof(int*));
+ for(a=0; a<k; a++)
+ R[a] = calloc(k2, sizeof(int));
+ int **RT = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ RT[a] = calloc(k, sizeof(int));
+ char shrinkResult;
+
+ /* for(a=k2; a<n2; a++) { */
+ /* alpha = dotProductMod(h2, GBar2[a], n2, 2); */
+
+ /* shrinkResult = shrink(n, &k, h, G, GBar, &Q1copy, &D1copy, &J1copy, GBar2[a], alpha); */
+ /* if(shrinkResult == 'E') {// if EMPTY */
+ /* free(h); */
+ /* deallocate_mem(&G, n); deallocate_mem(&GBar, n); */
+ /* free(D1copy); */
+ /* deallocate_mem(&J1copy, k); */
+ /* free(D2copy); */
+ /* deallocate_mem(&J2copy, k2); */
+ /* deallocate_mem(&R, k); */
+ /* deallocate_mem(&RT, k2); */
+ /* return(0.0+0.0*I); */
+ /* } */
+ /* } */
+
+ // Now K = K1 intersect K2 = (n1, k1, h1, G1, GBar1)
+ int *y = calloc(k2, sizeof(int));
+ int *tmpVec = calloc(n2, sizeof(int));
+ addVectorMod(h,h2,tmpVec,n2,2);
+ for(a=0; a<k2; a++) {
+ y[a] = dotProductMod(tmpVec,GBar2[a],n2,2);
+ for(b=0; b<k; b++) {
+ R[b][a] = dotProductMod(G[b],GBar2[a],n2,2);
+ }
+ }
+
+ // IF YOU UNCOMMENT THIS THEN YOU WILL BE CHANGING h2!!!
+ /* for(a=0; a<k2; a++) { */
+ /* for(b=0; b<n2; b++) { */
+ /* h2[b] = (h2[b] + y[a]*G2[a][b])%2; */
+ /* } */
+ /* } */
+ /* printf("h2 (should now match h):\n"); */
+ /* printVector(h2, n2); */
+ /* printf("h:\n"); */
+ /* printVector(h, n2); */
+
+ free(tmpVec);
+
+ // update Q2 and D2 using Eqs. 52 and 53 with y
+ for(a=0; a<k2; a++) {
+ Q2copy = Q2copy + D2copy[a]*y[a];
+ for(b=a+1; b<k2; b++) {
+ Q2copy = Q2copy + J2copy[a][b]*y[a]*y[b];
+ }
+ }
+ Q2copy = Q2copy%8;
+
+ for(a=0; a<k2; a++) {
+ for(b=0; b<k2; b++)
+ D2copy[a] = D2copy[a] + J2copy[a][b]*y[b];
+ D2copy[a] = D2copy[a]%8;
+ }
+
+ free(y);
+
+ // update D2copy and J2copy using Eqs. 49 and 50 with R
+ // (note that R is a kxk2 matrix so this may reduce the dimension of D2copy and J2copy)
+ y = calloc(k, sizeof(int));
+ for(c=0;c<k;c++) {
+ y[c] = 0;
+ for(a=0;a<k2;a++) {
+ y[c] = y[c] + R[c][a]*D2copy[a];
+ for(b=a+1;b<k2;b++)
+ y[c] = y[c] + J2copy[a][b]*R[c][a]*R[c][b];
+ }
+ y[c] = y[c]%8;
+ }
+ free(D2copy);
+ D2copy = calloc(k, sizeof(int));
+ memcpy(D2copy, y, sizeof(int)*k);
+ free(y);
+
+ transp(R,RT,k,k2);
+ int **tmpMatrix = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ tmpMatrix[a] = calloc(k, sizeof(int));
+ multMatrixMod(J2copy,RT,tmpMatrix,k2,k2,k2,k,8);
+ deallocate_mem(&J2copy, k2);
+ J2copy = calloc(k, sizeof(int*));
+ for(a=0;a<k;a++)
+ J2copy[a] = calloc(k, sizeof(int));
+ multMatrixMod(R,tmpMatrix,J2copy,k,k2,k2,k,8);
+
+ deallocate_mem(&tmpMatrix, k2);
+
+ deallocate_mem(&R, k);
+ deallocate_mem(&RT, k2);
+
+ // Now (Q1, D1, J1) and (Q2, D2, J2) are defined in the same basis
+
+ int Q = (Q1copy - Q2copy)%8 < 0 ? (Q1copy - Q2copy)%8 + 8 : (Q1copy - Q2copy)%8;
+ int *D = calloc(k, sizeof(int));
+ int **J = calloc(k, sizeof(int*));
+
+ for(a=0; a<k; a++) { // again, if k=0 then no need to fill these
+ J[a] = calloc(k, sizeof(int));
+ D[a] = (D1copy[a] - D2copy[a])%8 < 0 ? (D1copy[a] - D2copy[a])%8 + 8 : (D1copy[a] - D2copy[a])%8;
+ for(b=0; b<k; b++)
+ J[a][b] = (J1copy[a][b] - J2copy[a][b])%8 < 0 ? (J1copy[a][b] - J2copy[a][b])%8 + 8 : (J1copy[a][b] - J2copy[a][b])%8;
+
+ }
+
+ /* free(h); */
+ /* deallocate_mem(&G, n); deallocate_mem(&GBar, n); */
+
+ double complex amplitude = pow(2,-0.5*(k1+k2))*exponentialsum(&k, &Q, D, J);
+
+ free(D);
+ deallocate_mem(&J, k);
+
+ free(D1copy);
+ deallocate_mem(&J1copy, k);
+
+ free(D2copy);
+ deallocate_mem(&J2copy, k);
+
+ return(amplitude);
+
+}
--- /dev/null
+double complex innerproduct_equatorial(int n1, int k1, int *h1, int **G1, int **GBar1, int Q1, int *D1, int **J1, int n2, int k2, int *h2, int **G2, int **GBar2, int Q2, int *D2, int **J2);
--- /dev/null
+void deallocate_mem(int ***arr, int rows);
+void printMatrix(int** a, int rows, int cols);
+void printVector(int* a, int length);
+int** multMatrix(int **A, int **B, int ro1, int co1, int ro2, int co2);
+void multMatrixMod(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2, int mod);
+int** outerMatrix(int **A, int **B, int ro1, int co1, int ro2, int co2);
+void outerMatrixMod(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2, int mod);
+void outerVectorMod(int *A, int *B, int *C, int ro1, int ro2, int mod);
+void transp(int **a, int **b, int rows, int cols);
+int dotProductMod(int *a, int *b, int length, int mod);
+int trace(int **a, int rows, int cols);
+void scalarmultMatrix(int scalar, int **a, int **b, int rows, int cols);
+int** addMatrix(int **A, int **B, int rows, int cols);
+void addMatrixMod(int **A, int **B, int **C, int rows, int cols, int mod);
+void appendBlockMatrix(int **A, int **B, int **C, int ro1, int co1, int ro2, int co2);
+void addVectorMod(int *A, int *B, int *C, int length, int mod);
+void addSubMatrix(int **A, int **B, int ro1, int co1, int rooff1, int cooff1, int rooff2, int cooff2);
+void appendVector(int *A, int *B, int *C, int ro1, int ro2);
--- /dev/null
+
+double measurepauli(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int m, int *X, int *Z);
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <complex.h>
+#include <math.h>
+
+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<N; i++) {
+ IN = outerMatrix(IN,I2,pow(2,i),pow(2,i),2,2);
+ }
+
+ complex double (*(psiT[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
+ complex double (*(psi0[])) = { (double complex[]) {1.0+0.0*I}, (double complex[]) {0.0*I}}; // T gate magic state
+
+ complex double **psiN;
+ if(K > 0) {
+ psiN = psiT;
+ for(i=1; i<K; i++)
+ psiN = outerMatrix(psiN, psiT, pow(2,i), 1, 2, 1);
+ for(i=K; i<N; i++)
+ psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+ } else {
+ psiN = psi0;
+ for(i=1; i<N; i++)
+ psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+ }
+
+
+ int omega, alpha[N], beta[N], gamma[N], delta[N];
+
+ int Paulicounter = 0;
+
+ complex double **fullP; // full product: \prod_i 1/2*(1+P_i)
+ complex double **P; // P (P_i above) is made up of products of one-qubit Paulis, P1
+ complex double **P1[N]; // one-qubit Paulis
+
+ complex double tr;
+
+ printf("psiN:\n");
+ for(i=0; i<pow(2,N); i++) {
+ printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0]));
+ }
+
+
+ while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) { // go over the product of 1/2*(I+Paulis) that makes up the full projector
+
+ Paulicounter++;
+ if(Paulicounter > N) {
+ printf("Error: Number of Paulis is greater than N!\n");
+ return 1;
+ }
+
+
+ printf("%d\n", omega);
+ for(i=0; i<N; i++) {
+ printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
+ P1[i] = addMatrix(addMatrix(addMatrix(scalarmultMatrix(alpha[i],I2,2,2),scalarmultMatrix(beta[i],Z,2,2),2,2),scalarmultMatrix(gamma[i],X,2,2),2,2),scalarmultMatrix(delta[i],Y,2,2),2,2);
+ }
+
+ P = P1[0];
+ for(i=1; i<N; i++)
+ P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
+ P = scalarmultMatrix(cpow(I,omega),P,pow(2,N),pow(2,N));
+
+ if(Paulicounter == 1)
+ fullP = scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N));
+ else {
+ fullP = multMatrix(scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N)),fullP,pow(2,N),pow(2,N),pow(2,N),pow(2,N));
+ }
+ deallocate_mem(&P, pow(2,N));
+
+ }
+
+ complex double **psiNfinal = multMatrix(fullP,psiN,pow(2,N),pow(2,N),pow(2,N),1);
+
+ printf("fullP:\n");
+ for(i=0; i<pow(2,N); i++) {
+ for(j=0; j<pow(2,N); j++) {
+ printf("%lf+%lfI ", creal(fullP[i][j]), cimag(fullP[i][j]));
+ }
+ printf("\n");
+ }
+ printf("psiNfinal:\n");
+ for(i=0; i<pow(2,N); i++) {
+ printf("%d: %lf+%lfI\n", i, creal(psiNfinal[i][0]), cimag(psiNfinal[i][0]));
+ }
+
+ tr = 0.0 + 0.0*I;
+ //printf("tr:\n");
+ for(i=0; i<pow(2,N); i++) {
+ tr += conj(psiN[i][0])*psiNfinal[i][0];
+ //printf("%d: %lf+%lfI\n", i, creal(tr), cimag(tr));
+ }
+
+ if(creal(tr+0.00000001)>0)
+ 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<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = A[i][j] + B[i][j];
+
+ return C;
+}
+
+complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
+{
+ int i, j;
+
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = scalar*a[i][j];
+
+ return C;
+}
+
+complex double trace(complex double **a, int rows, int cols)
+{
+ int i;
+ complex double tr = 0.0*I;
+
+ for(i=0; i<rows; i++)
+ tr += a[i][i];
+
+ return tr;
+}
+
+complex double** transp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = a[j][i];
+ }
+
+ return C;
+}
+
+complex double** conjtransp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = conj(a[j][i]);
+ }
+
+ return C;
+}
+
+void printMatrix(complex double** a, int rows, int cols)
+{
+ int i, j;
+ printf("Matrix[%d][%d]\n", rows, cols);
+ for(i=0; i<rows; i++) {
+ for(j=0; j<cols; j++) {
+ printf("%lf+%lfI ", creal(a[i][j]), cimag(a[i][j]));
+ }
+ printf("\n");
+ }
+}
+
+complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k;
+ complex double **C;
+ C = calloc(ro1, sizeof(complex double*));
+ for(i=0; i<ro1; i++)
+ C[i] = calloc(co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++) {
+ for(j=0; j<co2; j++) {
+ C[i][j] = 0;
+ for(k=0; k<co1; k++)
+ C[i][j] += A[i][k] * B[k][j];
+ }
+ }
+
+ return C;
+}
+
+complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k, l;
+ complex double **C;
+ C = calloc(ro1*ro2, sizeof(complex double*));
+ for(i=0; i<ro1*ro2; i++)
+ C[i] = calloc(co1*co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ for(k=0; k<co1; k++)
+ for(l=0; l<co2; l++)
+ C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
+
+ return C;
+}
+
+
+
+void deallocate_mem(complex double ***arr, int rows)
+{
+ int i;
+ for(i=0; i<rows; i++)
+ free((*arr)[i]);
+ free(*arr);
+}
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int i;
+
+ if(scanf("%d", omega) != EOF) {
+ for(i=0; i<numqubits; i++) {
+ if(scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]) == EOF) {
+ printf("Error: Too few input coeffs!\n");
+ exit(0);
+ }
+ if(alpha[i]+beta[i]+gamma[i]+delta[i] > 1) {
+ printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+ exit(0);
+ }
+ }
+ return 1;
+ } else
+ return 0;
+
+}
+
+
+
//d = rand()%(n+1);
//printf("d=%d total=%f randPd=%f\n", d, total, randPd);
//printf("!!!d=%d\n", d);
+ //d = 0;
- *k = n; // we will get it to be k=n-d by caling shrinkstar d times below
+ printf("!!\n");
+ fflush(stdout);
+
+ *k = n; // we will get it to be k=n-d by calling shrinkstar d times below
randomX = calloc(n*d,sizeof(float));
randomXcopy = calloc(n*d,sizeof(float));
}
*G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
-
*h = calloc(n, sizeof(int));
for(i=0; i<n; i++) {
--- /dev/null
+void randomstabilizerstate(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J, double** Pd);
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <string.h>
+#include <math.h>
+
+#include "matrix.h"
+
+#include "randomstabilizerstate_equatorial.h"
+
+#define ZEROTHRESHOLD (0.00000001)
+
+// Note: Make sure you seed the random number generator before calling this!
+// i.e. srand((unsigned)time(NULL));
+void randomstabilizerstate_equatorial(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J)
+{
+ // vector and matrix pointers should all be passed unallocated! (no freeing of memory performed)
+
+ int d;
+
+ int i, j;
+
+ d = 0;
+
+ *k = n; // we will get it to be k=n-d by calling shrinkstar d times below
+
+ int info;
+ int numsingvals = -1;
+
+ *G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
+ *h = calloc(n, sizeof(int));
+
+ for(i=0; i<n; i++) {
+ (*G)[i] = calloc(n, sizeof(int));
+ (*GBar)[i] = calloc(n, sizeof(int));
+
+ (*G)[i][i] = 1;
+ (*GBar)[i][i] = 1;
+ }
+
+ // keep the equatorial states with h=0 since it shouldn't matter if they have equal amplitude across all computational basis states
+ for(i=0; i<n; i++)
+ (*h)[i] = 0;//rand()%2;
+
+ *Q = (rand()%8); // Q \in Z/8Z
+
+ *D = calloc(*k, sizeof(int));
+ for(i=0; i<*k; i++)
+ (*D)[i] = (rand()%4)*2; // D_a \in {0, 2, 4, 6}
+
+ *J = calloc(*k, sizeof(int*));
+ for(i=0; i<*k; i++) {
+ (*J)[i] = calloc(*k, sizeof(int));
+ for(j=(i+1); j<*k; j++)
+ (*J)[i][j] = (rand()%2)*4; // J_{a,b} \in {0, 4} for a!=b
+ (*J)[i][i] = (2*(*D)[i])%8;
+ }
+ for(i=0; i<*k; i++) {
+ for(j=(i+1); j<*k; j++) {
+ (*J)[j][i] = (*J)[i][j];
+ }
+ }
+
+ return;
+
+}
--- /dev/null
+void randomstabilizerstate_equatorial(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J);
--- /dev/null
+char shrink(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int *xi, int alpha);
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#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; a<n; a++)
+ G[setWalker->data][a] = (G[setWalker->data][a] + G[i][a])%2;
+
+ // update D and J
+ /*R[setWalker->data][i] = 1;
+ RT[i][setWalker->data] = 1;
+
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;*/
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ /* if(c!=i)
+ tmp += R[c][c]*((*D)[c]);
+ tmp += R[c][i]*((*D)[i]);
+ if(i>c)
+ tmp += ((*J)[c][i])*R[c][c]*R[c][i];
+ else if(i<c)
+ tmp += ((*J)[i][c])*R[c][i]*R[c][c];
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=i) {
+ if(d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ }
+ } else {
+ if(d!=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }
+ }
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp %= 8;
+ tmpMatrix[c][d] = tmp;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));*/
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+
+ if(setWalker->data != i) {
+ R[setWalker->data][i] = 0; // reset R
+ RT[i][setWalker->data] = 0;
+ }
+
+ // update GBar rows, gbar
+ for(a=0; a<n; a++)
+ GBar[i][a] = (GBar[i][a] + GBar[setWalker->data][a])%2;
+
+ setWalker = setWalker->next;
+ }
+
+ // Swap 'i' and 'k'
+ R[*k-1][*k-1] = 0; R[i][i] = 0;
+ R[i][*k-1] = 1; R[*k-1][i] = 1;
+ RT[*k-1][*k-1] = 0; RT[i][i] = 0;
+ RT[i][*k-1] = 1; RT[*k-1][i] = 1;
+
+ tmpVec = G[i];
+ G[i] = G[*k-1];
+ G[*k-1] = tmpVec;
+ tmpVec = GBar[i];
+ GBar[i] = GBar[*k-1];
+ GBar[*k-1] = tmpVec;
+
+ // update D and J
+ //(Note: can't quite directly use previous faster algorithm since R[c][c]!=1 forall c. However, this can still be brought down by a factor of N by the same algorithm appropriately tweaked.)
+ /* if(i!=(*k-1)) {
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;*/
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) { */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* } */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ /* tmp += R[c][c]*((*D)[c]); // iff c!=i and c!=(*k-1) this can be non-zero
+ tmp += R[c][*k-1]*((*D)[*k-1]); // iff c=i this can be non-zero
+ tmp += R[c][i]*((*D)[i]); // iff c=*k-1 this can be non-zero
+ */
+ /* if(i>c) */
+ /* tmp += ((*J)[c][i])*R[c][c]*R[c][i]; */
+ /* else if(i<c) */
+ /* tmp += ((*J)[i][c])*R[c][i]*R[c][c]; */
+ /* if((*k-1)>c) */
+ /* tmp += ((*J)[c][*k-1])*R[c][c]*R[c][*k-1]; */
+ /* else if((*k-1)<c) // pretty sure this can't happen */
+ /* tmp += ((*J)[*k-1][c])*R[c][*k-1]*R[c][c]; */
+ /* if((*k-1)>i) */
+ /* tmp += ((*J)[i][*k-1])*R[c][i]*R[c][*k-1]; */
+ /* else if((*k-1)<i) */
+ /* tmp += ((*J)[*k-1][i])*R[c][*k-1]*R[c][i]; */
+ /* tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=(*k-1) && c!=i) {
+ if(d!=(*k-1) && d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][c]*((*J)[c][*k-1])*R[d][*k-1];
+ }
+ } else {
+ if(d!=(*k-1) && d !=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ tmp += R[c][*k-1]*((*J)[*k-1][d])*R[d][d];
+ } else {
+ tmp += R[c][i]*((*J)[i][*k-1])*R[d][*k-1];
+ tmp += R[c][*k-1]*((*J)[*k-1][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp += R[c][*k-1]*((*J)[*k-1][*k-1])*R[d][*k-1];
+ }
+ }
+ tmpMatrix[c][d] = tmp%8;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));*/
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+ /*}*/
+ // reset R & RT
+ R[*k-1][*k-1] = 1; R[i][i] = 1;
+ R[i][*k-1] = 0; R[*k-1][i] = 0;
+ RT[*k-1][*k-1] = 1; RT[i][i] = 1;
+ RT[i][*k-1] = 0; RT[*k-1][i] = 0;
+
+ for(a=0; a<n; a++)
+ h[a] = (h[a] + beta*G[*k-1][a])%2;
+
+ // update Q & D using Eqs. 52 & 53 with y = beta*G[*k-1][:]
+ // since h is expressed in F_2^n, in F_2^n, y_i = beta*delta_{i,*k-1}
+ // (and so the second sum of Eq. 52 is always over zero terms)
+ /**Q = (*Q + ((*D)[*k-1])*beta)%8;
+
+ for(a=0; a<*k; a++) {
+ (*D)[a] = ((*D)[a] + (*J)[a][*k-1]*beta)%8;
+ }*/
+
+ // remove the kth row & column from J
+ // remove the ith element from D
+
+ // we only need tmpMatrix to be (k-1)x(k-1) and it is kxk
+ /*for(a=0; a<(*k-1); a++) {
+ //tmpMatrix[a] = calloc(*k-1, sizeof(int)); // no need!
+ for(b=0; b<(*k-1); b++)
+ tmpMatrix[a][b] = (*J)[a][b];
+ }
+ deallocate_mem(J, *k);
+ *J = calloc(*k-1, sizeof(int*));
+ for(a=0; a<(*k-1); a++) {
+ (*J)[a] = calloc(*k-1, sizeof(int));
+ for(b=0; b<(*k-1); b++)
+ (*J)[a][b] = tmpMatrix[a][b];
+ }
+ deallocate_mem(&tmpMatrix, *k);*/
+
+ /*tmpVec = calloc(*k-1, sizeof(int));
+ for(a=0; a<(*k-1); a++)
+ tmpVec[a] = (*D)[a];
+ free(*D);
+ *D = calloc(*k-1, sizeof(int));
+ memcpy(*D, tmpVec, sizeof(int)*(*k-1));
+ free(tmpVec);*/
+
+ deallocate_mem(&R, *k);
+ deallocate_mem(&RT, *k);
+
+ // decrement 'k'
+ *k = *k - 1;
+
+ freeList(&setS);
+
+ return 'U'; // SUCCESS
+
+}
+
+
+void append(struct Node** head_ref, int new_data)
+{
+ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
+ struct Node *last = *head_ref;
+
+ new_node->data = new_data;
+ new_node->next = NULL;
+
+ // if the linked list is empty, then make the new node as head
+ if (*head_ref == NULL) {
+ *head_ref = new_node;
+ return;
+ }
+
+ // else traverse till the last node
+ while (last->next != NULL) {
+ last = last->next;
+ }
+ last->next = new_node;
+
+ return;
+}
+
+/*void freeList(struct Node** head_ref)
+{
+
+ struct Node *second_last_walker = *head_ref;
+ struct Node *walker = *head_ref;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ while(walker->next != NULL) {
+ while(walker->next != NULL) {
+ //printf("!%d\n", walker->data);
+ second_last_walker = walker;
+ walker = walker->next;
+ //printf("%d\n", walker->data);
+ //printf("!%d\n", second_last_walker->data);
+ }
+ free(walker);
+ second_last_walker->next = NULL;
+ walker = *head_ref;
+ //printf("!!%d\n", second_last_walker->data);
+ }
+ }
+ free(walker);
+
+ return;
+ }*/
+
+void freeList(struct Node** head_ref)
+{
+
+ struct Node *walker = *head_ref;
+ struct Node *walker2;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ walker2 = walker->next;
+ while(walker2 != NULL) {
+ free(walker);
+ walker = walker2;
+ walker2 = walker->next;
+ }
+ }
+ free(walker);
+
+ return;
+}
+
+void printLinkedList(struct Node *node)
+{
+ if(node == NULL)
+ printf("NULL\n");
+ else {
+ while (node != NULL)
+ {
+ printf(" %d ", node->data);
+ node = node->next;
+ }
+ }
+ printf("\n");
+}
--- /dev/null
+char shrinkstar(int n, int *k, int *h, int **G, int **GBar, int *xi, int alpha);
#include "matrix.h"
#include "supplement.h"
+#include "supplement2.h"
#include "sparsify.h"
/****************************************************************************
// Note: Make sure you seed the random number generator before calling this!
// i.e. srand((unsigned)time(NULL));
-int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling) {
+int sparsify(long **stabStateIndices, int *numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor) {
int i, j, k, l;
double L1Norm = pow(sqrt(1-sin(phi)) + sqrt(1-cos(phi)),T);
//printf("L1Norm squared (= stabilizer extent) = %lf\n", pow(L1Norm,2.0));
+ // fraction of optimal T supplemental states
+ //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/((double)T);
+ //double alpha = 2.0*sqrt(2.0)*delta*pow(L1Norm,2.0)/(4.0+pow(delta,2.0))/((double)T);
+ double alpha;
+
+ // fraction of T supplemental states that we actually need for optimal scaling
+ // if first term is kept from Seddon et al.
+ //double fractionSupplement = alpha*0.25*delta*pow(L1Norm,2.0)/((double)T);
+ // if first term isn't kept
+ double fractionSupplement;
+
printf("delta = %lf\n", delta);
- //printf("L1Norm = %lf\n", L1Norm);
+ printf("L1Norm = %lf\n", L1Norm);
double gamma;
int increment;
- if(coherentSampling) {
+ if(coherentSampling!=0) {
+ if(coherentSampling==2) {
+ alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(2*T-1));
+ fractionSupplement = alpha*1.0;
+ *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(2*T-1))+1)))*(floor(fractionSupplement*(2*T-1))+1);
+ increment = floor(fractionSupplement*(2*T-1))+1;
+ } else if(coherentSampling==1) {
+ alpha = 10.0*pow(delta,2.0)*1.0*pow(L1Norm,2.0)/delta/((double)(T));
+ fractionSupplement = alpha*1.0;
+ *numStabStates = ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta/((double)(floor(fractionSupplement*(T))+1)))*(floor(fractionSupplement*(T))+1);
+ increment = floor(fractionSupplement*(T))+1;
+ }
+
printf("Performing coherent sampling!\n");
- gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
- // since we will be taking a multiple of (T+1) samples, we round numStabStates up to the nearest multiple of (T+1)
- *numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
- increment = T+1;
- } else {
+ printf("Using an optimal fraction of T supplemental states!\n");
+ printf("Optimal fraction = %lf\n", fractionSupplement);
+ // since we will be taking a multiple of (fractionSupplement*T+1) samples, we round numStabStates up to the nearest multiple of (fractionSupplement*T+1)
+ // if first term is kept from Seddon et al.
+ //*numStabStates = ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+ // if first term isn't kept
+ //*numStabStates = ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+ //*numStabStates = ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))/((double)(ceil(fractionSupplement*T)+1)))*(ceil(fractionSupplement*T)+1);
+
+ printf("increment= %d\n", increment);
+ //printf("prebinned numStabStates=%d\n", (int)ceil((8.0/(8.0*delta-sqrt(2.0)*sqrt(16.0-16.0*alpha+pow(alpha,3.0)*pow(delta,4.0))))*pow(L1Norm,2.0)));
+ //printf("prebinned numStabStates=%d\n", (int)ceil(0.0005*(sqrt(2.0)*pow(L1Norm,2.0)/delta-0.5*fractionSupplement*T/pow(delta,2.0)+pow(fractionSupplement*T/L1Norm,2.0)/8.0/sqrt(2.0)/pow(delta,3.0))));
+ //printf("prebinned numStabStates=%d\n", (int)ceil(sqrt(2.0)*(2.0+pow(delta,2.0))*pow(L1Norm,2.0)/delta/(4.0+pow(delta,2.0))));
+ printf("prebinned numStabStates=%d\n", (int)ceil(samplePrefactor*0.05*pow(L1Norm,2.0)/delta));
+ fflush(stdout);
+ if(fractionSupplement > 1.0) {
+ printf("Too many supplemental states requested!\n");
+ printf("Note: the correlated sampling is currently only implemented for a maximum of (2t-1) supplemental correlated states!\n");
+ return 1;
+ }
+ } else if(coherentSampling==0) {
//gamma = 1.0 + (1.0 - 1.0/pow(2.0,0.02*T))*T;
+ // SPARSIFY
//*numStabStates = ceil(((pow(L1Norm,2.0)-gamma)/pow(delta,2.0))/((double)(T+1)))*(T+1);
- gamma = 1.0;
- *numStabStates = ceil((pow(L1Norm,2.0)-gamma)/pow(delta,2.0));
+ // if first term is kept from Seddon et al.
+ //*numStabStates = ceil((2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
+ // if first term isn't kept
+ *numStabStates = ceil(samplePrefactor*(1.0*2.0+sqrt(2.0))*(pow(L1Norm,2.0))/pow(delta,1.0));
+ // multiply the factor of 2.0 by 1.0 to get the prior state-of-the-art
+ // multiply the factor of 2.0 by 0.0 to get the tighter bound on the prior state-of-the-art
increment = 1;
+ } else {
+ printf("Error: Unknown value for \"coherent sampling\" flag");
+ return 1;
}
}
(*stabStateIndices)[i] = index;
//printf("first index = %d\n", index);
- if(coherentSampling)
- for(k=1; k<=T; k++)
+ if(coherentSampling==2)
+ //for(k=1; k<=T; k++)
+ //(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
+ supplement2(&(*stabStateIndices)[i], T, fractionSupplement);
+ else if(coherentSampling==1)
+ //for(k=1; k<=T; k++)
//(*stabStateIndices)[i+k] = index ^ (long)pow(2,k-1); // flip the (k-1)th bit index
- supplement(&(*stabStateIndices)[i], T);
+ supplement(&(*stabStateIndices)[i], T, fractionSupplement);
break;
}
//printf("%lf\n", sum);
--- /dev/null
+#define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
+
+double binomcoeff(int m, int n);
+int sparsify(long** stabStateIndices, int* numStabStates, int T, double phi, double delta, int coherentSampling, double samplePrefactor);
}
-void supplement(long stabStateIndices[], int T)
+void supplement(long stabStateIndices[], int T, double fractionSupplement)
{
int K = round(log(T)/log(2)); // T = 2^K
unsigned long bitstring;
int bitstringCounter = 1; // bitstringCounter starts at 1 since 0 is the given bitstring
+ int bitstringCounterMax = floor(fractionSupplement*T);
int i;
+ if(bitstringCounter > bitstringCounterMax) return;
+
// first added bitstring is alpha ... alpha for alpha = 01
bitstring = 1;
for(i=1; i<T/2; i++) {
//printf("%lu\n", bitstring);
stabStateIndices[bitstringCounter] = given^(bitstring);
bitstringCounter++;
+ if(bitstringCounter > bitstringCounterMax) return;
// second added bitstring is beta ... beta for beta = 10
bitstring = 2;
//printf("%lu\n", bitstring);
stabStateIndices[bitstringCounter] = given^bitstring;
bitstringCounter++;
+ if(bitstringCounter > bitstringCounterMax) return;
maskList[0] = round(pow(2,pow(2,K)))-1;
//printf("new bitstring=");
//printBinary(stabStateIndices[bitstringCounter], T);
bitstringCounter++;
+ if(bitstringCounter > bitstringCounterMax) return;
}
maskList[treeDepth] = levelMask;
stabStateIndices[0] = (int)(pow(2,T))-1;
- supplement(&stabStateIndices[0], T);
+ supplement(&stabStateIndices[0], T, 1.0);
int i;
for(i=0; i<T+1; i++)
--- /dev/null
+void supplement(long stabStateIndices[], int T, double fractionSupplement);
+void printBinary(unsigned long number, int T);
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+void printBinary(unsigned long number, int T) {
+
+ int i;
+
+ for(i=0; i<T; i++)
+ printf("%d", (int)(number/pow(2,i))%2);
+ printf("\n");
+
+}
+
+void supplement2(long stabStateIndices[], int T, double fractionSupplement)
+{
+
+ int bitstringCounterMax = floor(fractionSupplement*(2*T-1));
+ if(bitstringCounterMax < 1)
+ return;
+
+ int K = round(log(T)/log(2)); // T = 2^K
+
+ //printf("K=%d\n", K);
+
+ unsigned long given = (unsigned long)stabStateIndices[0]; // the given bitstring
+
+ //printf("given bitstring =");
+ //printBinary(stabStateIndices[0], T);
+
+ unsigned long long precurserBitstrings[2*T-1];
+ unsigned long long tmpBitstrings[2*T-1];
+
+ int precurserBitstringCounter = 0;
+
+ int i, j, k;
+
+ // first added bitstring is alpha = 10
+ precurserBitstrings[0] = 2;
+ //printf("added precurser=");
+ //printBinary(precurserBitstrings[0], T);
+ // second added bitstring is beta = 01
+ precurserBitstrings[1] = 1;
+ //printf("added precurser=");
+ //printBinary(precurserBitstrings[1], T);
+ // third added bitstring is gamma = 11
+ precurserBitstrings[2] = 0;
+ //printf("added precurser=");
+ //printBinary(precurserBitstrings[2], T);
+ precurserBitstringCounter+=3;
+
+ long mask;
+ for(i=1; i<=K-1; i++) {
+ mask = 0;
+ for(j=0;j<pow(2,i-1);j++)
+ mask += 3*pow(2,2*j+pow(2,i));
+ //printf("mask=");
+ //printBinary(mask, T);
+ for(j=0; j<precurserBitstringCounter; j++) {
+ tmpBitstrings[j] = precurserBitstrings[j]+pow(2,pow(2,i))*precurserBitstrings[j];
+ //printf("added precurser=");
+ //printBinary(tmpBitstrings[j], T);
+ tmpBitstrings[j+precurserBitstringCounter] = precurserBitstrings[j]+((long)(pow(2,pow(2,i))*(precurserBitstrings[j]))^mask);
+ //printf("!!added precurser=");
+ //printBinary(tmpBitstrings[j+precurserBitstringCounter], T);
+ }
+ tmpBitstrings[2*precurserBitstringCounter] = (long)((tmpBitstrings[2])^((long)(mask/pow(2,pow(2,i)))));
+ //printf("!added precurser=");
+ //printBinary(tmpBitstrings[2*precurserBitstringCounter], T);
+ precurserBitstringCounter = pow(2,i+2)-1;
+ memcpy(precurserBitstrings,tmpBitstrings,precurserBitstringCounter*sizeof(unsigned long));
+ }
+
+ int bitstringCounter = 1; // bitstringCounter starts at 1 since 0 is the given bitstring
+
+ mask = mask+mask/pow(2,pow(2,K-1));
+ for(i=0; i<2*T-1; i++) {
+ stabStateIndices[bitstringCounter] = given^(precurserBitstrings[i]^((long)mask));
+ //stabStateIndices[bitstringCounter] = given^(precurserBitstrings[i]);
+ //printf("new bitstring=");
+ //printBinary(stabStateIndices[bitstringCounter], T);
+ bitstringCounter++;
+ if(bitstringCounter > bitstringCounterMax) return;
+ }
+
+}
+
+
+int main(int argc, char *argv[])
+{
+
+ int T = (int)pow(2,4);
+ long* stabStateIndices;
+
+ stabStateIndices = calloc(2*T+1,sizeof(long));
+
+ stabStateIndices[0] = (int)(pow(2,T))-1;
+
+ printBinary(stabStateIndices[0],T);
+ printf("!!\n");
+
+ supplement2(&stabStateIndices[0], T, 1.0);
+
+ int i;
+ for(i=0; i<2*T; i++)
+ printBinary(stabStateIndices[i], T);
+ printf("\n");
+
+ free(stabStateIndices);
+
+ return 0;
+
+}
--- /dev/null
+void supplement2(long stabStateIndices[], int T, double fractionSupplement);
+void printBinary(unsigned long number, int T);
--- /dev/null
+#!/bin/bash
+# Bash script to get some timing analysis
+
+# Runtime for different t's
+
+# number of runs to run in parallel
+maxthreads=10
+# NOTE: that this is not the same as running weaksim in parallel! For that you need to set OMP_NUM_THREADS appropriately (see README.txt)
+# Note that setting OMP_NUM_THREADS>1 will lead to more pessimistic runtime analysis, due to parallelization overhead, though the runs should complete faster!
+
+## For the most accurate time analysis, run with
+## export OMP_NUM_THREADS=1
+## export OMP_PROC_BIND=false
+#export OMP_NUM_THREADS=6
+#export OMP_PROC_BIND=true
+export OMP_NUM_THREADS=1
+export OMP_PROC_BIND=false
+
+delta=0.1
+deltasquared=`bc -l <<< $delta^2`
+
+numruns=10 #1000
+
+#for t in 4 8
+for t in 8
+#for t in 16
+#for t in 32
+#for t in 2 4 8 #16 #32
+#for t in 8 16 32
+do
+ #delta=`bc -l <<< 10^-$deltapower`
+ echo "t=$t"
+ # alpha is the fraction of t supplemental states that we will take to keep the cost the same
+ # alpha = delta*xi^t/t
+ #echo "0.25*$delta*1.17^$t/$t"
+ #alpha=`bc -l <<< 0.25*$delta*1.17^$t/$t`
+ #%echo "alpha=$alpha"
+ TIMEFORMAT='%3R';
+ #number of runs
+ for (( i=1; i<=$numruns; i++ ))
+ do
+ # Generate two t-Pauli measurements
+ ./randominputcommutingHermitianPauli2 $t $t 1000 | head -n $[2*$t+4] > inputPauli_$[i].txt
+ sleep 1
+ done
+ for (( i=1; i<=$numruns; i++ ))
+ do
+ # Exact measurements
+ ./multipauli < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 1 > tmp_$[t]_exact_$[i].txt &
+ #seed=`od -A n -t d -N 4 /dev/urandom |tr -d ' '|tr -d '-'`
+ #./weaksim 0.001 0.25 2 $seed < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 3 > tmp_$[t]_exact_$[i].txt &
+ #./weaksim 0.01 0.25 2 $seed < inputPauli_$[i].txt | tail -1 | cut -d ' ' -f 3 > tmp_$[t]_exact_$[i].txt &
+ # minus one below so as not to count grep search job
+ numthreads=$[`ps -aux | grep multipauli | wc -l`-1]
+ #numthreads=$[`ps -aux | grep weaksim | wc -l`-1]
+ while [ $numthreads -gt $[$maxthreads-1] ]
+ do
+ sleep 10
+ numthreads=$[`ps -aux | grep multipauli | wc -l`-1]
+ #numthreads=$[`ps -aux | grep weaksim | wc -l`-1]
+ done
+ done
+
+ for (( i=1; i<=$numruns; i++ ))
+ do
+ echo "run: $i"
+
+ approxerror=1.0
+ approxcoherror=1.0
+
+ samplePrefactor=0.001
+ seed=`od -A n -t d -N 4 /dev/urandom |tr -d ' '|tr -d '-'`
+ while (( $(echo "$approxerror > $deltasquared" |bc -l) )); do
+ echo "samplePrefactor=$samplePrefactor"
+ ( time ( ./weaksim $delta 0.25 0 $seed $samplePrefactor < inputPauli_$[i].txt >> tmp_$[t]_iid_$[i].txt ) ) 2> tmp_$[t]_iid_$[i]_timing.txt
+
+ exactanswer=`tail -1 tmp_$[t]_exact_$[i].txt`
+ approxanswer=`tail -1 tmp_$[t]_iid_$[i].txt | cut -d ' ' -f 3`
+ approxerror=`bc -l <<< "sqrt(($exactanswer-$approxanswer)^2)"`
+ echo $approxerror
+ samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor*1.1)/1"`
+ done
+ samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor/1.1)/1"`
+ echo "samplePrefactor used=$samplePrefactor"
+ echo $approxerror >> tmp_$[t]_iid_$[i]_error.txt
+
+ samplePrefactor=0.001
+ while (( $(echo "$approxcoherror > $deltasquared" |bc -l) )); do
+ echo "samplePrefactor=$samplePrefactor"
+ ( time ( ./weaksim $delta 0.25 2 $seed $samplePrefactor < inputPauli_$[i].txt >> tmp_$[t]_t_$[i].txt ) ) 2> tmp_$[t]_t_$[i]_timing.txt
+
+ exactanswer=`tail -1 tmp_$[t]_exact_$[i].txt`
+ approxcoherentanswer=`tail -1 tmp_$[t]_t_$[i].txt | cut -d ' ' -f 3`
+ approxcoherror=`bc -l <<< "sqrt(($exactanswer-$approxcoherentanswer)^2)"`
+ echo $approxcoherror
+ samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor*1.1)/1"`
+ done
+ samplePrefactor=`bc -l <<< "scale=4;($samplePrefactor/1.1)/1"`
+ echo "samplePrefactor used=$samplePrefactor"
+ echo $approxcoherror >> tmp_$[t]_t_$[i]_error.txt
+
+ done
+
+ for (( i=1; i<=$numruns; i++ ))
+ do
+ rm tmp_$[t]_exact_$[i].txt;
+ rm tmp_$[t]_iid_$[i].txt;
+ rm tmp_$[t]_t_$[i].txt;
+ rm inputPauli_$[i].txt;
+ #echo "*******************" >> tmp_$[t]_iid.txt
+ #echo "*******************" >> tmp_$[t]_t.txt
+ # sleep 1
+ done
+ echo "Longest runtime:" >> tmp_$[t]_iid.txt
+ cat tmp_$[t]_iid_[1-9]*_timing.txt | sort -n | tail -2 >> tmp_$[t]_iid.txt
+ rm tmp_$[t]_iid_[1-9]*_timing.txt
+ echo "Longest runtime:" >> tmp_$[t]_t.txt
+ cat tmp_$[t]_t_[1-9]*_timing.txt | sort -n | tail -2 >> tmp_$[t]_t.txt
+ rm tmp_$[t]_t_[1-9]*_timing.txt
+ echo "Largest error:" >> tmp_$[t]_iid.txt
+ cat tmp_$[t]_iid_[1-9]*_error.txt | sort -n | tail -1 >> tmp_$[t]_iid.txt
+ rm tmp_$[t]_iid_[1-9]*_error.txt
+ echo "Largest error:" >> tmp_$[t]_t.txt
+ cat tmp_$[t]_t_[1-9]*_error.txt | sort -n | tail -1 >> tmp_$[t]_t.txt
+ rm tmp_$[t]_t_[1-9]*_error.txt
+ echo "delta=$delta done!"
+done
projects / weak_simulation_stab_extent.git / commitdiff