diff --git a/bench/stim-compare/.gitignore b/bench/stim-compare/.gitignore new file mode 100644 index 000000000..d57c36f39 --- /dev/null +++ b/bench/stim-compare/.gitignore @@ -0,0 +1,5 @@ +target*/ +prof/ +__pycache__/ +# Regenerated by sweep.py. +circuits/ diff --git a/bench/stim-compare/Cargo.lock b/bench/stim-compare/Cargo.lock new file mode 100644 index 000000000..9830648e9 --- /dev/null +++ b/bench/stim-compare/Cargo.lock @@ -0,0 +1,1270 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "ahash" +version = "0.8.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a15f179cd60c4584b8a8c596927aadc462e27f2ca70c04e0071964a73ba7a75" +dependencies = [ + "cfg-if", + "getrandom 0.3.4", + "once_cell", + "version_check", + "zerocopy", +] + +[[package]] +name = "aho-corasick" +version = "1.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c982642fa9e8606056828ee9a8505737230110bb1099153c79efe865c59d12ba" +dependencies = [ + "memchr", +] + +[[package]] +name = "allocator-api2" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923" + +[[package]] +name = "anstyle" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "940b3a0ca603d1eade50a4846a2afffd5ef57a9feac2c0e2ec2e14f9ead76000" + +[[package]] +name = "anyhow" +version = "1.0.104" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "330a5ed07fa54e4702c9d6c4174f74427fc0ef6e214bbd677ae50a5099946470" + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "ar_archive_writer" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73cd58deff2140a0a8eae87e417bd01db68a33e148aa93d1e8cd837e55e312b6" +dependencies = [ + "object", +] + +[[package]] +name = "ariadne" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "36f5e3dca4e09a6f340a61a0e9c7b61e030c69fc27bf29d73218f7e5e3b7638f" +dependencies = [ + "unicode-width 0.1.14", + "yansi", +] + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bitflags" +version = "2.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06" + +[[package]] +name = "bitvec" +version = "1.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ddcec3d12c579d40898fe0a9a358a803c23e9c52ca3c425707f81c9436211837" +dependencies = [ + "funty", + "radium", + "tap", + "wyz", +] + +[[package]] +name = "bnum" +version = "0.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "119771309b95163ec7aaf79810da82f7cd0599c19722d48b9c03894dca833966" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "bon" +version = "3.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60eafe0d77c3a2fc292c1d1346c3041b33c0a108085a2afabf672b70f69dbbc9" +dependencies = [ + "bon-macros", +] + +[[package]] +name = "bon-macros" +version = "3.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd0f9631d8aaaee112c41985d675ef269e02acbd4f33122836af4f0c5f699ff6" +dependencies = [ + "darling", + "ident_case", + "prettyplease", + "proc-macro2", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "bumpalo" +version = "3.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cc" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f360145194ee8e21db5ee7f3fcd4fe52210864c75c985dae33218202c8bbe040" +dependencies = [ + "find-msvc-tools", + "shlex", +] + +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "chacha20" +version = "0.10.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65c35e4b699c7e15ccbe7ee35c005e4fc0a278d22238a2857e6ce2dadeda1b06" +dependencies = [ + "cfg-if", + "cpufeatures", + "rand_core", +] + +[[package]] +name = "chumsky" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ba4a05c9ce83b07de31b31c874e87c069881ac4355db9e752e3a55c11ec75a6" +dependencies = [ + "hashbrown 0.15.5", + "regex-automata", + "serde", + "stacker", + "unicode-ident", + "unicode-segmentation", +] + +[[package]] +name = "clap" +version = "4.6.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa8876b300ab35ba921adea3dfd70157a46249b33f95c9084ae5709785478946" +dependencies = [ + "clap_builder", +] + +[[package]] +name = "clap_builder" +version = "4.6.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0797fb7aeb1406c84efac526901f7ec3ead2124f946b494e72879d4b54704d" +dependencies = [ + "anstyle", + "clap_lex", + "strsim", +] + +[[package]] +name = "clap_lex" +version = "1.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c133bc6a41be0d194c306b5506d15e6feeea7b1d6604bd3f8310dfb2ca96486" + +[[package]] +name = "codespan-reporting" +version = "0.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "af491d569909a7e4dee0ad7db7f5341fef5c614d5b8ec8cf765732aba3cff681" +dependencies = [ + "serde", + "termcolor", + "unicode-width 0.2.2", +] + +[[package]] +name = "cpufeatures" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5ca28b0ae3115b884660db4118d803791fd6756b6e88f39c0f3f7859060d7566" +dependencies = [ + "libc", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "622f3fc73690be383c7214310406f28a90e6edeadc3cea882f9d71e495b9711a" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc74980687109a3b14c72fd458107bf0baa1da1a1a805e178d15501ba9b86d9d" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a31eee39dddec8330830986fcd7625edb5a24ec90ea038215273bbc3adb08ac6" + +[[package]] +name = "cxx" +version = "1.0.202" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13f6de320895f42e6e081abb5c7983bedcf0b6d0ff9323de0d33f620c8ac1199" +dependencies = [ + "cc", + "cxx-build", + "cxxbridge-cmd", + "cxxbridge-flags", + "cxxbridge-macro", + "foldhash 0.2.0", + "link-cplusplus", +] + +[[package]] +name = "cxx-build" +version = "1.0.202" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fde53ca86b9704a943fef0f1e1d836239a6aedca5de3c65fa9f97ec0bd46d39" +dependencies = [ + "cc", + "codespan-reporting", + "indexmap", + "proc-macro2", + "quote", + "scratch", + "syn 3.0.6", +] + +[[package]] +name = "cxxbridge-cmd" +version = "1.0.202" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "07bae89236c811fd4d08ed3441759ac4ac98d752cbfd3de341315ba16ad20ec3" +dependencies = [ + "clap", + "codespan-reporting", + "indexmap", + "proc-macro2", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "cxxbridge-flags" +version = "1.0.202" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "49045042e5fced01b80742aba5508de82aa4f13677ed3fa2b4cda709c40c5918" + +[[package]] +name = "cxxbridge-macro" +version = "1.0.202" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b181252e2e3d3b5d183afbdc033a3958b445eb3e1a0ae65fc3d4f5259f5da6fd" +dependencies = [ + "indexmap", + "proc-macro2", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "darling" +version = "0.24.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed17f5901b6630b993ca003def43f2f8ef4014fc13b047b57aad617ff32bc2ec" +dependencies = [ + "darling_core", + "darling_macro", +] + +[[package]] +name = "darling_core" +version = "0.24.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6837e2cf7485aaae18f86181d2f0e9a7ed297a025e220aeabf63fdebd3a2ddff" +dependencies = [ + "ident_case", + "proc-macro2", + "quote", + "strsim", + "syn 3.0.6", +] + +[[package]] +name = "darling_macro" +version = "0.24.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2ac7135c3ef02b2f7833bbeb1be5ba7f966dcde8a87c6b87f65a778d71a02785" +dependencies = [ + "darling_core", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "dashmap" +version = "6.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6361d5c062261c78a176addb82d4c821ae42bed6089de0e12603cd25de2059c" +dependencies = [ + "cfg-if", + "crossbeam-utils", + "hashbrown 0.14.5", + "lock_api", + "once_cell", + "parking_lot_core", + "rayon", +] + +[[package]] +name = "either" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "252afb9ae5eaa683babdc6a068b3f5726eb19e05070c731f9b2a23a7c3e8ed34" + +[[package]] +name = "equivalent" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f" + +[[package]] +name = "find-msvc-tools" +version = "0.1.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aedcfb3409746eddb02b9e19ebda1c3394f759a152e48ee875a0844d1b955484" + +[[package]] +name = "foldhash" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2" + +[[package]] +name = "foldhash" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "77ce24cb58228fbb8aa041425bb1050850ac19177686ea6e0f41a70416f56fdb" + +[[package]] +name = "funty" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6d5a32815ae3f33302d95fdcb2ce17862f8c65363dcfd29360480ba1001fc9c" + +[[package]] +name = "fxhash" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c31b6d751ae2c7f11320402d34e41349dd1016f8d5d45e48c4312bc8625af50c" +dependencies = [ + "byteorder", +] + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi 5.3.0", + "wasip2", +] + +[[package]] +name = "getrandom" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "300e883d756b2e4ec94e02791f39b04b522276138852cfc41d9fb7e904106099" +dependencies = [ + "cfg-if", + "js-sys", + "libc", + "r-efi 6.0.0", + "rand_core", + "wasm-bindgen", +] + +[[package]] +name = "gxhash" +version = "3.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f3ce1bab7aa741d4e7042b2aae415b78741f267a98a7271ea226cd5ba6c43d7d" +dependencies = [ + "rustversion", +] + +[[package]] +name = "hashbrown" +version = "0.14.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" + +[[package]] +name = "hashbrown" +version = "0.15.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9229cfe53dfd69f0609a49f65461bd93001ea1ef889cd5529dd176593f5338a1" +dependencies = [ + "allocator-api2", + "equivalent", + "foldhash 0.1.5", +] + +[[package]] +name = "hashbrown" +version = "0.17.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed5909b6e89a2db4456e54cd5f673791d7eca6732202bbf2a9cc504fe2f9b84a" +dependencies = [ + "allocator-api2", + "equivalent", + "foldhash 0.2.0", +] + +[[package]] +name = "ident_case" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" + +[[package]] +name = "indexmap" +version = "2.14.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cc4e190f5d26ca7051642629da2c52fc03bde85a03197c99408dcd291734c855" +dependencies = [ + "equivalent", + "hashbrown 0.17.1", +] + +[[package]] +name = "itertools" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b192c782037fadd9cfa75548310488aabdbf3d2da73885b31bd0abd03351285" +dependencies = [ + "either", +] + +[[package]] +name = "js-sys" +version = "0.3.106" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7883d941dae510fb2d978fc3fe018c71c9e2892fd38854de3e8b92c2e5ad9cc5" +dependencies = [ + "cfg-if", + "wasm-bindgen", +] + +[[package]] +name = "libc" +version = "0.2.189" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2" + +[[package]] +name = "link-cplusplus" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f78c730aaa7d0b9336a299029ea49f9ee53b0ed06e9202e8cb7db9bae7b8c82" +dependencies = [ + "cc", +] + +[[package]] +name = "lock_api" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "224399e74b87b5f3557511d98dff8b14089b3dadafcab6bb93eab67d3aace965" +dependencies = [ + "scopeguard", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f607c237553f086e7043417a51df26b2eb899d3caff94e6a67592ff992fedc7" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "memchr" +version = "2.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98" + +[[package]] +name = "ndarray" +version = "0.17.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "520080814a7a6b4a6e9070823bb24b4531daac8c4627e08ba5de8c5ef2f2752d" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "portable-atomic", + "portable-atomic-util", + "rawpointer", +] + +[[package]] +name = "num" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "35bd024e8b2ff75562e5f34e7f4905839deb4b22955ef5e73d2fea1b9813cb23" +dependencies = [ + "num-bigint", + "num-complex", + "num-integer", + "num-iter", + "num-rational", + "num-traits", +] + +[[package]] +name = "num-bigint" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-iter" +version = "0.1.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c92800bd69a1eac91786bcfe9da64a897eb72911b8dc3095decbd07429e8048b" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "object" +version = "0.39.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2e5a6c098c7a3b6547378093f5cc30bc54fd361ce711e05293a5cc589562739b" +dependencies = [ + "memchr", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "parking_lot_core" +version = "0.9.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2621685985a2ebf1c516881c026032ac7deafcda1a2c9b7850dc81e3dfcb64c1" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-link", +] + +[[package]] +name = "portable-atomic" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05c8b63e8d9609db387f0324918f81d68fe27748f084ef092fb35954d0539a85" + +[[package]] +name = "portable-atomic-util" +version = "0.2.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "10ab3eb7f3becc3a1cbc4f2c6f20267996cfc1a6467a873763411b136a122715" +dependencies = [ + "portable-atomic", +] + +[[package]] +name = "ppvm-lossy-pauli-word-2" +version = "0.1.0" +dependencies = [ + "bitvec", + "bytemuck", + "fxhash", + "ppvm-pauli-word-2", + "ppvm-traits-2", +] + +[[package]] +name = "ppvm-pauli-sum" +version = "0.1.0" +dependencies = [ + "approx", + "bon", + "dashmap", + "fxhash", + "gxhash", + "indexmap", + "itertools", + "num", + "ppvm-pauli-word", + "ppvm-traits", +] + +[[package]] +name = "ppvm-pauli-sum-2" +version = "0.1.0" +dependencies = [ + "approx", + "hashbrown 0.17.1", + "indexmap", + "num", + "ppvm-lossy-pauli-word-2", + "ppvm-pauli-word-2", + "ppvm-traits-2", + "rand", +] + +[[package]] +name = "ppvm-pauli-word" +version = "0.1.0" +dependencies = [ + "anyhow", + "bitvec", + "bytemuck", + "fxhash", + "gxhash", + "itertools", + "num", + "ppvm-traits", +] + +[[package]] +name = "ppvm-pauli-word-2" +version = "0.1.0" +dependencies = [ + "bitvec", + "bytemuck", + "fxhash", + "num", + "ppvm-traits-2", +] + +[[package]] +name = "ppvm-stim" +version = "0.1.0" +dependencies = [ + "bitvec", + "itertools", + "num", + "ppvm-pauli-sum", + "ppvm-tableau", + "ppvm-tableau-2", + "rand", + "smallvec", + "stim-parser", + "thiserror", +] + +[[package]] +name = "ppvm-tableau" +version = "0.1.0" +dependencies = [ + "bitvec", + "bnum", + "fxhash", + "getrandom 0.4.3", + "itertools", + "num", + "ppvm-pauli-sum", + "ppvm-pauli-word", + "ppvm-traits", + "rand", + "smallvec", +] + +[[package]] +name = "ppvm-tableau-2" +version = "0.1.0" +dependencies = [ + "bnum", + "bytemuck", + "fxhash", + "getrandom 0.4.3", + "gxhash", + "num", + "ppvm-pauli-sum-2", + "ppvm-pauli-word-2", + "ppvm-traits-2", + "rand", +] + +[[package]] +name = "ppvm-traits" +version = "0.1.0" +dependencies = [ + "ahash", + "bitvec", + "bytemuck", + "dashmap", + "fxhash", + "gxhash", + "indexmap", + "num", + "rayon", +] + +[[package]] +name = "ppvm-traits-2" +version = "0.1.0" +dependencies = [ + "num", + "rand", +] + +[[package]] +name = "prettyplease" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2bfe0f4c752e450fc2faf62654f1c134747922825d5b04ca717b8874f41a40c0" +dependencies = [ + "proc-macro2", + "syn 3.0.6", +] + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "psm" +version = "0.1.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4dcd034599e63b970727f70d79e02d62390a4a84f7c6b827c27c46d5ac3fa622" +dependencies = [ + "ar_archive_writer", + "cc", +] + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + +[[package]] +name = "r-efi" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8dcc9c7d52a811697d2151c701e0d08956f92b0e24136cf4cf27b57a6a0d9bf" + +[[package]] +name = "radium" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc33ff2d4973d518d823d61aa239014831e521c75da58e3df4840d3f47749d09" + +[[package]] +name = "rand" +version = "0.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "65c9fb96cbc91e3478eaae79a69fcd3f1ae4ad052e471fe6732fff548984b4af" +dependencies = [ + "chacha20", + "getrandom 0.4.3", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "63b8176103e19a2643978565ca18b50549f6101881c443590420e4dc998a3c69" + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "rayon" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "redox_syscall" +version = "0.5.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed2bf2547551a7053d6fdfafda3f938979645c44812fbfcda098faae3f1a362d" +dependencies = [ + "bitflags", +] + +[[package]] +name = "regex-automata" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "59b23e92ee4318893fa3fe3e6fb365258efbfe6ac6ab30f090cdcbb7aa37efa9" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.7.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dbb5fb1acd8a1a18b3dd5be62d25485eb770e05afb408a9627d14d451bae12da" + +[[package]] +name = "rustversion" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf54715a573b99ac80df0bc206da022bcd442c974952c7b9720069370852e21f" + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "scratch" +version = "1.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d68f2ec51b097e4c1a75b681a8bec621909b5e91f15bb7b840c4f2f7b01148b2" + +[[package]] +name = "serde" +version = "1.0.229" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4148590afebada386688f18773da617792bf2ef03ffc1e4cbd2b1d45b023e0ba" +dependencies = [ + "serde_core", + "serde_derive", +] + +[[package]] +name = "serde_core" +version = "1.0.229" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "67dca2c9c51e58a4791a4b1ed58308b39c64224d349a935ab5039aa360942a48" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.229" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e7a5d71263a5a7d47b41f6b3f06ba276f10cc18b0931f1799f710578e2309348" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "shlex" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8fadd59c855ef2080decdef8ff161eb6661b86933c9d82e5ba29dc602a55aba" + +[[package]] +name = "smallvec" +version = "1.16.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f9395f0f0eee849a9b707b2f06bb92a6a422090e2123bb2ef8e87a0e61892a8e" + +[[package]] +name = "stacker" +version = "0.1.25" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "707f49d46706bacf8a2b00d51dace3f9de527c13eec3778f570c411f89e69967" +dependencies = [ + "cc", + "cfg-if", + "libc", + "psm", + "windows-sys", +] + +[[package]] +name = "stim" +version = "0.4.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c2bb10313d1124650826574b6a4ba38b4e0e2a069ceaf396d71f84451bb8885" +dependencies = [ + "cxx", + "ndarray", + "num-complex", + "stim-cxx", +] + +[[package]] +name = "stim-compare" +version = "0.0.0" +dependencies = [ + "bnum", + "ppvm-stim", + "rand", + "stim", +] + +[[package]] +name = "stim-cxx" +version = "0.4.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "957eb60a924d9e6060deb5952584adfa51ac5796225e6059e16176a4ae0967e7" +dependencies = [ + "cxx", + "cxx-build", +] + +[[package]] +name = "stim-parser" +version = "0.1.0" +dependencies = [ + "ariadne", + "chumsky", +] + +[[package]] +name = "strsim" +version = "0.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f" + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8593e8e72159ed2257d083c7a454a85cbf854f37a0966d8d483aff8c8a3ebcee" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "tap" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "55937e1799185b12863d447f42597ed69d9928686b8d88a1df17376a097d8369" + +[[package]] +name = "termcolor" +version = "1.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "06794f8f6c5c898b3275aebefa6b8a1cb24cd2c6c79397ab15774837a0bc5755" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "thiserror" +version = "2.0.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09e52cb86a36cede5cb101bf8908837b3e4c6e5e59fe7fd85c23fb56200d189e" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "2.0.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fe5197923287db20a58125f0bc85c062f7f2c892de97b18c356f9efb14b28524" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.6", +] + +[[package]] +name = "unicode-ident" +version = "1.0.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d245f478577f809a851594d02313b640fb437e0bb33866753cff937863096954" + +[[package]] +name = "unicode-segmentation" +version = "1.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c6f5d3c3b1bf09027a88a6bc961fc00497d651009560b5463668dc81b0fa87a8" + +[[package]] +name = "unicode-width" +version = "0.1.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7dd6e30e90baa6f72411720665d41d89b9a3d039dc45b8faea1ddd07f617f6af" + +[[package]] +name = "unicode-width" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4ac048d71ede7ee76d585517add45da530660ef4390e49b098733c6e897f254" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + +[[package]] +name = "wasm-bindgen" +version = "0.2.129" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9bb54f33acc68fd454578d9820b0bde1a1a3d17aa17bb7b6595806d02886d409" +dependencies = [ + "cfg-if", + "once_cell", + "rustversion", + "wasm-bindgen-macro", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-macro" +version = "0.2.129" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2e29d0c35b16e224a7eeb5cd2d25e3e1968fbd65604117b44d3b789d00ee8535" +dependencies = [ + "quote", + "wasm-bindgen-macro-support", +] + +[[package]] +name = "wasm-bindgen-macro-support" +version = "0.2.129" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6f501a8bc3719dba86ef8ae4728879c08001bea749eb1333ac5b91e040e2a6b7" +dependencies = [ + "bumpalo", + "proc-macro2", + "quote", + "syn 3.0.6", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-shared" +version = "0.2.129" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23f0c9c52aa7cd7d77769a4cfe2a9adb1b331f489a41d912ce14513d5ab995c6" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + +[[package]] +name = "wyz" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05f360fc0b24296329c78fda852a1e9ae82de9cf7b27dae4b7f62f118f77b9ed" +dependencies = [ + "tap", +] + +[[package]] +name = "yansi" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cfe53a6657fd280eaa890a3bc59152892ffa3e30101319d168b781ed6529b049" + +[[package]] +name = "zerocopy" +version = "0.8.59" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6df92bf3d9227be3d53173901ddbffac2babc27ae50f397776ffd6dc33f800cb" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.8.59" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac4f328cf2f05d084e496c3e9c3f33ed0a183656a16e1fcec4d464d8373aec82" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] diff --git a/bench/stim-compare/Cargo.toml b/bench/stim-compare/Cargo.toml new file mode 100644 index 000000000..c46cc560a --- /dev/null +++ b/bench/stim-compare/Cargo.toml @@ -0,0 +1,23 @@ +[package] +name = "stim-compare" +version = "0.0.0" +edition = "2024" +publish = false +description = "Wall-clock comparison of ppvm's tableau against Stim's TableauSimulator on surface_d30" + +# Standalone: keeps the C++ `stim` dependency out of the ppvm workspace. +[workspace] + +[features] +default = ["traits-2"] +traits-2 = ["ppvm-stim/traits-2"] +legacy = ["ppvm-stim/legacy"] + +[dependencies] +ppvm-stim = { path = "../../crates/ppvm-stim", default-features = false } +bnum = { version = "0.13", features = ["numtraits"] } +rand = "0.10.1" +stim = "=0.4.5" + +[profile.release] +debug = "line-tables-only" diff --git a/bench/stim-compare/README.md b/bench/stim-compare/README.md new file mode 100644 index 000000000..ec48c6400 --- /dev/null +++ b/bench/stim-compare/README.md @@ -0,0 +1,348 @@ +# PR 204 tableau vs Stim + +Experiments comparing the `ppvm-tableau-2` generalized tableau (PR 204, +`codex/traits-2-impl`) with Stim's `TableauSimulator` on pure-Clifford circuits, +and the measurement changes on `david/batch-mr-reset` that came out of them. + +All numbers: Apple M-series (arm64), single thread, release builds, time per +shot of the circuit execution only (parsing and simulator construction are +outside the timer). On arm64 Stim uses plain 64-bit words (`bitword<64>`; it has +no NEON tableau path), so none of these comparisons involve SIMD. + +## Setup + +- **ppvm**: `stim-compare` (this crate, its own `[workspace]` so the C++ dependency + stays out of ppvm's). Runs `execute_validated_with_rng` on a fresh + `GeneralizedTableau` per shot. +- **Stim reference**: `stim_bench.py`, official PyPI `stim` 1.15.0 from + `ppvm-python/.venv`, one `TableauSimulator.do_circuit` call per shot. Use this, + not the Rust side's `stim` crate (0.4.5): that crate compiles Stim without + `NDEBUG`. In practice the difference was small (18.8 vs 18.2 ms on + `surface_d30`), but the PyPI build is what Stim users run. +- **surface_d30**: `crates/ppvm-stim/examples/surface_d30.stim` (1889 qubits, + 27 870 measurements, `REPEAT 29`). The harness also derives variants by + deleting instructions: `full-r1` (one repeated round), `no-noise`, + `no-measure` (drops `M`/`MR`/`R` and the records that reference them), + `gates-only`. + +```bash +STIM_RS_BUILD_FROM_SOURCE=1 cargo build --release +./target/release/stim-compare 20 ppvm [variant] # surface_d30 variants +./target/release/stim-compare file 50 # any circuit +../../ppvm-python/.venv/bin/python stim_bench.py 20 [variant] +../../ppvm-python/.venv/bin/python stim_bench.py file 50 +../../ppvm-python/.venv/bin/python sweep.py 3 # regenerates circuits/, writes sweep.csv +../../ppvm-python/.venv/bin/python attribute.py ... # writes attribution.csv +../../ppvm-python/.venv/bin/python nonclifford.py ... # ppvm only; writes nonclifford.csv +../../ppvm-python/.venv/bin/python gates.py 50 # per-gate cost; see gates.log +``` + +`file` mode takes an optional fifth argument, the qubit count; with `only = ppvm` +it then skips Stim's parse, so programs with `T` gates and rotations run. + +`attribute.py` expects one harness binary per commit in ``, named by +commit, each built from this crate with the workspace checked out at that commit. + +## Starting point: why PR 204 was 2× slower than Stim + +On `surface_d30` PR 204 took 37.6 ms/shot against Stim's 18.2 ms. Removing the +measurements made ppvm *faster* than Stim (10.1 vs 12.1 ms), so gates, noise +and SIMD were not the cause. The whole gap was measurement and reset, and +mostly the ~900 random ones (≈27 µs each against Stim's ≈7 µs). + +PR 204 already had Stim's complexity: deterministic measurements are O(n/W) +through the inverse tableau, random ones O(n²/W) in the worst case. The cost +was memory access. Profiled with samply: + +- `MR`/`R` ran one qubit at a time, column-major, so every random projection + gathered single generators bit by bit across all 1889 columns + (`project_inverse` → `gather_row`, ≈11.5 ms) plus a strided per-qubit pass in + the forward projection (≈6 ms). +- The batched `M` took a row guard up front, so even its deterministic targets + read their columns strided (`gather_column`, ≈6 ms). + +## Attribution + +Each commit on `david/batch-mr-reset`, timed in one session, alternating all +binaries and Stim over 3 rounds (median of round medians; PR 204 is `3befaf1b`). +This is the third such run; per-commit numbers agree with the earlier two +within a few percent. +Stim = official 1.15.0. Raw data: `attribution.csv`. + +| Workload | PR 204 | `2fa24ce4` | `d7d23c0e` | `4e8aeb5a` | `b4ea1ce9` | `295aaa8b` | `cfd25143` | `6e8bc36a` | Stim | +|---|---|---|---|---|---|---|---|---|---| +| surface_d30 full | 37.86 ms | 37.61 | **26.29** | **22.38** | 20.47 | **16.72** | 16.62 | **15.49** | 18.26 | +| surface_d30 1 round | 25.45 ms | 25.42 | **13.65** | **9.57** | 9.08 | **6.10** | 6.32 | 6.07 | 6.65 | +| surface_d30 no measure | 9.99 ms | 9.99 | 10.02 | 10.02 | 10.05 | 10.01 | 10.00 | **8.90** | 12.06 | +| surface d7 (118 q) | 160.1 µs | 160.2 | 140.8 | 138.1 | **110.9** | **84.4** | 83.8 | **72.2** | 61.4 | +| surface d11 (274 q) | 734.0 µs | 733.2 | 610.9 | 582.2 | **482.1** | **371.6** | 374.5 | **326.3** | 234.9 | +| surface d19 (778 q) | 6.32 ms | 6.32 | **4.51** | 4.05 | 3.58 | **2.75** | 2.76 | **2.38** | 2.45 | +| repetition d75 (149 q) | 972.0 µs | 975.7 | 983.7 | 1050.0 | **675.5** | **499.0** | 502.5 | **445.1** | 207.1 | +| color d31 (1081 q) | 59.27 ms | 59.51 | 52.35 | 52.18 | 50.46 | 50.02 | **8.11** | 7.13 | 7.08 | +| color d43 (2080 q) | 343.40 ms | 342.69 | 323.06 | 320.34 | 314.09 | 278.39 | **31.42** | 28.95 | 31.70 | + +Every commit leaves outcomes and RNG draws unchanged on these circuits; for +`295aaa8b`, `surface_d30` records were checked bit-identical against +`b4ea1ce9` over 22 seeded shots. As agreed, seeded results *can* change for +noisy `MR(p)` and multi-amplitude states (draw order), never the distribution. + +### `2fa24ce4` feat(tableau-2): Stim-style `measure_batch` + +API only: `GeneralizedTableau::measure_batch` checks which targets are random +in the canonical orientation (contiguous), collapses those under one row guard, +then measures the deterministic ones column-major — Stim's `collapse_z`. Nothing +calls it yet, so no change. + +### `d7d23c0e` perf(stim): batch `M`, `MR`, `R` — surface_d30 37.6 → 26.3 ms + +The traits-2 executor routes noise-free `M`, `MR` and `R`/`RZ` through +`measure_batch`, applying the `X` resets afterwards (gates need column-major; +`X_q` commutes with `Z_p`, so deferring is exact; repeated targets keep the old +loop). Random measurements now amortize one transpose per instruction and +contiguous projections; all-deterministic instructions never transpose. This is +the largest single win on `surface_d30` (−11.3 ms; −11.8 ms of the 1-round +variant). + +**Dead end, not committed:** batching under the *eager* row guard (transpose +every `MR`/`R` regardless) made `surface_d30` 123.6 ms/shot. Under the guard, +every measurement's column reads become strided, including the 29 all- +deterministic rounds. Checking determinism first, as Stim does, is what makes +batching pay. + +### `4e8aeb5a` perf(tableau-2): gather columns once — 26.3 → 22.4 ms + +A random measurement read the same two X columns at the measured qubit three +times (decomposition, `project_inverse`, `project_row_major`), each a strided +pass under the guard. Now gathered once and passed down. −3.9 ms on +`surface_d30`. `repetition_d75` (deterministic measurements only) is ~6% +slower with this commit, reproduced in three attribution runs; not +investigated (b4ea1ce9 more than recovers it). + +### `b4ea1ce9` perf(tableau-2): reuse buffers, skip unused masks — 22.4 → 20.5 ms + +Column buffers live in `MeasureScratch` instead of being allocated per +measurement, and a deterministic outcome on a state whose amplitudes are all at +index 0 (every Clifford run) skips widening the masks into the 2048-bit branch +index (one full-width shift and OR per set bit). Small on `surface_d30` +(−1.9 ms) but the biggest win for many-measurement, small-n circuits: +`repetition_d75` −36%, surface d7 −20%. + +### `295aaa8b` perf(tableau-2): Stim-style collapse on stabilizer states — 20.5 → 16.7 ms + +On a stabilizer state (one amplitude at index 0, inverse signs valid) +`measure_batch` now collapses like Stim's `collapse_qubit_z`: CX appends from +the pivot to the other anticommuting stabilizers, `S` if the accumulated +destabilizer anticommutes with `Z`, `H`, `X` if the outcome needs it. The +destabilizer column is never read and no other destabilizer is multiplied; a +deterministic outcome is one inverse-sign read. The new stabilizer is `±Z` +times other stabilizers rather than `±Z` itself, so the frame differs from the +textbook projection's while describing the same state. This is what puts ppvm +ahead of Stim on `surface_d30` (0.92×) and gives a further −24% on surface d7. + +### `cfd25143` perf(stim): batch `MX`, `MY`, `MRX`, `MRY`, noisy `M` — color d43 278.4 → 31.4 ms + +These still ran per qubit (`H`, one measurement, `H`) on the column-major +strided path, which profiling showed was 78% of color d31: Stim's color-code +generator ends with `MX`. The executor now rotates every target onto Z at once +around one `measure_batch` (`measure_in_basis`; a repeated target keeps the +per-target order), and `StimTableau::measure_noisy_many` flips each record +after the batch. Color d31 −84%, d43 −89%, which takes the large color codes +level with Stim; circuits without these instructions are unchanged. + +### `6e8bc36a` perf(tableau-2): trim gate kernels to live words — every circuit −4% to −14% + +A column's stride is rounded up to a whole 4-word block, so below 256 qubits +every gate kernel and the inverse-sign row product swept 4 words where only +`n.div_ceil(64)` can hold a set bit. The padding is zero and every kernel maps +zero words to zero words, so `gate1_mut` / `gate2_mut` and `inv_pair_phase` now +borrow only the live words, and the per-gate `get_disjoint_mut` overlap check +becomes a debug check on ranges that are disjoint by layout. `CX` at n ≤ 128 +goes from ~25 to ~15.5 ns (see "Per-gate cost"), surface d7 −14%, surface d19 +−14% (now ahead of Stim), `surface_d30` −7%, and the gates-only `no measure` +variant −11%. Tableaus are unchanged: mean outcomes over 200 seeded shots are +identical to `cfd25143` on five Clifford and three non-Clifford circuits. + +### `f5942998` perf(tableau-2): fuse `CX`'s forward update with its inverse-sign products + +A `CX` fetched its columns twice: once for the inverse-sign row products +(`ix_c·ix_t`, `iz_c·iz_t`), once per half for the forward update. Those are the +same eight columns, so `TableauData::cnot_fused` borrows them and both phase +planes once, reads the two `g`-rule terms, then applies the forward kernel. The +loops stay separate so they vectorize (one interleaved loop was 48% *slower* +at n=512), except at one live word (n ≤ 64), where a single scalar pass is +cheaper. Tableaus are unchanged: mean outcomes over 100–300 seeded shots are +identical to `6e8bc36a` on 8 Clifford and 4 non-Clifford circuits. + +Two-build run (`6e8bc36a` vs `f5942998` vs Stim, 3 alternating rounds) on a +machine that was not fully idle — Stim itself measured 3–5% slower than in the +runs above, so read the ratios, not the absolute times. Raw data: +`attribution-f5942998.csv`, `attribution-f5942998.log`. + +| Workload | `6e8bc36a` | `f5942998` | change | Stim | `f5942998` / Stim | +|---|---|---|---|---|---| +| surface_d30 full | 15.90 ms | 15.47 | −2.7% | 19.09 | 0.81× | +| surface_d30 1 round | 6.29 ms | 6.25 | −0.6% | 7.01 | 0.89× | +| surface_d30 no measure | 9.29 ms | 8.75 | −5.8% | 12.52 | 0.70× | +| surface d7 (118 q) | 73.5 µs | 67.9 | −7.6% | 62.6 | 1.08× | +| surface d11 (274 q) | 328.9 µs | 302.4 | −8.1% | 236.1 | 1.28× | +| surface d19 (778 q) | 2.46 ms | 2.23 | −9.3% | 2.54 | 0.88× | +| repetition d75 (149 q) | 450.5 µs | 409.6 | −9.1% | 207.1 | 1.98× | +| color d31 (1081 q) | 7.39 ms | 6.62 | −10.4% | 7.24 | 0.91× | +| color d43 (2080 q) | 29.92 ms | 28.68 | −4.1% | 32.70 | 0.88× | + +### `5f840789` perf(stim): sample noise Stim-style — repetition d25 −28% + +`DEPOLARIZE1`, `DEPOLARIZE2`, `X/Y/Z_ERROR` and `PAULI_CHANNEL_1` drew one +random number per target (and `DEPOLARIZE2` built and scanned a 15-entry table +per pair), where Stim's `RareErrorIterator` skips to the next error with a +geometric gap. The traits-2 adapter now does the same and then picks the +error's Pauli; lost qubits behave as before. Same distributions (statistical +executor test per channel and the `DEPOLARIZE2` pair correlation); seeded +results change. + +Against `f5942998`, alternating the two builds (2 rounds, both agree): + +| Workload | `f5942998` | `5f840789` | change | +|---|---|---|---| +| repetition d25 / d75 / d675 | 37.5 µs / 411.4 µs / 67.2 ms | 27.1 µs / 311.2 µs / 58.8 ms | −28% / −25% / −13% | +| surface d5 / d11 / d19 | 23.8 µs / 301.1 µs / 2.23 ms | 19.9 µs / 264.0 µs / 2.01 ms | −16% / −12% / −10% | +| color d9 / d31 | 85.4 µs / 6.63 ms | 74.7 µs / 6.13 ms | −13% / −8% | +| surface_d30 | 14.96 ms | 14.31 ms | −4% | + +Mean `1` outcomes now agree with Stim's compiled sampler within ≈1σ on +repetition d9/d25, surface d5/d7 and color d5/d9 (20 000 ppvm shots, 200 000 +Stim samples). The color codes' earlier ≈2σ excess is gone (color d5 +0.75σ, +color d9 +0.37σ), which may point at the old per-target sampler or may have +been chance. + +### `771daf64` perf(tableau-2): carry-save phase counters — `CX` −4% to −24% + +Keeping the inverse signs current is about half of a `CX` (a throwaway build +without it ran gate-only `surface_d30` in 4.1 ms instead of 7.7 ms; lazy signs +would not recover that for QEC circuits, which read the signs every round). +The cost is the `g`-rule row products, which ran two popcounts per word. +`PhaseCounter` keeps a 2-bit counter per bit position instead (Stim's +`cnt1` / `cnt2` update, checked exhaustively against the `g` term) and +popcounts once at the end. + +What the measurements forced: + +- A single serial counter does not vectorize: +21% at n=512, +9% on + `surface_d30`. From 8 words on, four counters over 4-word chunks; below + that one serial counter. Indexing the lanes as `i % 4` instead of by chunk + was +73% at n=512. +- Out of line in `cnot_fused` it cost more than it saved (surface d19 +3%); + `#[inline(always)]` fixed that. + +Against `5f840789`, alternating builds (2–3 rounds): + +| Workload | `5f840789` | `771daf64` | change | +|---|---|---|---| +| `CX` n = 64 / 274 / 778 / 2048 | 113.8 µs / 886 µs / 4.23 ms / 26.4 ms | 108.8 µs / 841 µs / 3.53 ms / 20.1 ms | −4% / −5% / −17% / −24% | +| surface d7 / d11 / d19 / d23 | 56.9 µs / 262.1 µs / 2.02 ms / 4.37 ms | 52.8 µs / 257.6 µs / 1.96 ms / 4.31 ms | −7% / −2% / −3% / −1% | +| repetition d75 | 311.4 µs | 266.7 µs | −14% | +| color d9 | 74.2 µs | 68.9 µs | −7% | +| surface_d30 | 14.38 ms | 13.21 ms | −8% | + +Outcomes are identical to `5f840789` on 8 Clifford and 4 non-Clifford circuits. + +## Size sweep + +`sweep.py` generates Stim memory circuits (`rounds = d`, all four noise +parameters 0.001): rotated surface code, repetition code, and color code +(`C_XYZ`, which ppvm-stim rejects, rewritten as `H` then `SQRT_X_DAG` for both +simulators), plus clifft-bench's `pure_surface_d7_r7`. Branch = `771daf64`. +Raw data: `sweep.csv`, `sweep.log`. Ratios are time / Stim time. + +| Surface | d=3 | 5 | 7 | 9 | 11 | 15 | 19 | 23 | 27 | 31 | +|---|---|---|---|---|---|---|---|---|---|---| +| qubits | 26 | 64 | 118 | 188 | 274 | 494 | 778 | 1126 | 1538 | 2014 | +| branch / Stim | 0.76 | 0.77 | 0.87 | 0.97 | 1.09 | 0.83 | 0.78 | 0.81 | 0.74 | 0.72 | +| PR 204 / Stim | 1.82 | 2.25 | 2.62 | 2.91 | 3.14 | 2.62 | 2.56 | 2.46 | 2.23 | 2.07 | + +| Repetition | d=3 | 9 | 25 | 75 | 225 | 675 | +|---|---|---|---|---|---|---| +| qubits | 5 | 17 | 49 | 149 | 449 | 1349 | +| branch / Stim | 0.50 | 0.78 | 1.10 | 1.28 | 0.96 | 0.95 | +| PR 204 / Stim | 1.33 | 2.96 | 4.54 | 4.74 | 3.39 | 2.30 | + +| Color | d=3 | 5 | 7 | 9 | 13 | 17 | 21 | 25 | 31 | 37 | 43 | +|---|---|---|---|---|---|---|---|---|---|---|---| +| qubits | 10 | 28 | 55 | 91 | 190 | 325 | 496 | 703 | 1081 | 1540 | 2080 | +| branch / Stim | 0.61 | 1.02 | 1.36 | 1.25 | 1.14 | 1.14 | 0.89 | 0.86 | 0.80 | 0.79 | 0.74 | +| PR 204 / Stim | 1.24 | 2.83 | 4.16 | 3.71 | 4.59 | 5.13 | 3.03 | 7.18 | 8.39 | 9.26 | 10.74 | + +Before `cfd25143` the color row was non-monotonic (d21 at 1.34× between d17 at +3.59× and d25 at 5.75×, `295aaa8b`): that was the unbatched `MX`. + +clifft-bench `pure_surface_d7_r7`: branch 52.5 µs, Stim 60.5 µs (0.87×). +Mean `1` outcomes per shot track Stim's on every circuit (e.g. surface d19 +2228.5 vs 2228.1, color d31 4929.1 vs 4926.2). + +## Non-Clifford regression check + +The fast collapse only fires on stabilizer states, but `4e8aeb5a` and +`b4ea1ce9` touch the general measurement path too. `nonclifford.py` times PR 204 +against `cfd25143` (3 alternating rounds) on ppvm's `cultivation_d5` and +clifft-bench's non-Clifford circuits, with clifft's half-turn `R_X(a)` / `U3` +rewritten as ppvm's `I[R_X(theta=a*pi)]` / `I[U3(...)]` tags for both builds. +Raw data: `nonclifford.csv`, `nonclifford.log`. + +| Program | qubits | PR 204 | `cfd25143` | ratio | +|---|---|---|---|---| +| cultivation_d5 | 42 | 5.98 ms | 5.99 ms | 1.00× | +| msc d3 | 15 | 83.8 µs | 82.0 µs | 0.98× | +| msc d5 | 42 | 5.99 ms | 6.01 ms | 1.00× | +| distillation | 85 | 201.9 µs | 193.7 µs | 0.96× | +| coherent d3 r1 | 26 | 285.3 µs | 280.6 µs | 0.98× | +| coherent d3 r3 | 26 | 4.29 ms | 4.30 ms | 1.00× | +| quantum volume q10 | 10 | 140.1 ms | 140.6 ms | 1.00× | + +No regressions. (`coherent_d5` and `quantum_volume_q20` take seconds to minutes +per shot and were left out.) + +## Per-gate cost + +Profiles of surface d11, repetition d75 and color d9 put `CX` at 39–55% of the +time at small n, with measurement down to 22–25%. `gates.py` isolates it: +`REPEAT 200` of a `CX` brickwork or an `H` layer, no noise or measurement. +Raw data: `gates.log`. + +| n | `CX` `cfd25143` | `6e8bc36a` | `f5942998` | Stim | `f5942998` / Stim | `H` `cfd25143` | `f5942998` | Stim | `f5942998` / Stim | +|---|---|---|---|---|---|---|---|---|---| +| 32 | 25.0 ns | 15.6 | 8.6 | 7.6 | 1.13× | 8.0 ns | 4.9 | 3.7 | 1.34× | +| 64 | 24.7 | 15.6 | 8.8 | 7.4 | 1.19× | 7.9 | 4.5 | 3.6 | 1.27× | +| 128 | 23.8 | 15.4 | 12.0 | 9.3 | 1.30× | 7.3 | 4.6 | 5.3 | 0.86× | +| 274 | 27.3 | 19.4 | 16.7 | 14.5 | 1.15× | 10.2 | 9.0 | 7.4 | 1.22× | +| 512 | 27.8 | 20.2 | 18.0 | 25.0 | 0.72× | 11.3 | 10.8 | 9.6 | 1.12× | +| 1024 | 40.7 | 34.5 | 33.7 | 41.2 | 0.82× | 15.6 | 13.1 | 14.9 | 0.88× | +| 2048 | 72.5 | 65.2 | 65.5 | 71.6 | 0.91× | 22.1 | 22.6 | 19.3 | 1.17× | + +The `f5942998` column (`gates.log`) was measured on the not-fully-idle machine; +`H` is unaffected by that commit, so its `f5942998` column doubles as a noise +check against `6e8bc36a` (4.6 / 4.5 / 4.4 / 8.4 / 10.6 / 12.9 / 21.7 ns there). +Before `6e8bc36a` ppvm's `CX` had a ~25 ns floor up to ~512 qubits: a fixed per-gate cost (the 4-word stride padding, swept twice per +`CX`, plus the per-gate overlap check). What remains at small n is mostly +structural: every gate updates the forward generator phases *and* the inverse +signs (`inv_pair_phase`, two row-phase products per `CX`), where Stim keeps only +its inverse tableau. The forward phases carry the multi-amplitude state +(`odd_phase_destabilizer_mask` in branching, case-a merge and expectations) and +the fallback when the inverse goes stale, so dropping either side is a redesign. + +## Open gaps + +1. **Mid-size circuits.** As of `771daf64` the branch is ahead of Stim on most + of the sweep; the remaining band is ≈50–330 qubits (1–6 words per column): + surface d11 1.09×, repetition d25/d75 1.10×/1.28×, color d7–d17 1.14–1.36×. + Fixed per-gate and per-measurement costs dominate there; the rest of the + per-gate gap is the dual phase bookkeeping (see "Per-gate cost"). `CZ` / `CY` + could get `CX`'s fusion, but only distillation, cultivation and MSC use `CZ` + and no benchmark uses `CY`. +2. **Measurement overhead at small n**: `measure_batch_one`, the per-target + 2048-bit "all amplitudes at index 0" check (`memcmp`), the determinism + pre-scan, `has_repeats`. +3. `4e8aeb5a`'s ~6% slowdown on deterministic-only measurement + (`repetition_d75`); more than recovered by `b4ea1ce9`. +4. Color codes sat ≈2σ above Stim's mean `1` count with the old per-target noise + sampler and ≈0.4–0.75σ with `5f840789`'s; a targeted test of the old + sampler would tell whether that was a bias or chance. diff --git a/bench/stim-compare/attribute.py b/bench/stim-compare/attribute.py new file mode 100644 index 000000000..a7e6ee552 --- /dev/null +++ b/bench/stim-compare/attribute.py @@ -0,0 +1,67 @@ +"""Per-commit attribution: time harness binaries built at each commit (see +README) against official Stim on a fixed workload set, alternating over rounds. + +Usage: `attribute.py ... [--rounds N]`. Writes `attribution.csv`. +""" + +import argparse +import csv +import pathlib +import re +import statistics +import subprocess + +HERE = pathlib.Path(__file__).resolve().parent +PYTHON = "/Users/david/git/ppvm/ppvm-python/.venv/bin/python" +# (label, harness args after the binary, stim_bench.py args) +WORKLOADS = [ + ("surface_d30 full", ["20", "ppvm", "full"], ["20", "full"]), + ("surface_d30 1 round", ["20", "ppvm", "full-r1"], ["20", "full-r1"]), + ("surface_d30 no measure", ["20", "ppvm", "no-measure"], ["20", "no-measure"]), + *( + (name, ["file", str(HERE / f"circuits/{name}.stim"), shots], ["file", str(HERE / f"circuits/{name}.stim"), shots]) + for name, shots in [ + ("surface_d7", "500"), + ("surface_d11", "500"), + ("surface_d19", "200"), + ("repetition_d75", "500"), + ("color_d31", "20"), + ("color_d43", "10"), + ] + ), +] + + +def execute_us(out: str) -> float: + m = re.search(r"execute median\s+([\d.]+)(ns|µs|ms|s)\b", out) + return float(m.group(1)) * {"ns": 1e-3, "µs": 1.0, "ms": 1e3, "s": 1e6}[m.group(2)] + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("bin_dir") + ap.add_argument("commits", nargs="+") + ap.add_argument("--rounds", type=int, default=3) + args = ap.parse_args() + tools = {c: [str(pathlib.Path(args.bin_dir) / c)] for c in args.commits} + tools["stim"] = [PYTHON, str(HERE / "stim_bench.py")] + + rows = [] + for label, harness, stim_args in WORKLOADS: + times = {t: [] for t in tools} + for _ in range(args.rounds): + for tool, cmd in tools.items(): + extra = stim_args if tool == "stim" else harness + out = subprocess.run([*cmd, *extra], capture_output=True, text=True, check=True).stdout + times[tool].append(execute_us(out)) + med = {t: statistics.median(v) for t, v in times.items()} + rows.append({"workload": label, **{f"{t}_us": round(v, 1) for t, v in med.items()}}) + print(label, " ".join(f"{t}={v:.1f}" for t, v in med.items()), flush=True) + with open(HERE / "attribution.csv", "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + +if __name__ == "__main__": + main() diff --git a/bench/stim-compare/attribution-f5942998.csv b/bench/stim-compare/attribution-f5942998.csv new file mode 100644 index 000000000..b8983dcbe --- /dev/null +++ b/bench/stim-compare/attribution-f5942998.csv @@ -0,0 +1,10 @@ +workload,6e8bc36a_us,f5942998_us,stim_us +surface_d30 full,15900.0,15470.0,19093.8 +surface_d30 1 round,6290.0,6250.0,7010.9 +surface_d30 no measure,9290.0,8750.0,12517.5 +surface_d7,73.5,67.9,62.6 +surface_d11,328.9,302.4,236.1 +surface_d19,2460.0,2230.0,2536.3 +repetition_d75,450.5,409.6,207.1 +color_d31,7390.0,6620.0,7237.4 +color_d43,29920.0,28680.0,32699.2 diff --git a/bench/stim-compare/attribution-f5942998.log b/bench/stim-compare/attribution-f5942998.log new file mode 100644 index 000000000..3f8512310 --- /dev/null +++ b/bench/stim-compare/attribution-f5942998.log @@ -0,0 +1,9 @@ +surface_d30 full 6e8bc36a=15900.0 f5942998=15470.0 stim=19093.8 +surface_d30 1 round 6e8bc36a=6290.0 f5942998=6250.0 stim=7010.9 +surface_d30 no measure 6e8bc36a=9290.0 f5942998=8750.0 stim=12517.5 +surface_d7 6e8bc36a=73.5 f5942998=67.9 stim=62.6 +surface_d11 6e8bc36a=328.9 f5942998=302.4 stim=236.1 +surface_d19 6e8bc36a=2460.0 f5942998=2230.0 stim=2536.3 +repetition_d75 6e8bc36a=450.5 f5942998=409.6 stim=207.1 +color_d31 6e8bc36a=7390.0 f5942998=6620.0 stim=7237.4 +color_d43 6e8bc36a=29920.0 f5942998=28680.0 stim=32699.2 diff --git a/bench/stim-compare/attribution.csv b/bench/stim-compare/attribution.csv new file mode 100644 index 000000000..10a9b7062 --- /dev/null +++ b/bench/stim-compare/attribution.csv @@ -0,0 +1,10 @@ +workload,3befaf1b_us,2fa24ce4_us,d7d23c0e_us,4e8aeb5a_us,b4ea1ce9_us,295aaa8b_us,cfd25143_us,6e8bc36a_us,stim_us +surface_d30 full,37860.0,37610.0,26290.0,22380.0,20470.0,16720.0,16620.0,15490.0,18260.6 +surface_d30 1 round,25450.0,25420.0,13650.0,9570.0,9080.0,6100.0,6320.0,6070.0,6654.1 +surface_d30 no measure,9990.0,9990.0,10020.0,10020.0,10050.0,10010.0,10000.0,8900.0,12058.3 +surface_d7,160.1,160.2,140.8,138.1,110.9,84.4,83.8,72.2,61.4 +surface_d11,734.0,733.2,610.9,582.2,482.1,371.6,374.5,326.3,234.9 +surface_d19,6320.0,6320.0,4510.0,4050.0,3580.0,2750.0,2760.0,2380.0,2454.0 +repetition_d75,972.0,975.7,983.7,1050.0,675.5,499.0,502.5,445.1,207.1 +color_d31,59270.0,59510.0,52350.0,52180.0,50460.0,50020.0,8110.0,7130.0,7080.0 +color_d43,343400.0,342690.0,323060.0,320340.0,314090.0,278390.0,31420.0,28950.0,31697.6 diff --git a/bench/stim-compare/attribution.log b/bench/stim-compare/attribution.log new file mode 100644 index 000000000..718217050 --- /dev/null +++ b/bench/stim-compare/attribution.log @@ -0,0 +1,9 @@ +surface_d30 full 3befaf1b=37860.0 2fa24ce4=37610.0 d7d23c0e=26290.0 4e8aeb5a=22380.0 b4ea1ce9=20470.0 295aaa8b=16720.0 cfd25143=16620.0 6e8bc36a=15490.0 stim=18260.6 +surface_d30 1 round 3befaf1b=25450.0 2fa24ce4=25420.0 d7d23c0e=13650.0 4e8aeb5a=9570.0 b4ea1ce9=9080.0 295aaa8b=6100.0 cfd25143=6320.0 6e8bc36a=6070.0 stim=6654.1 +surface_d30 no measure 3befaf1b=9990.0 2fa24ce4=9990.0 d7d23c0e=10020.0 4e8aeb5a=10020.0 b4ea1ce9=10050.0 295aaa8b=10010.0 cfd25143=10000.0 6e8bc36a=8900.0 stim=12058.3 +surface_d7 3befaf1b=160.1 2fa24ce4=160.2 d7d23c0e=140.8 4e8aeb5a=138.1 b4ea1ce9=110.9 295aaa8b=84.4 cfd25143=83.8 6e8bc36a=72.2 stim=61.4 +surface_d11 3befaf1b=734.0 2fa24ce4=733.2 d7d23c0e=610.9 4e8aeb5a=582.2 b4ea1ce9=482.1 295aaa8b=371.6 cfd25143=374.5 6e8bc36a=326.3 stim=234.9 +surface_d19 3befaf1b=6320.0 2fa24ce4=6320.0 d7d23c0e=4510.0 4e8aeb5a=4050.0 b4ea1ce9=3580.0 295aaa8b=2750.0 cfd25143=2760.0 6e8bc36a=2380.0 stim=2454.0 +repetition_d75 3befaf1b=972.0 2fa24ce4=975.7 d7d23c0e=983.7 4e8aeb5a=1050.0 b4ea1ce9=675.5 295aaa8b=499.0 cfd25143=502.5 6e8bc36a=445.1 stim=207.1 +color_d31 3befaf1b=59270.0 2fa24ce4=59510.0 d7d23c0e=52350.0 4e8aeb5a=52180.0 b4ea1ce9=50460.0 295aaa8b=50020.0 cfd25143=8110.0 6e8bc36a=7130.0 stim=7080.0 +color_d43 3befaf1b=343400.0 2fa24ce4=342690.0 d7d23c0e=323060.0 4e8aeb5a=320340.0 b4ea1ce9=314090.0 295aaa8b=278390.0 cfd25143=31420.0 6e8bc36a=28950.0 stim=31697.6 diff --git a/bench/stim-compare/gates.log b/bench/stim-compare/gates.log new file mode 100644 index 000000000..f0baa53be --- /dev/null +++ b/bench/stim-compare/gates.log @@ -0,0 +1,15 @@ +circuit gates ppvm ns/gate stim ns/gate ratio +CX_n32 6200 8.59 7.60 1.13 +CX_n64 12600 8.83 7.40 1.19 +CX_n128 25400 12.04 9.26 1.30 +CX_n274 54600 16.74 14.50 1.15 +CX_n512 102200 18.00 25.03 0.72 +CX_n1024 204600 33.72 41.21 0.82 +CX_n2048 409400 65.49 71.60 0.91 +H_n32 6400 4.93 3.67 1.34 +H_n64 12800 4.54 3.57 1.27 +H_n128 25600 4.56 5.32 0.86 +H_n274 54800 9.04 7.38 1.22 +H_n512 102400 10.84 9.64 1.12 +H_n1024 204800 13.13 14.93 0.88 +H_n2048 409600 22.56 19.27 1.17 diff --git a/bench/stim-compare/gates.py b/bench/stim-compare/gates.py new file mode 100644 index 000000000..d6b24ae52 --- /dev/null +++ b/bench/stim-compare/gates.py @@ -0,0 +1,56 @@ +"""Per-gate cost: `REPEAT 200` of a gate layer on n qubits, ppvm vs official Stim. + +`CX` is a brickwork (pairs (0,1),(2,3),... then (1,2),(3,4),...), `H` hits +every qubit. Prints ns per gate. Usage: `gates.py [shots]`. +""" + +import pathlib +import re +import subprocess +import sys + +HERE = pathlib.Path(__file__).resolve().parent +PYTHON = "/Users/david/git/ppvm/ppvm-python/.venv/bin/python" +SIZES = [32, 64, 128, 274, 512, 1024, 2048] +REPS = 200 + + +def layers(gate: str, n: int) -> tuple[str, int]: + """(circuit, gates per shot).""" + if gate == "CX": + even = list(range(0, n - n % 2)) + odd = list(range(1, n - 1 - (n - 1) % 2 + 1))[: (n - 1) // 2 * 2] + body = f"CX {' '.join(map(str, even))}\nCX {' '.join(map(str, odd))}\n" + count = (len(even) + len(odd)) // 2 + else: + body = f"H {' '.join(map(str, range(n)))}\n" + count = n + return f"REPEAT {REPS} {{\n{body}}}\n", REPS * count + + +def execute_us(out: str) -> float: + m = re.search(r"execute median\s+([\d.]+)(ns|µs|ms|s)\b", out) + return float(m.group(1)) * {"ns": 1e-3, "µs": 1.0, "ms": 1e3, "s": 1e6}[m.group(2)] + + +def main() -> None: + shots = sys.argv[1] if len(sys.argv) > 1 else "50" + out_dir = HERE / "circuits/gates" + out_dir.mkdir(parents=True, exist_ok=True) + print(f"{'circuit':10s} {'gates':>7s} {'ppvm ns/gate':>13s} {'stim ns/gate':>13s} {'ratio':>6s}") + for gate in ["CX", "H"]: + for n in SIZES: + text, count = layers(gate, n) + path = out_dir / f"{gate}_n{n}.stim" + path.write_text(text) + ppvm = execute_us(subprocess.run( + [str(HERE / "target/release/stim-compare"), "file", str(path), shots, "ppvm", str(n)], + capture_output=True, text=True, check=True).stdout) + stim = execute_us(subprocess.run( + [PYTHON, str(HERE / "stim_bench.py"), "file", str(path), shots], + capture_output=True, text=True, check=True).stdout) + print(f"{gate}_n{n:<6d} {count:7d} {1e3 * ppvm / count:13.2f} {1e3 * stim / count:13.2f} {ppvm / stim:6.2f}") + + +if __name__ == "__main__": + main() diff --git a/bench/stim-compare/nonclifford.csv b/bench/stim-compare/nonclifford.csv new file mode 100644 index 000000000..b63a77963 --- /dev/null +++ b/bench/stim-compare/nonclifford.csv @@ -0,0 +1,8 @@ +program,qubits,shots,3befaf1b_us,cfd25143_us,ratio +cultivation_d5,42,166,5980.0,5990.0,1.002 +msc_d3_inject_cultivate_p1e-3,15,500,83.8,82.0,0.979 +msc_d5_inject_cultivate_p1e-3,42,161,5990.0,6010.0,1.003 +distillation,85,500,201.9,193.7,0.959 +coherent_d3_r1,26,500,285.3,280.6,0.984 +coherent_d3_r3,26,230,4290.0,4300.0,1.002 +quantum_volume_q10_seed42,10,7,140060.0,140610.0,1.004 diff --git a/bench/stim-compare/nonclifford.log b/bench/stim-compare/nonclifford.log new file mode 100644 index 000000000..5422a300f --- /dev/null +++ b/bench/stim-compare/nonclifford.log @@ -0,0 +1,7 @@ +cultivation_d5 n= 42 shots=166 3befaf1b 5980.0µs cfd25143 5990.0µs cfd25143/3befaf1b 1.00x +msc_d3_inject_cultivate_p1e-3 n= 15 shots=500 3befaf1b 83.8µs cfd25143 82.0µs cfd25143/3befaf1b 0.98x +msc_d5_inject_cultivate_p1e-3 n= 42 shots=161 3befaf1b 5990.0µs cfd25143 6010.0µs cfd25143/3befaf1b 1.00x +distillation n= 85 shots=500 3befaf1b 201.9µs cfd25143 193.7µs cfd25143/3befaf1b 0.96x +coherent_d3_r1 n= 26 shots=500 3befaf1b 285.3µs cfd25143 280.6µs cfd25143/3befaf1b 0.98x +coherent_d3_r3 n= 26 shots=230 3befaf1b 4290.0µs cfd25143 4300.0µs cfd25143/3befaf1b 1.00x +quantum_volume_q10_seed42 n= 10 shots= 7 3befaf1b 140060.0µs cfd25143 140610.0µs cfd25143/3befaf1b 1.00x diff --git a/bench/stim-compare/nonclifford.py b/bench/stim-compare/nonclifford.py new file mode 100644 index 000000000..ccc6fa472 --- /dev/null +++ b/bench/stim-compare/nonclifford.py @@ -0,0 +1,110 @@ +"""Regression check on non-Clifford programs (ppvm only; Stim can't run them): +time harness binaries built at two commits, alternating over rounds. + +Usage: `nonclifford.py ... [--rounds N]`. Writes `nonclifford.csv`. +""" + +import argparse +import csv +import pathlib +import re +import statistics +import subprocess + +HERE = pathlib.Path(__file__).resolve().parent +CLIFFT = pathlib.Path.home() / "git/clifft-bench/workloads/circuits" +PROGRAMS = [ + HERE / "../../crates/ppvm-stim/tests/data/cultivation_d5.stim", + *(CLIFFT / f"{name}.stim" for name in [ + "msc_d3_inject_cultivate_p1e-3", + "msc_d5_inject_cultivate_p1e-3", + "distillation", + "coherent_d3_r1", + "coherent_d3_r3", + "quantum_volume_q10_seed42", + ]), +] +ANNOTATIONS = {"DETECTOR", "OBSERVABLE_INCLUDE", "QUBIT_COORDS", "SHIFT_COORDS", "TICK", "REPEAT"} + + +def to_ppvm_dialect(text: str) -> str: + """Clifft's `R_X(a)` / `U3(t, p, l)` (half-turns) as ppvm's `I[...]` tags (radians).""" + + def rot(m: re.Match) -> str: + return f"I[R_{m.group(1)}(theta={m.group(2).strip()}*pi)]" + + def u3(m: re.Match) -> str: + t, p, l = (a.strip() for a in m.group(1).split(",")) + return f"I[U3(theta={t}*pi, phi={p}*pi, lambda={l}*pi)]" + + text = re.sub(r"\bR_([XYZ])\(([^)]*)\)", rot, text) + return re.sub(r"\bU3\(([^)]*)\)", u3, text) + + +def n_qubits(text: str) -> int: + """Largest qubit target plus one, without a Stim parse (Stim rejects `T`).""" + top = 0 + for line in text.splitlines(): + line = line.split("#")[0].strip() + name = re.split(r"[\s(]", line, maxsplit=1)[0] + if not line or name in ANNOTATIONS or line == "}": + continue + targets = re.sub(r"\([^)]*\)", " ", line[len(name):]) + for token in re.split(r"[\s*]+", targets): + token = token.lstrip("!").lstrip("XYZ") + if token.isdigit(): + top = max(top, int(token) + 1) + return max(top, 1) + + +def execute_us(out: str) -> float: + m = re.search(r"execute median\s+([\d.]+)(ns|µs|ms|s)\b", out) + return float(m.group(1)) * {"ns": 1e-3, "µs": 1.0, "ms": 1e3, "s": 1e6}[m.group(2)] + + +def run(binary: pathlib.Path, program: pathlib.Path, shots: int, n: int) -> float: + out = subprocess.run( + [str(binary), "file", str(program), str(shots), "ppvm", str(n)], + capture_output=True, text=True, check=True, + ).stdout + return execute_us(out) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("bin_dir") + ap.add_argument("commits", nargs="+") + ap.add_argument("--rounds", type=int, default=3) + args = ap.parse_args() + bins = {c: pathlib.Path(args.bin_dir) / c for c in args.commits} + + out_dir = HERE / "circuits/nonclifford" + out_dir.mkdir(parents=True, exist_ok=True) + rows = [] + for source in PROGRAMS: + program = out_dir / source.name + program.write_text(to_ppvm_dialect(source.read_text())) + n = n_qubits(program.read_text()) + pilot = run(bins[args.commits[-1]], program, 2, n) + shots = max(3, min(500, int(1.0e6 / max(pilot, 1.0)))) + times = {c: [] for c in bins} + for _ in range(args.rounds): + for c, b in bins.items(): + times[c].append(run(b, program, shots, n)) + med = {c: statistics.median(v) for c, v in times.items()} + first, last = args.commits[0], args.commits[-1] + row = {"program": program.stem, "qubits": n, "shots": shots, + **{f"{c}_us": round(v, 1) for c, v in med.items()}, + "ratio": round(med[last] / med[first], 3)} + rows.append(row) + print(f"{program.stem:34s} n={n:3d} shots={shots:3d} " + + " ".join(f"{c} {v:10.1f}µs" for c, v in med.items()) + + f" {last}/{first} {row['ratio']:.2f}x", flush=True) + with open(HERE / "nonclifford.csv", "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + +if __name__ == "__main__": + main() diff --git a/bench/stim-compare/src/main.rs b/bench/stim-compare/src/main.rs new file mode 100644 index 000000000..33db347ed --- /dev/null +++ b/bench/stim-compare/src/main.rs @@ -0,0 +1,160 @@ +//! Time one surface_d30 shot on ppvm's `GeneralizedTableau` and on Stim's +//! `TableauSimulator`, interleaved so machine drift hits both sides equally. +//! +//! Usage: `stim-compare [shots] [only] [variant]` where `only` is `ppvm`, `stim` +//! or `both`, and `variant` is one of [`VARIANTS`] (default: all of them). +//! `stim-compare file [shots] [only] [n_qubits]` times one `.stim` file +//! instead; with `only = ppvm` and `n_qubits` given, Stim never parses it, so +//! non-Clifford programs work. + +use std::time::{Duration, Instant}; + +use bnum::types::U2048; +#[cfg(feature = "legacy")] +use ppvm_stim::backend::config::indexmap::ByteFxHashF64; +use ppvm_stim::backend::prelude::*; +use ppvm_stim::{ExtendedProgram, execute_validated_with_rng, parse_extended, validate}; +use rand::SeedableRng; + +#[cfg(feature = "legacy")] +type Tab = GeneralizedTableau, U2048>; +#[cfg(not(feature = "legacy"))] +type Tab = GeneralizedTableau; + +const SRC: &str = include_str!("../../../crates/ppvm-stim/examples/surface_d30.stim"); +const WARMUP: usize = 2; + +const MEASURE_OPS: &[&str] = &["M", "MR", "R", "DETECTOR", "OBSERVABLE_INCLUDE"]; +const NOISE_OPS: &[&str] = &["DEPOLARIZE1", "DEPOLARIZE2", "X_ERROR"]; + +/// (name, ops stripped from surface_d30, `REPEAT 29` count). Records go with +/// the measurements because `rec[-k]` targets would dangle. +const VARIANTS: &[(&str, &[&str], &[&str], u32)] = &[ + ("full", &[], &[], 29), + ("full-r1", &[], &[], 1), + ("no-noise", NOISE_OPS, &[], 29), + ("no-measure", MEASURE_OPS, &[], 29), + ("no-measure-r1", MEASURE_OPS, &[], 1), + ("gates-only", MEASURE_OPS, NOISE_OPS, 29), +]; + +/// Drop every line whose instruction name is in `a` or `b`; set the round count. +fn strip(src: &str, a: &[&str], b: &[&str], reps: u32) -> String { + src.lines() + .filter(|line| { + let name = line.trim_start().split(['(', ' ']).next().unwrap_or(""); + !a.contains(&name) && !b.contains(&name) + }) + .map(|line| match line { + "REPEAT 29 {" => format!("REPEAT {reps} {{\n"), + _ => format!("{line}\n"), + }) + .collect() +} + +/// Per shot: (construct time, execute time, number of `1` outcomes). +type Shot = (Duration, Duration, usize); + +fn ppvm_shot(prog: &ExtendedProgram, n_qubits: usize, seed: u64) -> Shot { + let t0 = Instant::now(); + let mut tab = Tab::new(n_qubits, 1e-10); + let mut rng = rand::rngs::SmallRng::seed_from_u64(seed); + let mut rec = Vec::with_capacity(prog.measurement_count()); + let t1 = Instant::now(); + execute_validated_with_rng(&prog.instructions, &mut tab, &mut rec, &mut rng); + let t2 = Instant::now(); + (t1 - t0, t2 - t1, rec.iter().filter(|&&b| b == Some(true)).count()) +} + +fn stim_shot(circuit: &stim::Circuit, n_qubits: usize, seed: u64) -> Shot { + let t0 = Instant::now(); + let mut sim = stim::TableauSimulator::with_seed(seed); + // Pre-size so do_circuit's ensure_large_enough_for_qubits is a no-op. + sim.set_num_qubits(n_qubits); + let t1 = Instant::now(); + sim.do_circuit(circuit); + let t2 = Instant::now(); + let ones = sim.current_measurement_record().iter().filter(|&&b| b).count(); + (t1 - t0, t2 - t1, ones) +} + +fn median(times: &mut [Duration]) -> Duration { + times.sort(); + times[times.len() / 2] +} + +fn stats(name: &str, shots: &[Shot], n_meas: usize) -> Duration { + let mut build: Vec<_> = shots.iter().map(|s| s.0).collect(); + let mut exec: Vec<_> = shots.iter().map(|s| s.1).collect(); + let mean_ones = shots.iter().map(|s| s.2).sum::() as f64 / shots.len() as f64; + let (b, e) = (median(&mut build), median(&mut exec)); + println!( + "{name:>5}: construct {b:>9.2?} execute median {e:>9.2?} (min {:>9.2?}, max {:>9.2?}) | mean ones/shot {mean_ones:.3} of {n_meas}", + exec[0], + exec[exec.len() - 1], + ); + e +} + +fn run_variant( + name: &str, + src: &str, + shots: usize, + run_ppvm: bool, + run_stim: bool, + n_qubits: Option, +) { + let prog = parse_extended(src).expect("ppvm parse"); + validate(&prog).expect("ppvm validate"); + let n_meas = prog.measurement_count(); + let circuit: Option = + (run_stim || n_qubits.is_none()).then(|| src.parse().expect("stim parse")); + let n_qubits = n_qubits.unwrap_or_else(|| circuit.as_ref().unwrap().num_qubits().max(1)); + if let Some(circuit) = &circuit { + assert_eq!(n_meas as u64, circuit.num_measurements()); + } + println!("\n[{name}] {n_qubits} qubits, {n_meas} measurements"); + + let (mut ps, mut ss) = (vec![], vec![]); + for i in 0..WARMUP + shots { + let seed = i as u64; + let p = run_ppvm.then(|| ppvm_shot(&prog, n_qubits, seed)); + let s = run_stim.then(|| stim_shot(circuit.as_ref().unwrap(), n_qubits, seed)); + if i >= WARMUP { + ps.extend(p); + ss.extend(s); + } + } + + let p = run_ppvm.then(|| stats("ppvm", &ps, n_meas)); + let s = run_stim.then(|| stats("stim", &ss, n_meas)); + if let (Some(p), Some(s)) = (p, s) { + println!("execute median ppvm / stim = {:.2}x", p.as_secs_f64() / s.as_secs_f64()); + } +} + +fn main() { + let args: Vec = std::env::args().collect(); + let backend = if cfg!(feature = "legacy") { "legacy" } else { "traits-2" }; + if args.get(1).map(String::as_str) == Some("file") { + let path = args.get(2).expect("file path"); + let shots: usize = args.get(3).map_or(20, |s| s.parse().expect("shots")); + let only = args.get(4).map_or("ppvm", String::as_str); + let n_qubits = args.get(5).map(|s| s.parse().expect("n_qubits")); + let src = std::fs::read_to_string(path).expect("read circuit"); + println!("{path}: ppvm backend = {backend}, {shots} shots"); + run_variant("file", &src, shots, only != "stim", only != "ppvm", n_qubits); + return; + } + let shots: usize = args.get(1).map_or(20, |s| s.parse().expect("shots")); + let only = args.get(2).map_or("both", String::as_str); + let variant = args.get(3).map(String::as_str); + println!("surface_d30: ppvm backend = {backend}, {shots} shots"); + + for &(name, a, b, reps) in VARIANTS { + if variant.is_none_or(|v| v == name) { + let src = strip(SRC, a, b, reps); + run_variant(name, &src, shots, only != "stim", only != "ppvm", None); + } + } +} diff --git a/bench/stim-compare/stim_bench.py b/bench/stim-compare/stim_bench.py new file mode 100644 index 000000000..8c73b7baf --- /dev/null +++ b/bench/stim-compare/stim_bench.py @@ -0,0 +1,82 @@ +"""Time Stim's TableauSimulator (official PyPI build) on the same surface_d30 +variants as the Rust harness, one shot per `do_circuit` call. + +Usage: `stim_bench.py [shots] [variant]`, or `stim_bench.py file [shots]`. +""" + +import pathlib +import statistics +import sys +import time + +import stim + +SRC = ( + pathlib.Path(__file__).parent / "../../crates/ppvm-stim/examples/surface_d30.stim" +).read_text() +WARMUP = 2 + +MEASURE_OPS = {"M", "MR", "R", "DETECTOR", "OBSERVABLE_INCLUDE"} +NOISE_OPS = {"DEPOLARIZE1", "DEPOLARIZE2", "X_ERROR"} + +# Same (name, stripped ops, REPEAT count) table as src/main.rs. +VARIANTS = [ + ("full", set(), 29), + ("full-r1", set(), 1), + ("no-noise", NOISE_OPS, 29), + ("no-measure", MEASURE_OPS, 29), + ("no-measure-r1", MEASURE_OPS, 1), + ("gates-only", MEASURE_OPS | NOISE_OPS, 29), +] + + +def strip(src: str, ops: set[str], reps: int) -> str: + out = [] + for line in src.splitlines(): + name = line.lstrip().replace("(", " ").split(" ")[0] + if name in ops: + continue + out.append(f"REPEAT {reps} {{" if line == "REPEAT 29 {" else line) + return "\n".join(out) + "\n" + + +def run_variant(name: str, src: str, shots: int) -> None: + circuit = stim.Circuit(src) + n_qubits = max(circuit.num_qubits, 1) + build, execute, ones = [], [], [] + for seed in range(WARMUP + shots): + t0 = time.perf_counter() + sim = stim.TableauSimulator(seed=seed) + sim.set_num_qubits(n_qubits) + t1 = time.perf_counter() + sim.do_circuit(circuit) + t2 = time.perf_counter() + if seed >= WARMUP: + build.append(t1 - t0) + execute.append(t2 - t1) + ones.append(sum(sim.current_measurement_record())) + ms = lambda s: f"{1e6 * s:10.1f}µs" + print(f"\n[{name}] {n_qubits} qubits, {circuit.num_measurements} measurements") + print( + f" stim: construct {ms(statistics.median(build))} execute median " + f"{ms(statistics.median(execute))} (min {ms(min(execute))}, max {ms(max(execute))})" + f" | mean ones/shot {statistics.mean(ones):.1f}" + ) + + +def main() -> None: + if len(sys.argv) > 2 and sys.argv[1] == "file": + shots = int(sys.argv[3]) if len(sys.argv) > 3 else 20 + print(f"{sys.argv[2]}: stim {stim.__version__}, {shots} shots") + run_variant("file", pathlib.Path(sys.argv[2]).read_text(), shots) + return + shots = int(sys.argv[1]) if len(sys.argv) > 1 else 20 + only = sys.argv[2] if len(sys.argv) > 2 else None + print(f"surface_d30: stim {stim.__version__}, {shots} shots") + for name, ops, reps in VARIANTS: + if only in (None, name): + run_variant(name, strip(SRC, ops, reps), shots) + + +if __name__ == "__main__": + main() diff --git a/bench/stim-compare/sweep.csv b/bench/stim-compare/sweep.csv new file mode 100644 index 000000000..0712e1997 --- /dev/null +++ b/bench/stim-compare/sweep.csv @@ -0,0 +1,29 @@ +family,d,qubits,measurements,shots,branch_us,pr204_us,stim_us,branch_over_stim,pr204_over_stim,branch_ones,pr204_ones,stim_ones +clifft_surface,7,118,385,500,52.5,157.7,60.5,0.868,2.606,115.62,113.2,113.3 +surface,3,26,33,500,5.3,12.7,7.0,0.756,1.816,10.594,11.0,10.7 +surface,5,64,145,500,19.2,56.1,25.0,0.77,2.245,43.802,44.3,43.2 +surface,7,118,385,500,52.0,157.6,60.1,0.866,2.622,115.62,113.2,113.3 +surface,9,188,801,500,119.2,356.2,122.4,0.974,2.91,239.086,239.1,238.4 +surface,11,274,1441,500,252.2,725.0,230.6,1.094,3.144,432.324,429.4,430.7 +surface,15,494,3585,500,665.3,2090.0,797.3,0.834,2.621,1086.514,1090.8,1083.6 +surface,19,778,7201,500,1900.0,6210.0,2424.0,0.784,2.562,2228.538,2230.1,2228.1 +surface,23,1126,12673,240,4160.0,12710.0,5156.9,0.807,2.465,3988.512,4009.5,4001.0 +surface,27,1538,20385,122,8120.0,24400.0,10943.5,0.742,2.23,6556.066,6574.7,6594.4 +surface,31,2014,30721,69,14440.0,41670.0,20137.1,0.717,2.069,10046.913,10164.4,10185.7 +repetition,3,5,9,500,0.8,2.0,1.5,0.5,1.333,0.066,0.1,0.1 +repetition,9,17,81,500,4.0,15.1,5.1,0.776,2.965,1.574,1.4,1.5 +repetition,25,49,625,500,26.3,108.2,23.8,1.105,4.545,27.908,26.6,28.2 +repetition,75,149,5625,500,261.5,968.0,204.2,1.281,4.74,634.412,639.6,639.2 +repetition,225,449,50625,364,2770.0,9810.0,2894.6,0.957,3.389,12564.223,12725.6,12624.1 +repetition,675,1349,455625,19,50380.0,121960.0,53024.1,0.95,2.3,180786.105,181673.2,180729.1 +color,3,10,16,500,4.0,8.0,6.5,0.615,1.237,6.626,6.6,6.6 +color,5,28,64,500,13.0,36.2,12.8,1.019,2.832,27.306,27.5,27.7 +color,7,55,163,500,32.3,98.9,23.8,1.358,4.155,65.686,65.6,65.5 +color,9,91,331,500,67.2,199.6,53.8,1.249,3.71,128.408,129.1,127.8 +color,13,190,946,500,213.1,859.0,187.1,1.139,4.591,373.976,376.5,375.9 +color,17,325,2053,500,582.2,2630.0,512.8,1.135,5.129,825.632,828.8,829.1 +color,21,496,3796,500,1090.0,3710.0,1225.4,0.89,3.028,1510.2,1519.0,1514.8 +color,25,703,6319,414,2380.0,19970.0,2781.1,0.856,7.181,2571.821,2572.6,2577.8 +color,31,1081,11881,178,5620.0,59310.0,7068.0,0.795,8.391,4929.101,4930.0,4926.2 +color,37,1540,20008,77,12560.0,146490.0,15826.9,0.794,9.256,8393.558,8416.1,8397.3 +color,43,2080,31186,42,23670.0,344410.0,32056.8,0.738,10.744,13334.0,13370.4,13222.0 diff --git a/bench/stim-compare/sweep.log b/bench/stim-compare/sweep.log new file mode 100644 index 000000000..61ca74712 --- /dev/null +++ b/bench/stim-compare/sweep.log @@ -0,0 +1,28 @@ +clifft_surface d= 7 n= 118 meas= 385 shots=500 branch 52.5µs pr204 157.7µs stim 60.5µs branch/stim 0.87x pr204/stim 2.61x +surface d= 3 n= 26 meas= 33 shots=500 branch 5.3µs pr204 12.7µs stim 7.0µs branch/stim 0.76x pr204/stim 1.82x +surface d= 5 n= 64 meas= 145 shots=500 branch 19.2µs pr204 56.1µs stim 25.0µs branch/stim 0.77x pr204/stim 2.25x +surface d= 7 n= 118 meas= 385 shots=500 branch 52.0µs pr204 157.6µs stim 60.1µs branch/stim 0.87x pr204/stim 2.62x +surface d= 9 n= 188 meas= 801 shots=500 branch 119.2µs pr204 356.2µs stim 122.4µs branch/stim 0.97x pr204/stim 2.91x +surface d= 11 n= 274 meas= 1441 shots=500 branch 252.2µs pr204 725.0µs stim 230.6µs branch/stim 1.09x pr204/stim 3.14x +surface d= 15 n= 494 meas= 3585 shots=500 branch 665.3µs pr204 2090.0µs stim 797.3µs branch/stim 0.83x pr204/stim 2.62x +surface d= 19 n= 778 meas= 7201 shots=500 branch 1900.0µs pr204 6210.0µs stim 2424.0µs branch/stim 0.78x pr204/stim 2.56x +surface d= 23 n= 1126 meas= 12673 shots=240 branch 4160.0µs pr204 12710.0µs stim 5156.9µs branch/stim 0.81x pr204/stim 2.46x +surface d= 27 n= 1538 meas= 20385 shots=122 branch 8120.0µs pr204 24400.0µs stim 10943.5µs branch/stim 0.74x pr204/stim 2.23x +surface d= 31 n= 2014 meas= 30721 shots= 69 branch 14440.0µs pr204 41670.0µs stim 20137.1µs branch/stim 0.72x pr204/stim 2.07x +repetition d= 3 n= 5 meas= 9 shots=500 branch 0.8µs pr204 2.0µs stim 1.5µs branch/stim 0.50x pr204/stim 1.33x +repetition d= 9 n= 17 meas= 81 shots=500 branch 4.0µs pr204 15.1µs stim 5.1µs branch/stim 0.78x pr204/stim 2.96x +repetition d= 25 n= 49 meas= 625 shots=500 branch 26.3µs pr204 108.2µs stim 23.8µs branch/stim 1.10x pr204/stim 4.54x +repetition d= 75 n= 149 meas= 5625 shots=500 branch 261.5µs pr204 968.0µs stim 204.2µs branch/stim 1.28x pr204/stim 4.74x +repetition d=225 n= 449 meas= 50625 shots=364 branch 2770.0µs pr204 9810.0µs stim 2894.6µs branch/stim 0.96x pr204/stim 3.39x +repetition d=675 n= 1349 meas= 455625 shots= 19 branch 50380.0µs pr204 121960.0µs stim 53024.1µs branch/stim 0.95x pr204/stim 2.30x +color d= 3 n= 10 meas= 16 shots=500 branch 4.0µs pr204 8.0µs stim 6.5µs branch/stim 0.61x pr204/stim 1.24x +color d= 5 n= 28 meas= 64 shots=500 branch 13.0µs pr204 36.2µs stim 12.8µs branch/stim 1.02x pr204/stim 2.83x +color d= 7 n= 55 meas= 163 shots=500 branch 32.3µs pr204 98.9µs stim 23.8µs branch/stim 1.36x pr204/stim 4.16x +color d= 9 n= 91 meas= 331 shots=500 branch 67.2µs pr204 199.6µs stim 53.8µs branch/stim 1.25x pr204/stim 3.71x +color d= 13 n= 190 meas= 946 shots=500 branch 213.1µs pr204 859.0µs stim 187.1µs branch/stim 1.14x pr204/stim 4.59x +color d= 17 n= 325 meas= 2053 shots=500 branch 582.2µs pr204 2630.0µs stim 512.8µs branch/stim 1.14x pr204/stim 5.13x +color d= 21 n= 496 meas= 3796 shots=500 branch 1090.0µs pr204 3710.0µs stim 1225.4µs branch/stim 0.89x pr204/stim 3.03x +color d= 25 n= 703 meas= 6319 shots=414 branch 2380.0µs pr204 19970.0µs stim 2781.1µs branch/stim 0.86x pr204/stim 7.18x +color d= 31 n= 1081 meas= 11881 shots=178 branch 5620.0µs pr204 59310.0µs stim 7068.0µs branch/stim 0.80x pr204/stim 8.39x +color d= 37 n= 1540 meas= 20008 shots= 77 branch 12560.0µs pr204 146490.0µs stim 15826.9µs branch/stim 0.79x pr204/stim 9.26x +color d= 43 n= 2080 meas= 31186 shots= 42 branch 23670.0µs pr204 344410.0µs stim 32056.8µs branch/stim 0.74x pr204/stim 10.74x diff --git a/bench/stim-compare/sweep.py b/bench/stim-compare/sweep.py new file mode 100644 index 000000000..0d5ceeb87 --- /dev/null +++ b/bench/stim-compare/sweep.py @@ -0,0 +1,114 @@ +"""Size sweep: Stim-generated Clifford memory circuits, timed on ppvm (this +branch and PR 204) and on official Stim, alternating tools over a few rounds. + +Usage: `sweep.py [rounds]`. Writes `circuits/*.stim` and `sweep.csv`. +""" + +import csv +import pathlib +import re +import statistics +import subprocess +import sys + +import stim + +HERE = pathlib.Path(__file__).resolve().parent +WORKTREES = HERE.parents[2] +PYTHON = "/Users/david/git/ppvm/ppvm-python/.venv/bin/python" +TOOLS = { + "branch": [str(HERE / "target/release/stim-compare"), "file"], + "pr204": [str(WORKTREES / "pr204-baseline/bench/stim-compare/target/release/stim-compare"), "file"], + "stim": [PYTHON, str(HERE / "stim_bench.py"), "file"], +} +FAMILIES = [ + ("surface", "surface_code:rotated_memory_z", [3, 5, 7, 9, 11, 15, 19, 23, 27, 31]), + ("repetition", "repetition_code:memory", [3, 9, 25, 75, 225, 675]), + ("color", "color_code:memory_xyz", [3, 5, 7, 9, 13, 17, 21, 25, 31, 37, 43]), +] +CLIFFT_D7 = pathlib.Path.home() / "git/clifft-bench/workloads/circuits/pure_surface_d7_r7_p1e-3.stim" + + +def rewrite_c_xyz(text: str) -> str: + """ppvm-stim has no `C_XYZ`; it equals `H` then `SQRT_X_DAG`.""" + + def expand(m: re.Match) -> str: + indent, targets = m.group(1), m.group(2) + return f"{indent}H{targets}\n{indent}SQRT_X_DAG{targets}" + + return re.sub(r"^(\s*)C_XYZ(\s.*)$", expand, text, flags=re.M) + + +def circuits() -> list[tuple[str, int, pathlib.Path]]: + out_dir = HERE / "circuits" + out_dir.mkdir(exist_ok=True) + out = [("clifft_surface", 7, CLIFFT_D7)] + for family, task, distances in FAMILIES: + for d in distances: + c = stim.Circuit.generated( + task, + distance=d, + rounds=d, + after_clifford_depolarization=1e-3, + before_round_data_depolarization=1e-3, + before_measure_flip_probability=1e-3, + after_reset_flip_probability=1e-3, + ) + path = out_dir / f"{family}_d{d}.stim" + path.write_text(rewrite_c_xyz(str(c)) + "\n") + out.append((family, d, path)) + return out + + +def run(tool: str, path: pathlib.Path, shots: int) -> tuple[float, float, int, int]: + """(execute median in µs, mean ones, qubits, measurements).""" + out = subprocess.run( + [*TOOLS[tool], str(path), str(shots)], capture_output=True, text=True, check=True + ).stdout + head = re.search(r"\[file\] (\d+) qubits, (\d+) measurements", out) + m = re.search(r"execute median\s+([\d.]+)(ns|µs|ms|s)\b", out) + scale = {"ns": 1e-3, "µs": 1.0, "ms": 1e3, "s": 1e6}[m.group(2)] + ones = float(re.search(r"mean ones/shot ([\d.]+)", out).group(1)) + return float(m.group(1)) * scale, ones, int(head.group(1)), int(head.group(2)) + + +def main() -> None: + rounds = int(sys.argv[1]) if len(sys.argv) > 1 else 3 + rows = [] + for family, d, path in circuits(): + pilot, *_ = run("branch", path, 3) + shots = max(10, min(500, int(1.0e6 / max(pilot, 1.0)))) + times = {tool: [] for tool in TOOLS} + ones = {} + for _ in range(rounds): + for tool in TOOLS: + t, o, n_qubits, n_meas = run(tool, path, shots) + times[tool].append(t) + ones[tool] = o + med = {tool: statistics.median(ts) for tool, ts in times.items()} + row = { + "family": family, + "d": d, + "qubits": n_qubits, + "measurements": n_meas, + "shots": shots, + **{f"{tool}_us": round(med[tool], 1) for tool in TOOLS}, + "branch_over_stim": round(med["branch"] / med["stim"], 3), + "pr204_over_stim": round(med["pr204"] / med["stim"], 3), + **{f"{tool}_ones": ones[tool] for tool in TOOLS}, + } + rows.append(row) + print( + f"{family:14s} d={d:3d} n={n_qubits:5d} meas={n_meas:7d} shots={shots:3d} " + f"branch {med['branch']:10.1f}µs pr204 {med['pr204']:10.1f}µs stim {med['stim']:10.1f}µs " + f"branch/stim {row['branch_over_stim']:.2f}x pr204/stim {row['pr204_over_stim']:.2f}x", + flush=True, + ) + with open(HERE / "sweep.csv", "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + +if __name__ == "__main__": + main() diff --git a/crates/ppvm-stim/src/executor/adapter/mod.rs b/crates/ppvm-stim/src/executor/adapter/mod.rs index 16647c073..5d81ed2aa 100644 --- a/crates/ppvm-stim/src/executor/adapter/mod.rs +++ b/crates/ppvm-stim/src/executor/adapter/mod.rs @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2026 The PPVM Authors // SPDX-License-Identifier: Apache-2.0 +use super::helpers::measure_reset_z; + mod sealed { pub trait Sealed {} @@ -113,6 +115,81 @@ pub trait StimTableau: sealed::Sealed { rng: &mut R, ) -> Option; fn flip_with_prob(&mut self, bit: bool, p: f64, rng: &mut R) -> bool; + + /// `R` on every target, in order. Backends may batch the measurements. + fn reset_many(&mut self, q: &[usize], rng: &mut R) + where + Self: Sized, + { + for &q in q { + self.reset(q, rng); + } + } + + /// `DEPOLARIZE1(p)` on every target. Backends may sample the errors jointly. + fn depolarize1_many(&mut self, q: &[usize], p: f64, rng: &mut R) + where + Self: Sized, + { + for &q in q { + self.depolarize1(q, p, rng); + } + } + + /// `DEPOLARIZE2(p)` on consecutive target pairs. Backends may sample the + /// errors jointly. + fn depolarize2_many(&mut self, q: &[usize], p: f64, rng: &mut R) + where + Self: Sized, + { + for &[a, b] in q.as_chunks::<2>().0 { + self.depolarize2(a, b, p, rng); + } + } + + /// `X`/`Y`/`Z` with probabilities `p` on every target (`X_ERROR`, `Y_ERROR`, + /// `Z_ERROR`, `PAULI_CHANNEL_1`). Backends may sample the errors jointly. + fn pauli_error_many(&mut self, q: &[usize], p: [f64; 3], rng: &mut R) + where + Self: Sized, + { + for &q in q { + self.pauli_error(q, p, rng); + } + } + + /// `M(noise)` on every target, in order, pushing each recorded bit onto + /// `results`. Backends may batch the measurements. + fn measure_noisy_many( + &mut self, + q: &[usize], + noise: f64, + rng: &mut R, + results: &mut Vec>, + ) where + Self: Sized, + { + for &q in q { + results.push(self.measure_noisy(q, noise, rng)); + } + } + + /// `MR(noise)` on every target, in order, pushing each recorded bit onto + /// `results`. Backends may batch the measurements. + fn measure_reset_many( + &mut self, + q: &[usize], + noise: f64, + rng: &mut R, + results: &mut Vec>, + ) where + Self: Sized, + { + for &q in q { + results.push(measure_reset_z(self, q, noise, rng)); + } + } + fn measurement_record(&self) -> &[Option]; fn append_measurement_record(&mut self, result: Option); fn overwrite_last_measurement_record(&mut self, result: Option); diff --git a/crates/ppvm-stim/src/executor/adapter/traits_2.rs b/crates/ppvm-stim/src/executor/adapter/traits_2.rs index 718dcbf6c..38d1832be 100644 --- a/crates/ppvm-stim/src/executor/adapter/traits_2.rs +++ b/crates/ppvm-stim/src/executor/adapter/traits_2.rs @@ -8,6 +8,8 @@ use ppvm_tableau_2::prelude::{ }; use super::StimTableau; +use crate::executor::helpers::{has_repeats, measure_reset_z}; +use rand::RngExt; macro_rules! unary { ($name:ident, $trait:ident) => { @@ -31,6 +33,77 @@ macro_rules! batch { }; } +/// Flip each of the last `outcomes.len()` records with probability `noise`, +/// pushing the recorded bits onto `results`. A lost qubit's `None` stays. +fn record_with_noise( + tab: &mut GeneralizedTableau, + outcomes: &[Option], + noise: f64, + rng: &mut R, + results: &mut Vec>, +) { + let base = tab.measurement_record.len() - outcomes.len(); + for (k, outcome) in outcomes.iter().enumerate() { + let recorded = outcome.map(|b| GeneralizedTableau::::flip_with_prob(b, noise, rng)); + tab.measurement_record[base + k] = recorded; + results.push(recorded); + } +} + +/// Call `hit(i, rng)` for each index in `0..n` that an independent event of +/// probability `p` lands on, in order. Gaps between hits are geometric, so this +/// draws once per hit rather than once per index — Stim's `RareErrorIterator`. +fn for_each_hit( + n: usize, + p: f64, + rng: &mut R, + mut hit: impl FnMut(usize, &mut R), +) { + if p <= 0.0 { + return; + } + if p >= 1.0 { + (0..n).for_each(|i| hit(i, rng)); + return; + } + let log_miss = (-p).ln_1p(); + let mut i = 0; + while i < n { + // `1 - u` is in (0, 1], so the gap is finite and non-negative. + let gap = (1.0 - rng.random::()).ln() / log_miss; + if gap >= (n - i) as f64 { + return; + } + i += gap as usize; + hit(i, rng); + i += 1; + } +} + +/// Pauli `1 = X`, `2 = Y`, `3 = Z` on `q`; `0` is the identity. A lost qubit is +/// left alone by the gates themselves. +fn apply_pauli(tab: &mut GeneralizedTableau, q: usize, pauli: usize) { + match pauli { + 1 => Clifford::x(tab, q), + 2 => Clifford::y(tab, q), + 3 => Clifford::z(tab, q), + _ => {} + } +} + +/// `X` on every target whose true outcome was `1`. +fn reset_ones( + tab: &mut GeneralizedTableau, + q: &[usize], + outcomes: &[Option], +) { + for (&q, &outcome) in q.iter().zip(outcomes) { + if outcome == Some(true) { + Clifford::x(tab, q); + } + } +} + impl StimTableau for GeneralizedTableau where I: Bitstring, @@ -124,7 +197,7 @@ where q: &[usize], rng: &mut R, ) -> Vec> { - Measure::measure_many(self, q, rng) + GeneralizedTableau::measure_batch(self, q, rng) } fn measure_noisy( &mut self, @@ -137,6 +210,82 @@ where fn flip_with_prob(&mut self, bit: bool, p: f64, rng: &mut R) -> bool { GeneralizedTableau::::flip_with_prob(bit, p, rng) } + // Batch the measurements, then apply the `X` resets: `X_q` commutes with + // `Z_p` for `p != q`, so deferring it is exact. A repeated target can't defer. + fn reset_many(&mut self, q: &[usize], rng: &mut R) { + if has_repeats(q) { + q.iter().for_each(|&q| Reset::reset(self, q, rng)); + return; + } + let outcomes = self.measure_batch(q, rng); + let kept = self.measurement_record.len() - q.len(); + self.measurement_record.truncate(kept); + reset_ones(self, q, &outcomes); + } + fn measure_reset_many( + &mut self, + q: &[usize], + noise: f64, + rng: &mut R, + results: &mut Vec>, + ) { + if has_repeats(q) { + q.iter() + .for_each(|&q| results.push(measure_reset_z(self, q, noise, rng))); + return; + } + let outcomes = self.measure_batch(q, rng); + record_with_noise(self, &outcomes, noise, rng, results); + reset_ones(self, q, &outcomes); + } + // A repeated target needs no fallback: its second measurement is simply + // deterministic, and each record's flip is independent. + fn measure_noisy_many( + &mut self, + q: &[usize], + noise: f64, + rng: &mut R, + results: &mut Vec>, + ) { + let outcomes = self.measure_batch(q, rng); + record_with_noise(self, &outcomes, noise, rng, results); + } + + // Noise draws once per error (geometric gaps), not once per target; each + // error then picks its Pauli. Same distribution as the per-target channels. + fn depolarize1_many(&mut self, q: &[usize], p: f64, rng: &mut R) { + for_each_hit(q.len(), p, rng, |i, rng| { + apply_pauli(self, q[i], rng.random_range(1..4)); + }); + } + fn depolarize2_many(&mut self, q: &[usize], p: f64, rng: &mut R) { + let pairs = q.len() / 2; + for_each_hit(pairs, p, rng, |i, rng| { + let (a, b) = (q[2 * i], q[2 * i + 1]); + if self.is_lost[a] || self.is_lost[b] { + return; + } + // One of the 15 non-identity pairs `(k / 4, k % 4)`, as `depolarize2`. + let k = rng.random_range(1..16); + apply_pauli(self, a, k / 4); + apply_pauli(self, b, k % 4); + }); + } + fn pauli_error_many(&mut self, q: &[usize], p: [f64; 3], rng: &mut R) { + let total: f64 = p.iter().sum(); + for_each_hit(q.len(), total, rng, |i, rng| { + let r = rng.random::() * total; + let pauli = if r < p[0] { + 1 + } else if r < p[0] + p[1] { + 2 + } else { + 3 + }; + apply_pauli(self, q[i], pauli); + }); + } + fn measurement_record(&self) -> &[Option] { self.current_measurement_record() } diff --git a/crates/ppvm-stim/src/executor/gates.rs b/crates/ppvm-stim/src/executor/gates.rs index 946d181c8..e02146bde 100644 --- a/crates/ppvm-stim/src/executor/gates.rs +++ b/crates/ppvm-stim/src/executor/gates.rs @@ -13,11 +13,7 @@ pub(super) fn execute( ) { let GateOp { name, targets, .. } = op; match name { - GateName::Reset | GateName::ResetZ => { - for &target in targets { - tab.reset(qubit(target), rng); - } - } + GateName::Reset | GateName::ResetZ => tab.reset_many(&qubits(targets), rng), GateName::ResetX => { for &target in targets { let q = qubit(target); diff --git a/crates/ppvm-stim/src/executor/helpers.rs b/crates/ppvm-stim/src/executor/helpers.rs index 6707be6d3..fcef6ddee 100644 --- a/crates/ppvm-stim/src/executor/helpers.rs +++ b/crates/ppvm-stim/src/executor/helpers.rs @@ -38,6 +38,13 @@ pub(super) fn record_bit(record: &[Option], k: usize) -> bool { .unwrap_or(false) } +/// Whether any qubit appears more than once in `qs`. +pub(super) fn has_repeats(qs: &[usize]) -> bool { + let mut sorted: SmallVec<[usize; TARGETS_INLINE]> = qs.into(); + sorted.sort_unstable(); + sorted.windows(2).any(|w| w[0] == w[1]) +} + pub(super) fn measure_reset_z( tab: &mut T, q: usize, diff --git a/crates/ppvm-stim/src/executor/measure.rs b/crates/ppvm-stim/src/executor/measure.rs index 37d77456d..1279ee4b0 100644 --- a/crates/ppvm-stim/src/executor/measure.rs +++ b/crates/ppvm-stim/src/executor/measure.rs @@ -4,7 +4,7 @@ use stim_parser::prelude::{MeasureName, MeasureOp, MppOp, PauliAxis}; use super::StimTableau; -use super::helpers::measure_reset_z; +use super::helpers::has_repeats; pub(super) fn execute( op: &MeasureOp, @@ -22,49 +22,31 @@ pub(super) fn execute( match name { MeasureName::M | MeasureName::MZ => { if noise > 0.0 { - for &q in targets { - results.push(tab.measure_noisy(q, noise, rng)); - } + tab.measure_noisy_many(targets, noise, rng, results); } else { results.extend(tab.measure_many(targets, rng)); } } - MeasureName::MR => { - for &q in targets { - results.push(measure_reset_z(tab, q, noise, rng)); - } - } - MeasureName::MX => { - for &q in targets { - tab.h(q); - results.push(tab.measure_noisy(q, noise, rng)); - tab.h(q); - } - } - MeasureName::MY => { - for &q in targets { - tab.s_dag(q); - tab.h(q); - results.push(tab.measure_noisy(q, noise, rng)); - tab.h(q); - tab.s(q); - } - } - MeasureName::MRX => { - for &q in targets { - tab.h(q); - results.push(measure_reset_z(tab, q, noise, rng)); - tab.h(q); - } + MeasureName::MR => tab.measure_reset_many(targets, noise, rng, results), + MeasureName::MX | MeasureName::MY => { + let axis = if *name == MeasureName::MX { + PauliAxis::X + } else { + PauliAxis::Y + }; + measure_in_basis(tab, targets, axis, results, |tab, q, results| { + tab.measure_noisy_many(q, noise, rng, results) + }); } - MeasureName::MRY => { - for &q in targets { - tab.s_dag(q); - tab.h(q); - results.push(measure_reset_z(tab, q, noise, rng)); - tab.h(q); - tab.s(q); - } + MeasureName::MRX | MeasureName::MRY => { + let axis = if *name == MeasureName::MRX { + PauliAxis::X + } else { + PauliAxis::Y + }; + measure_in_basis(tab, targets, axis, results, |tab, q, results| { + tab.measure_reset_many(q, noise, rng, results) + }); } MeasureName::MXX | MeasureName::MYY | MeasureName::MZZ | MeasureName::MPP => { unreachable!("unsupported measure {name:?} should have been rejected by validate") @@ -97,6 +79,29 @@ pub(super) fn execute_mpp( } } +/// Measure `targets` along `axis` with `measure`, a Z-basis batch. Distinct +/// targets rotate onto Z all at once around one batch; a repeated target keeps +/// the per-target order, where the rotations between its measurements matter. +fn measure_in_basis( + tab: &mut T, + targets: &[usize], + axis: PauliAxis, + results: &mut Vec>, + mut measure: impl FnMut(&mut T, &[usize], &mut Vec>), +) { + if has_repeats(targets) { + for &q in targets { + basis_to_z(tab, axis, q); + measure(tab, &[q], results); + basis_from_z(tab, axis, q); + } + return; + } + targets.iter().for_each(|&q| basis_to_z(tab, axis, q)); + measure(tab, targets, results); + targets.iter().for_each(|&q| basis_from_z(tab, axis, q)); +} + fn basis_to_z(tab: &mut T, axis: PauliAxis, q: usize) { match axis { PauliAxis::X => tab.h(q), diff --git a/crates/ppvm-stim/src/executor/noise.rs b/crates/ppvm-stim/src/executor/noise.rs index 35fe64efb..21a5d6221 100644 --- a/crates/ppvm-stim/src/executor/noise.rs +++ b/crates/ppvm-stim/src/executor/noise.rs @@ -19,22 +19,9 @@ pub(super) fn execute( .. } = op; match name { - NoiseName::Depolarize1 => { - for &q in targets { - tab.depolarize1(q, args[0], rng); - } - } - NoiseName::Depolarize2 => { - for (a, b) in targets.iter().copied().tuples() { - tab.depolarize2(a, b, args[0], rng); - } - } - NoiseName::PauliChannel1 => { - let p = [args[0], args[1], args[2]]; - for &q in targets { - tab.pauli_error(q, p, rng); - } - } + NoiseName::Depolarize1 => tab.depolarize1_many(targets, args[0], rng), + NoiseName::Depolarize2 => tab.depolarize2_many(targets, args[0], rng), + NoiseName::PauliChannel1 => tab.pauli_error_many(targets, [args[0], args[1], args[2]], rng), NoiseName::PauliChannel2 => { debug_assert!(targets.len().is_even()); let p = std::array::from_fn(|i| args[i]); @@ -50,9 +37,7 @@ pub(super) fn execute( NoiseName::ZError => [zero, zero, args[0]], _ => unreachable!(), }; - for &q in targets { - tab.pauli_error(q, p, rng); - } + tab.pauli_error_many(targets, p, rng); } NoiseName::IError | NoiseName::HeraldedErase diff --git a/crates/ppvm-stim/tests/executor.rs b/crates/ppvm-stim/tests/executor.rs index 31cbd3ee5..7e59c598d 100644 --- a/crates/ppvm-stim/tests/executor.rs +++ b/crates/ppvm-stim/tests/executor.rs @@ -592,3 +592,124 @@ fn sample_keeps_non_sync_factories_on_serial_builds() { assert_eq!(shots.len(), 3); assert_eq!(calls.get(), 3); } + +/// Batched `MR`/`R` must reset every target, keep the record in target order, +/// and see each earlier reset of a repeated target. +#[test] +fn batched_measure_reset_and_reset_over_many_targets() { + // Bell pair on (0, 1), |1> on 2, |+> on 3; the random targets follow + // deterministic ones in the target list. + let src = "H 0\nCX 0 1\nX 2\nH 3\nMR 2 1 0 3\nM 0 1 2 3"; + for seed in 0..16 { + let r = run_seeded(src, 4, seed); + assert_eq!(r[0], Some(true)); + assert_eq!(r[1], r[2], "the Bell pair must agree"); + assert_eq!(&r[4..], &[Some(false); 4], "MR must reset every target"); + } + + let (results, _) = run("X 0\nX 1\nR 1 0\nM 0 1", 2); + assert_eq!(results, vec![Some(false), Some(false)]); +} + +#[test] +fn batched_measure_reset_repeated_target_sees_the_reset() { + for seed in 0..16 { + let r = run_seeded("H 0\nMR 0 0", 1, seed); + assert_eq!(r[1], Some(false), "the second MR of 0 follows its reset"); + } + let (results, _) = run("X 0\nMR(1.0) 0 1", 2); + assert_eq!( + results, + vec![Some(false), Some(true)], + "MR(1) flips each record" + ); +} + +/// Batched `MX`/`MY`/`MRX`/`MRY` and noisy `M` must keep per-target semantics: +/// the right basis, the record order, the resets and the per-record noise. +#[test] +fn batched_basis_and_noisy_measurements() { + let (results, _) = run("H 0\nH 1\nH 2\nS 2\nMX 1 0\nMY 2", 3); + assert_eq!( + results, + vec![Some(false); 3], + "|+> in X and |+i> in Y are 0" + ); + + for seed in 0..16 { + // X⊗X = +1 on the Bell pair; the |+> qubit 2 sits between them. + let r = run_seeded("H 0\nCX 0 1\nH 2\nMX 0 2 1\nMRX 0 1\nMX 0 1", 3, seed); + assert_eq!(r[0], r[2], "the Bell pair agrees in X"); + assert_eq!(r[1], Some(false)); + assert_eq!(r[3], r[4], "MRX measures the same collapsed pair"); + assert_eq!(&r[5..], &[Some(false); 2], "MRX resets to |+>"); + + let r = run_seeded("MRX 0 0\nMRY 1 1", 2, seed); + assert_eq!( + (r[1], r[3]), + (Some(false), Some(false)), + "a repeat sees the reset" + ); + } + + let (results, _) = run("X 0\nM(1.0) 0 1\nH 2\nMX(1.0) 2", 3); + assert_eq!( + results, + vec![Some(false), Some(true), Some(true)], + "p = 1 flips each record" + ); +} + +/// Fraction of `1`s per measured qubit over `shots` seeded runs of `src`, and +/// the fraction of shots where qubits 0 and 1 both read `1`. +fn flip_rates(src: &str, shots: u64) -> (f64, f64) { + let (mut ones, mut both, mut total) = (0usize, 0usize, 0usize); + for seed in 0..shots { + let r = run_seeded(src, 8, 1000 + seed); + ones += r.iter().filter(|&&b| b == Some(true)).count(); + total += r.len(); + both += usize::from(r[0] == Some(true) && r[1] == Some(true)); + } + (ones as f64 / total as f64, both as f64 / shots as f64) +} + +/// The batched noise channels must keep the per-target distributions: each +/// target independently, each Pauli with its probability, and `DEPOLARIZE2`'s +/// 15 pair errors correlated across the pair. +#[test] +fn noise_channels_keep_their_distributions() { + let all = "0 1 2 3 4 5 6 7"; + let cases = [ + (format!("X_ERROR(0.1) {all}\nM {all}"), 0.1), + (format!("Y_ERROR(0.1) {all}\nM {all}"), 0.1), + ( + format!("H {all}\nZ_ERROR(0.1) {all}\nH {all}\nM {all}"), + 0.1, + ), + // X and Y flip a Z measurement: 2/3 of p. + (format!("DEPOLARIZE1(0.3) {all}\nM {all}"), 0.2), + ( + format!("PAULI_CHANNEL_1(0.05, 0.1, 0.2) {all}\nM {all}"), + 0.15, + ), + // 8 of the 15 pair errors put X or Y on a given qubit. + (format!("DEPOLARIZE2(0.3) {all}\nM {all}"), 0.16), + ]; + for (src, expected) in &cases { + let (rate, _) = flip_rates(src, 600); + assert!( + (rate - expected).abs() < 0.012, + "{src}: {rate} vs {expected}" + ); + } + + // 4 of the 15 pair errors flip both qubits of a pair. + let (_, both) = flip_rates(&format!("DEPOLARIZE2(0.3) {all}\nM {all}"), 4000); + assert!( + (both - 0.08).abs() < 0.012, + "DEPOLARIZE2 pair correlation: {both}" + ); + + assert_eq!(flip_rates(&format!("X_ERROR(1) {all}\nM {all}"), 4).0, 1.0); + assert_eq!(flip_rates(&format!("X_ERROR(0) {all}\nM {all}"), 4).0, 0.0); +} diff --git a/crates/ppvm-tableau-2/src/clifford.rs b/crates/ppvm-tableau-2/src/clifford.rs index 16f71bedb..1e49bda32 100644 --- a/crates/ppvm-tableau-2/src/clifford.rs +++ b/crates/ppvm-tableau-2/src/clifford.rs @@ -355,15 +355,10 @@ impl Clifford for Tableau { #[inline] fn cnot(&mut self, control: usize, target: usize) { - sweep2!( - self, - control, - target, - inverse: prepend_cnot, - |xc: &mut [u64], zc: &mut [u64], xt: &mut [u64], zt: &mut [u64], ph: &mut [u64]| { - blocks::cnot(xc, zc, xt, zt, ph) - } - ); + self.invalidate_hash(); + // One pass for the bits, the forward phases and the inverse-sign products. + let (g_x, g_z) = self.data.cnot_fused(control, target); + self.prepend_cnot(control, target, g_x, g_z); } #[inline] diff --git a/crates/ppvm-tableau-2/src/data.rs b/crates/ppvm-tableau-2/src/data.rs index cffcb6d21..6225ffd6c 100644 --- a/crates/ppvm-tableau-2/src/data.rs +++ b/crates/ppvm-tableau-2/src/data.rs @@ -30,6 +30,9 @@ use ppvm_traits_2::{Indexable, Pauli, Scale, Support}; use crate::storage::{BITS_PER_WORD, HALVES, Half, Orientation, Plane, TableauData, blocks}; +/// One bit column per half, indexed by [`Half`]. +pub(crate) type HalfColumns = [Vec; 2]; + /// Bit-string index type addressing one branch of the amplitude vector. /// /// Blanket-implemented for every primitive (and `bnum`-style) unsigned integer @@ -441,6 +444,29 @@ impl Tableau { /// one of the two columns at `addr0`, or their `XOR` for `Y` — contiguous in /// the canonical orientation, where the replaced code probed the same site /// in each of `n` separately addressed rows. + /// [`Self::anticommutation_column`] of both halves for `Z_addr0`, indexed by + /// [`Half`]. A measurement reads these three times, so it gathers them once. + pub(crate) fn z_anticommutation_columns(&self, addr0: usize) -> HalfColumns { + let mut columns = HalfColumns::default(); + self.z_anticommutation_columns_into(addr0, &mut columns); + columns + } + + /// The stabilizer half of [`Self::z_anticommutation_columns`], into `column`. + pub(crate) fn z_anticommuting_stabilizers_into(&self, addr0: usize, column: &mut Vec) { + column.resize(self.data.stride(), 0); + self.data.gather_column(Half::Stab, Plane::X, addr0, column); + } + + /// [`Self::z_anticommutation_columns`] into reused buffers. + pub(crate) fn z_anticommutation_columns_into(&self, addr0: usize, columns: &mut HalfColumns) { + for half in HALVES { + let column = &mut columns[half as usize]; + column.resize(self.data.stride(), 0); + self.data.gather_column(half, Plane::X, addr0, column); + } + } + pub(crate) fn anticommutation_column( &self, half: Half, @@ -544,6 +570,19 @@ impl Tableau { addr0: usize, q_idx: usize, outcome: bool, + ) { + let columns = self.z_anticommutation_columns(addr0); + self.update_tableau_with_columns(addr0, q_idx, outcome, &columns); + } + + /// [`Self::update_tableau_according_to_outcome`] with `addr0`'s + /// [`Self::z_anticommutation_columns`] already gathered. + pub(crate) fn update_tableau_with_columns( + &mut self, + addr0: usize, + q_idx: usize, + outcome: bool, + columns: &HalfColumns, ) { let n = self.n_qubits(); let stride = self.data.stride(); @@ -557,9 +596,9 @@ impl Tableau { // Both paths are exercised by the conformance differentials. if self.data.orientation() == Orientation::RowMajor { if self.data.inverse_valid() { - self.project_inverse(addr0, q_idx, outcome); + self.project_inverse(addr0, q_idx, outcome, columns); } - self.project_row_major(addr0, q_idx, outcome); + self.project_row_major(addr0, q_idx, outcome, columns); return; } @@ -574,11 +613,11 @@ impl Tableau { let selected = blocks::count_set(self.data.major(Half::Stab, Plane::X, addr0)) + blocks::count_set(self.data.major(Half::Destab, Plane::X, addr0)); if blocks::prefer_gather(selected, n) { - self.project_inverse(addr0, q_idx, outcome); + self.project_inverse(addr0, q_idx, outcome, columns); } else { self.enter_row_major(); - self.project_inverse(addr0, q_idx, outcome); - self.project_row_major(addr0, q_idx, outcome); + self.project_inverse(addr0, q_idx, outcome, columns); + self.project_row_major(addr0, q_idx, outcome, columns); self.exit_row_major(); return; } @@ -588,10 +627,8 @@ impl Tableau { // x_g[addr0]`, so each selector is one contiguous column; `q_idx` is // cleared because the pivot is not multiplied into itself — the // replaced loop's `if i == q_idx { continue }`. - let mut sel = [vec![0u64; stride], vec![0u64; stride]]; + let mut sel = columns.clone(); for half in HALVES { - self.data - .gather_column(half, Plane::X, addr0, &mut sel[half as usize]); TableauData::set_bit(&mut sel[half as usize], q_idx, false); } @@ -669,15 +706,79 @@ impl Tableau { /// generators, versus the column-wise form's fixed sweep over all `n` qubit /// columns. Cheaper exactly when the frame is dense, which is when a caller /// bothered to take the guard. - fn project_row_major(&mut self, addr0: usize, q_idx: usize, outcome: bool) { + /// Collapse onto outcome `outcome` of `Z_addr0` the way Stim's + /// `collapse_qubit_z` does: eliminate over the stabilizers only, so the + /// destabilizer column is never needed. `stab_column` selects the + /// stabilizers anticommuting with `Z_addr0`; `pivot` is one of them. + /// + /// As appends `U ↦ U·V`: `CX(p, k)` per selected `k` (`s_k ← s_k·s_p`, + /// `d_p ← d_p·d_k`), `S(p)` if `d_p` then anticommutes with `Z_addr0` + /// (`d_p ← i·d_p·s_p`), `H(p)` (swap the pair), and `X(p)` (negate `s_p`) + /// if the outcome needs it. The new `s_p` is `±Z_addr0` times other + /// stabilizers rather than `±Z_addr0` itself. Row-major only, inverse valid. + pub(crate) fn collapse_z( + &mut self, + addr0: usize, + pivot: usize, + outcome: bool, + stab_column: &[u64], + ) { + debug_assert_eq!(self.data.orientation(), Orientation::RowMajor); + self.invalidate_hash(); + let flip = self.collapse_inverse(addr0, pivot, outcome, stab_column, None); + let n = self.n_qubits(); let stride = self.data.stride(); - let mut stab_selector = vec![0u64; stride]; - let mut destab_selector = vec![0u64; stride]; - self.data - .gather_column(Half::Stab, Plane::X, addr0, &mut stab_selector); - self.data - .gather_column(Half::Destab, Plane::X, addr0, &mut destab_selector); + let data = &mut self.data; + let mut s_p = ScratchRow::zeroed(stride); + s_p.x + .copy_from_slice(data.major(Half::Stab, Plane::X, pivot)); + s_p.z + .copy_from_slice(data.major(Half::Stab, Plane::Z, pivot)); + s_p.phase = data.phase_of(Half::Stab, pivot); + let mut d_p = ScratchRow::zeroed(stride); + d_p.x + .copy_from_slice(data.major(Half::Destab, Plane::X, pivot)); + d_p.z + .copy_from_slice(data.major(Half::Destab, Plane::Z, pivot)); + d_p.phase = data.phase_of(Half::Destab, pivot); + + let mut src = ScratchRow::zeroed(stride); + for k in (0..n).filter(|&k| k != pivot && TableauData::bit(stab_column, k)) { + data.multiply_row_by(Half::Stab, k, &s_p.x, &s_p.z, s_p.phase); + d_p.mul_generator(data, Half::Destab, k, &mut src); + } + if TableauData::bit(&d_p.x, addr0) { + let g = blocks::row_multiply(&mut d_p.x, &mut d_p.z, &s_p.x, &s_p.z); + d_p.add_phase(g + s_p.phase + 1); + } + + data.major_mut(Half::Destab, Plane::X, pivot) + .copy_from_slice(&s_p.x); + data.major_mut(Half::Destab, Plane::Z, pivot) + .copy_from_slice(&s_p.z); + data.set_phase_of(Half::Destab, pivot, s_p.phase); + data.major_mut(Half::Stab, Plane::X, pivot) + .copy_from_slice(&d_p.x); + data.major_mut(Half::Stab, Plane::Z, pivot) + .copy_from_slice(&d_p.z); + data.set_phase_of( + Half::Stab, + pivot, + (d_p.phase + if flip { 2 } else { 0 }) % 4, + ); + } + + fn project_row_major( + &mut self, + addr0: usize, + q_idx: usize, + outcome: bool, + columns: &HalfColumns, + ) { + let n = self.n_qubits(); + let stride = self.data.stride(); + let [destab_selector, stab_selector] = columns; let data = &mut self.data; @@ -696,10 +797,10 @@ impl Tableau { if i == q_idx { continue; } - if TableauData::bit(&stab_selector, i) { + if TableauData::bit(stab_selector, i) { data.multiply_row_by(Half::Stab, i, &pivot.x, &pivot.z, pivot.phase); } - if TableauData::bit(&destab_selector, i) { + if TableauData::bit(destab_selector, i) { data.multiply_row_by(Half::Destab, i, &pivot.x, &pivot.z, pivot.phase); } } @@ -1286,6 +1387,36 @@ impl GeneralizedTableau { /// branch takes; nothing logically mutates, and the frame is byte-identical /// on return. pub fn compute_decomposition(&mut self, addr0: usize, pauli: Pauli) -> (u8, I, I) { + let anticomm = HALVES.map(|half| self.tableau.anticommutation_column(half, addr0, pauli)); + self.compute_decomposition_with(addr0, pauli, &anticomm) + } + + /// [`Self::compute_decomposition`] with both halves' anticommutation + /// columns (`Tableau::anticommutation_column`) already gathered. + pub(crate) fn compute_decomposition_with( + &mut self, + addr0: usize, + pauli: Pauli, + anticomm: &HalfColumns, + ) -> (u8, I, I) { + let n = self.n_qubits(); + let phase = self.decomposition_phase_with(addr0, pauli, anticomm); + let [destab_anticomm, stab_anticomm] = anticomm; + ( + phase, + bits_to_index::(stab_anticomm, n), + bits_to_index::(destab_anticomm, n), + ) + } + + /// The phase of [`Self::compute_decomposition_with`], without widening the + /// masks into branch indices. + pub(crate) fn decomposition_phase_with( + &mut self, + addr0: usize, + pauli: Pauli, + anticomm: &HalfColumns, + ) -> u8 { debug_assert_ne!(pauli, Pauli::I); let n = self.n_qubits(); let stride = self.tableau.data.stride(); @@ -1295,14 +1426,7 @@ impl GeneralizedTableau { // replaced code probed the same site on all `2n` separately addressed // rows; the *values*, the visit order and the accumulated phase below // are unchanged. - let destab_anticomm = self - .tableau - .anticommutation_column(Half::Destab, addr0, pauli); - let stab_anticomm = self - .tableau - .anticommutation_column(Half::Stab, addr0, pauli); - let destab_anticomm_bits = bits_to_index::(&destab_anticomm, n); - let stab_anticomm_bits = bits_to_index::(&stab_anticomm, n); + let [destab_anticomm, stab_anticomm] = anticomm; // The visit order — all selected stabilizers ascending, then all // selected destabilizers ascending — is a genuine convention, not a free @@ -1313,7 +1437,7 @@ impl GeneralizedTableau { p_word.set(addr0, pauli); let mut src = ScratchRow::zeroed(stride); for i in 0..n { - if TableauData::bit(&destab_anticomm, i) { + if TableauData::bit(destab_anticomm, i) { // The stabilizer is its own inverse up to its phase; rather // than inverting we multiply and divide out the phase // squared. @@ -1323,7 +1447,7 @@ impl GeneralizedTableau { } } for i in 0..n { - if TableauData::bit(&stab_anticomm, i) { + if TableauData::bit(stab_anticomm, i) { let phase = data.phase_of(Half::Destab, i); p_word.mul_generator(data, Half::Destab, i, &mut src); p_word.add_phase(8 - 2 * phase); @@ -1340,20 +1464,18 @@ impl GeneralizedTableau { if self.tableau.inverse_readable(pauli) { let phase = self.tableau - .decomposition_phase(addr0, pauli, &destab_anticomm, &stab_anticomm); + .decomposition_phase(addr0, pauli, destab_anticomm, stab_anticomm); debug_assert_eq!(phase, fold(&self.tableau.data)); - return (phase, stab_anticomm_bits, destab_anticomm_bits); + return phase; } - let selected = blocks::count_set(&destab_anticomm) + blocks::count_set(&stab_anticomm); - let phase = if blocks::prefer_gather(selected, n) { + let selected = blocks::count_set(destab_anticomm) + blocks::count_set(stab_anticomm); + if blocks::prefer_gather(selected, n) { fold(&self.tableau.data) } else { let guard = TransposedTableau::new(&mut self.tableau); fold(guard.data()) - }; - - (phase, stab_anticomm_bits, destab_anticomm_bits) + } } /// Multi-qubit generalization of [`Self::compute_decomposition`]: conjugate diff --git a/crates/ppvm-tableau-2/src/inverse.rs b/crates/ppvm-tableau-2/src/inverse.rs index 6b91ba7d6..4cc07cd06 100644 --- a/crates/ppvm-tableau-2/src/inverse.rs +++ b/crates/ppvm-tableau-2/src/inverse.rs @@ -63,8 +63,8 @@ use ppvm_traits_2::Pauli; -use crate::data::Tableau; -use crate::storage::{HALVES, Half, InvRow, Orientation, Plane, TableauData, blocks}; +use crate::data::{HalfColumns, Tableau}; +use crate::storage::{Half, InvRow, Orientation, TableauData, blocks}; impl Tableau { /// Whether the inverse-row signs can be read. @@ -233,16 +233,17 @@ impl Tableau { // ─── Two-qubit rules ────────────────────────────────────────────────── /// `CNOT†`: `X_c ↦ X_cX_t`, `Z_t ↦ Z_cZ_t`, the other two fixed. - pub(crate) fn prepend_cnot(&mut self, control: usize, target: usize) { + /// + /// `g_x` / `g_z` are the `g`-rule terms of the products `ix_c·ix_t` and + /// `iz_c·iz_t`, which [`TableauData::cnot_fused`] reads in the same pass as + /// the forward update. + pub(crate) fn prepend_cnot(&mut self, control: usize, target: usize, g_x: u8, g_z: u8) { if !self.inverse_valid() { return; } - let nx = self - .data - .inv_pair_phase((InvRow::X, control), (InvRow::X, target)); - let nz = self - .data - .inv_pair_phase((InvRow::Z, control), (InvRow::Z, target)); + let sign = |row| self.data.inv_sign(row, control) + self.data.inv_sign(row, target); + let nx = (sign(InvRow::X) + g_x) % 4; + let nz = (sign(InvRow::Z) + g_z) % 4; self.data.set_inv_sign(InvRow::X, control, nx); self.data.set_inv_sign(InvRow::Z, target, nz); } @@ -363,7 +364,29 @@ impl Tableau { /// Pauli's *bits* alone, so the frame's `ℤ/4` phases — which the forward /// projection is busy folding `g`-rules into — never enter here. The two /// bookkeepings are independent computations of the same frame. - pub(crate) fn project_inverse(&mut self, addr0: usize, pivot: usize, outcome: bool) { + pub(crate) fn project_inverse( + &mut self, + addr0: usize, + pivot: usize, + outcome: bool, + columns: &HalfColumns, + ) { + let [destab, stab] = columns; + self.collapse_inverse(addr0, pivot, outcome, stab, Some(destab)); + } + + /// The appends behind [`Self::project_inverse`]. With `destab_column` `None` + /// the destabilizer `CZ`s are skipped — Stim's `collapse_qubit_z`, which + /// leaves `s_p` a product of `±Z_a` and other stabilizers. Returns whether + /// the final `append_X(p)` ran. + pub(crate) fn collapse_inverse( + &mut self, + addr0: usize, + pivot: usize, + outcome: bool, + stab_column: &[u64], + destab_column: Option<&[u64]>, + ) -> bool { debug_assert!(self.data.inverse_valid()); let n = self.n_qubits(); let stride = self.data.stride(); @@ -381,8 +404,12 @@ impl Tableau { // `ω(Z_a, g) = x_g[a]`: one column per half, minus the pivot, which is // not multiplied into itself. let mut selected = [destab_sel, stab_sel]; - for (half, out) in HALVES.into_iter().zip(selected.iter_mut()) { - self.data.gather_column(half, Plane::X, addr0, out); + match destab_column { + Some(column) => selected[Half::Destab as usize].copy_from_slice(column), + None => selected[Half::Destab as usize].fill(0), + } + selected[Half::Stab as usize].copy_from_slice(stab_column); + for out in selected.iter_mut() { TableauData::set_bit(out, pivot, false); } @@ -479,9 +506,10 @@ impl Tableau { self.data.inv_sign_plane_mut(InvRow::Z), ); debug_assert!( - p.stab.0.iter().all(|&w| w == 0) - && blocks::first_set(p.stab.1) == Some(addr0) - && blocks::count_set(p.stab.1) == 1, + destab_column.is_none() + || p.stab.0.iter().all(|&w| w == 0) + && blocks::first_set(p.stab.1) == Some(addr0) + && blocks::count_set(p.stab.1) == 1, "the appends must leave the pivot stabilizer at ±Z_{addr0}" ); @@ -490,11 +518,13 @@ impl Tableau { // `iz_a = U'†Z_aU' = (−1)^r Z_p`, so that row's sign is the outcome. If // it disagrees, `append_X(p)` flips the sign of every inverse row with a // `Z` at site `p` — `iz_a` among them — which is exactly negating `s_p`. - if (self.data.inv_sign(InvRow::Z, addr0) == 2) != outcome { + let flip = (self.data.inv_sign(InvRow::Z, addr0) == 2) != outcome; + if flip { blocks::pauli_x(p.destab.1, self.data.inv_sign_plane_mut(InvRow::X)); blocks::pauli_x(p.destab.0, self.data.inv_sign_plane_mut(InvRow::Z)); } self.data.restore_inv_scratch(scratch); + flip } } diff --git a/crates/ppvm-tableau-2/src/inverse_tests.rs b/crates/ppvm-tableau-2/src/inverse_tests.rs index 2305ab75c..5b906c77a 100644 --- a/crates/ppvm-tableau-2/src/inverse_tests.rs +++ b/crates/ppvm-tableau-2/src/inverse_tests.rs @@ -218,3 +218,30 @@ fn scramble_frame(tab: &mut Tableau, rng: &mut SmallRng, gates: usize) { apply_indexed_gate(tab, rng.random_range(0..limit), a, b); } } + +/// Stim's collapse in `measure_batch` must keep the inverse signs exact and +/// leave the same state as the textbook projection: equal outcomes now, and +/// equal outcomes after a further shared circuit. +#[test] +fn stabilizer_collapse_matches_the_projection() { + for (n, seed) in [(1usize, 50u64), (6, 51), (64, 52), (70, 53)] { + let mut rng = SmallRng::seed_from_u64(seed); + let mut a = GeneralizedTableau::::new(n, 1e-12); + let mut b = GeneralizedTableau::::new(n, 1e-12); + let half: Vec = (0..n).step_by(2).collect(); + let all: Vec = (0..n).collect(); + for round in 0..3 { + let circuit = rng.random::(); + scramble_frame(&mut a.tableau, &mut SmallRng::seed_from_u64(circuit), 8 * n); + scramble_frame(&mut b.tableau, &mut SmallRng::seed_from_u64(circuit), 8 * n); + let shot = rng.random::(); + let (mut ar, mut br) = (SmallRng::seed_from_u64(shot), SmallRng::seed_from_u64(shot)); + let targets = if round == 2 { &all } else { &half }; + let ra = a.measure_batch(targets, &mut ar); + let rb: Vec<_> = targets.iter().map(|&q| b.measure(q, &mut br)).collect(); + assert_eq!(ra, rb, "n={n} round={round}"); + a.tableau.assert_inverse_consistent(); + assert!(a.tableau.inverse_valid()); + } + } +} diff --git a/crates/ppvm-tableau-2/src/measure.rs b/crates/ppvm-tableau-2/src/measure.rs index 12074eef3..5821a94be 100644 --- a/crates/ppvm-tableau-2/src/measure.rs +++ b/crates/ppvm-tableau-2/src/measure.rs @@ -50,9 +50,10 @@ use ppvm_traits_2::{Clifford, Measure, Pauli, Reset}; use rand::{Rng, RngExt}; use crate::data::{ - Bitstring, COMPLEX_PHASE_CONVERSION, GeneralizedTableau, Tableau, + Bitstring, COMPLEX_PHASE_CONVERSION, GeneralizedTableau, HalfColumns, Tableau, bits_to_index, compute_phase_with_mask_static, symplectic_inner, }; +use crate::storage::{Half, blocks}; /// The pure Clifford measurement procedure. /// @@ -156,6 +157,7 @@ pub struct MeasureScratch { a: Vec<(I, Complex64)>, bt: Vec<(I, Complex64)>, merged: Vec<(I, Complex64)>, + columns: HalfColumns, } impl MeasureScratch { @@ -170,6 +172,7 @@ impl MeasureScratch { a: Vec::new(), bt: Vec::new(), merged: Vec::new(), + columns: HalfColumns::default(), } } } @@ -200,11 +203,7 @@ impl Measure for GeneralizedTableau { return None; } - let decomposition = self.compute_decomposition(qubit, Pauli::Z); - - self.with_scratch(|s, scratch| { - s.measure_with_scratch(qubit, scratch, decomposition, true, rng) - }) + self.with_scratch(|s, scratch| s.measure_z_with_scratch(qubit, scratch, true, rng)) } /// Override the trait default (a per-target `measure` loop) with one scratch @@ -333,8 +332,123 @@ impl GeneralizedTableau { self.measurement_record.push(None); return None; } - let decomposition = self.compute_decomposition(idx, Pauli::Z); - self.measure_with_scratch(idx, scratch, decomposition, true, rng) + self.measure_z_with_scratch(idx, scratch, true, rng) + } + + /// Measure `indices` in the Z basis the way Stim's `collapse_z` does: the + /// random ones first, under one row guard, then the deterministic ones in + /// the canonical orientation. Records are pushed in the caller's order. + /// + /// A target is random when some stabilizer anticommutes with `Z`, a cheap + /// contiguous check before any transpose. Collapsing one target can make a + /// later one deterministic but never the reverse, so the split is exact and + /// a frame with no random target never re-orients. + /// + /// Unlike [`Self::measure_many`], the RNG-draw order is not that of a + /// per-target loop: random targets draw first. Measurements on distinct + /// qubits commute, so the outcome distribution is the same. On a stabilizer + /// state the collapse is Stim's, so the frame afterwards can differ from + /// [`Self::measure_many`]'s while describing the same state. + pub fn measure_batch( + &mut self, + indices: &[usize], + rng: &mut R, + ) -> Vec> { + self.with_scratch(|s, scratch| s.measure_batch_with_scratch(indices, scratch, rng)) + } + + /// [`Self::measure_batch`] with a caller-supplied scratch. + pub fn measure_batch_with_scratch( + &mut self, + indices: &[usize], + scratch: &mut MeasureScratch, + rng: &mut R, + ) -> Vec> { + let random: Vec = indices + .iter() + .map(|&q| !self.is_lost[q] && self.tableau.find_z_anticommuting_stabilizer(q).is_some()) + .collect(); + let mut outcomes = vec![None; indices.len()]; + if random.contains(&true) { + self.with_row_major(|s| { + for (k, &q) in indices.iter().enumerate().filter(|&(k, _)| random[k]) { + outcomes[k] = s.measure_batch_one(q, scratch, rng); + } + }); + } + for (k, &q) in indices.iter().enumerate().filter(|&(k, _)| !random[k]) { + outcomes[k] = self.measure_batch_one(q, scratch, rng); + } + self.measurement_record.extend_from_slice(&outcomes); + outcomes + } + + /// One unrecorded target of [`Self::measure_batch`]. On a stabilizer state + /// (one amplitude, at index 0) this is Stim's collapse, which never reads + /// the destabilizer column; otherwise the general kernel. + /// + /// Both draw the same `random::() < 0.5` for a random outcome, so the + /// outcomes match; only the frame chosen for the post-measurement state + /// differs. The amplitude is untouched: the frame alone carries the state. + fn measure_batch_one( + &mut self, + idx: usize, + scratch: &mut MeasureScratch, + rng: &mut R, + ) -> Option { + if self.is_lost[idx] { + return None; + } + let stabilizer_state = + self.tableau.inverse_valid() && self.coefficients.iter().all(|&(_, i)| i == I::zero()); + if !stabilizer_state { + return self.measure_z_with_scratch(idx, scratch, false, rng); + } + let mut column = std::mem::take(&mut scratch.columns[Half::Stab as usize]); + self.tableau + .z_anticommuting_stabilizers_into(idx, &mut column); + let n = self.n_qubits(); + let outcome = match blocks::first_set(&column).filter(|&p| p < n) { + None => self.tableau.inverse_outcome(idx), + Some(pivot) => { + let outcome = rng.random::() < 0.5; + self.tableau.collapse_z(idx, pivot, outcome, &column); + scratch.odd_phase_mask = None; + outcome + } + }; + scratch.columns[Half::Stab as usize] = column; + Some(outcome) + } + + /// Measure `Z_qubit` on a qubit that is not lost. Its anticommutation + /// columns are gathered once, for the decomposition and the projection. + pub(crate) fn measure_z_with_scratch( + &mut self, + qubit: usize, + scratch: &mut MeasureScratch, + record: bool, + rng: &mut R, + ) -> Option { + let mut columns = std::mem::take(&mut scratch.columns); + self.tableau + .z_anticommutation_columns_into(qubit, &mut columns); + let phase = self.decomposition_phase_with(qubit, Pauli::Z, &columns); + // A deterministic outcome reads the masks only through amplitude + // indices; with every index zero (a Clifford run) skip widening them. + let [destab, stab] = &columns; + let decomposition = if stab.iter().all(|&w| w == 0) + && self.coefficients.iter().all(|&(_, idx)| idx == I::zero()) + { + (phase, I::zero(), I::zero()) + } else { + let n = self.n_qubits(); + (phase, bits_to_index(stab, n), bits_to_index(destab, n)) + }; + let outcome = + self.measure_with_scratch(qubit, scratch, decomposition, &columns, record, rng); + scratch.columns = columns; + outcome } /// The coefficient-aware measurement kernel. @@ -372,6 +486,7 @@ impl GeneralizedTableau { addr0: usize, scratch: &mut MeasureScratch, decomposition: (u8, I, I), + columns: &HalfColumns, record: bool, rng: &mut R, ) -> Option { @@ -578,7 +693,7 @@ impl GeneralizedTableau { self.coefficients.normalize(); self.tableau - .update_tableau_according_to_outcome(addr0, q_idx, outcome); + .update_tableau_with_columns(addr0, q_idx, outcome, columns); // Destabilizer phases just changed; invalidate the cached mask. scratch.odd_phase_mask = None; if record { diff --git a/crates/ppvm-tableau-2/src/noise.rs b/crates/ppvm-tableau-2/src/noise.rs index 8de07cefe..28573bb79 100644 --- a/crates/ppvm-tableau-2/src/noise.rs +++ b/crates/ppvm-tableau-2/src/noise.rs @@ -29,7 +29,7 @@ use num::Zero; use ppvm_traits_2::{ AsymmetricLossChannel, Clifford, CorrelatedLossChannel, Depolarizing, Depolarizing2, - LossChannel, Pauli, PauliError, ResetLossChannel, TwoQubitPauliError, + LossChannel, PauliError, ResetLossChannel, TwoQubitPauliError, }; use rand::{Rng, RngExt}; @@ -251,8 +251,7 @@ impl GeneralizedTableau { let outcome = if self.is_lost[qubit] { None } else { - let decomposition = self.compute_decomposition(qubit, Pauli::Z); - self.measure_with_scratch(qubit, &mut MeasureScratch::new(), decomposition, false, rng) + self.measure_z_with_scratch(qubit, &mut MeasureScratch::new(), false, rng) }; if let Some(true) = outcome { Clifford::x(self, qubit); @@ -323,9 +322,8 @@ impl AsymmetricLossChannel for GeneralizedTableau { // scratch cannot serve this site again. Run the same projection kernel // with ephemeral buffers, as the legacy path did, while keeping the // intentionally observable record append. - let decomposition = self.compute_decomposition(qubit, Pauli::Z); if let Some(true) = - self.measure_with_scratch(qubit, &mut MeasureScratch::new(), decomposition, true, rng) + self.measure_z_with_scratch(qubit, &mut MeasureScratch::new(), true, rng) { Clifford::x(self, qubit); } diff --git a/crates/ppvm-tableau-2/src/storage/blocks.rs b/crates/ppvm-tableau-2/src/storage/blocks.rs index b25b8b08e..94005f8fe 100644 --- a/crates/ppvm-tableau-2/src/storage/blocks.rs +++ b/crates/ppvm-tableau-2/src/storage/blocks.rs @@ -134,13 +134,18 @@ pub(crate) fn sqrt_y_dag(x: &mut [u64], z: &mut [u64], ph: &mut [u64]) { #[inline] pub(crate) fn cnot(xc: &[u64], zc: &mut [u64], xt: &mut [u64], zt: &[u64], ph: &mut [u64]) { for i in 0..ph.len() { - let (a, b, c, d) = (xc[i], zc[i], xt[i], zt[i]); - ph[i] ^= a & d & !(c ^ b); - zc[i] = b ^ d; - xt[i] = c ^ a; + cnot_word(xc[i], &mut zc[i], &mut xt[i], zt[i], &mut ph[i]); } } +/// One word of [`cnot`]. +#[inline(always)] +pub(crate) fn cnot_word(xc: u64, zc: &mut u64, xt: &mut u64, zt: u64, ph: &mut u64) { + *ph ^= xc & zt & !(*xt ^ *zc); + *zc ^= zt; + *xt ^= xc; +} + /// `CZ`: `z_a ^= x_b`, `z_b ^= x_a`, sign flips where `x_a & x_b & (z_a ^ z_b)`. #[inline] pub(crate) fn cz(xa: &[u64], za: &mut [u64], xb: &[u64], zb: &mut [u64], ph: &mut [u64]) { @@ -184,18 +189,26 @@ pub(crate) fn row_multiply( src_x: &[u64], src_z: &[u64], ) -> u8 { - let mut sign_count = 0u32; - let mut imag_count = 0u32; - for i in 0..dst_x.len() { - let (a, b, c, d) = (dst_x[i], dst_z[i], src_x[i], src_z[i]); - let sign = (a & b & c & !d) | (a & !b & !c & d) | (!a & b & c & d); - let imag = (a & !b & d) | (a & !c & d) | (!a & b & c) | (b & c & !d); - sign_count += sign.count_ones(); - imag_count += imag.count_ones(); - dst_x[i] = a ^ c; - dst_z[i] = b ^ d; + // Short rows: one fused serial pass. Longer: phase first, then the bits — + // two simple passes vectorize, one fused does not. + if dst_x.len() < 8 { + let mut phase = PhaseCounter::default(); + for i in 0..dst_x.len() { + let (a, b, c, d) = (dst_x[i], dst_z[i], src_x[i], src_z[i]); + phase.add(a, b, c, d); + dst_x[i] = a ^ c; + dst_z[i] = b ^ d; + } + return phase.total(); + } + let phase = row_multiply_phase(dst_x, dst_z, src_x, src_z); + for (d, s) in dst_x.iter_mut().zip(src_x) { + *d ^= s; + } + for (d, s) in dst_z.iter_mut().zip(src_z) { + *d ^= s; } - ((2 * sign_count + imag_count) % 4) as u8 + phase } /// [`row_multiply`]'s phase without its bit writes. @@ -204,18 +217,63 @@ pub(crate) fn row_multiply( /// phase of a product of two rows whose *bits* the forward gate already /// maintains, so there is nothing to write back — and the operands are borrowed /// forward majors, which a writing kernel could not take. -#[inline] +#[inline(always)] pub(crate) fn row_multiply_phase(a_x: &[u64], a_z: &[u64], b_x: &[u64], b_z: &[u64]) -> u8 { - let mut sign_count = 0u32; - let mut imag_count = 0u32; - for i in 0..a_x.len() { - let (a, b, c, d) = (a_x[i], a_z[i], b_x[i], b_z[i]); - let sign = (a & b & c & !d) | (a & !b & !c & d) | (!a & b & c & d); - let imag = (a & !b & d) | (a & !c & d) | (!a & b & c) | (b & c & !d); - sign_count += sign.count_ones(); - imag_count += imag.count_ones(); + // Below two chunks one serial counter is cheapest. From there, four + // counters over 4-word chunks are independent lanes the compiler can + // vectorize (a serial counter cannot be); the tail uses one more counter. + if a_x.len() < 8 { + let mut phase = PhaseCounter::default(); + for i in 0..a_x.len() { + phase.add(a_x[i], a_z[i], b_x[i], b_z[i]); + } + return phase.total(); + } + let (ax, ax_tail) = a_x.as_chunks::<4>(); + let (az, az_tail) = a_z.as_chunks::<4>(); + let (bx, bx_tail) = b_x.as_chunks::<4>(); + let (bz, bz_tail) = b_z.as_chunks::<4>(); + let mut lanes = [PhaseCounter::default(); 4]; + for (((a, b), c), d) in ax.iter().zip(az).zip(bx).zip(bz) { + for k in 0..4 { + lanes[k].add(a[k], b[k], c[k], d[k]); + } + } + let mut tail = PhaseCounter::default(); + for i in 0..ax_tail.len() { + tail.add(ax_tail[i], az_tail[i], bx_tail[i], bz_tail[i]); + } + let lanes: u32 = lanes.iter().map(|c| u32::from(c.total())).sum(); + ((lanes + u32::from(tail.total())) % 4) as u8 +} + +/// The `g`-rule phase of a row product, summed mod 4 with one 2-bit counter +/// per bit position (`lo` + 2·`hi`), so a product popcounts once at the end +/// rather than twice per word. This is Stim's `cnt1` / `cnt2` update; it adds +/// exactly the Aaronson–Gottesman `g` term at every bit (checked exhaustively +/// in the tests below). +#[derive(Clone, Copy, Default)] +pub(crate) struct PhaseCounter { + lo: u64, + hi: u64, +} + +impl PhaseCounter { + /// Add the term of `(a, b) · (c, d)` — `x`/`z` bits of the left and right + /// factor — at every bit. + #[inline(always)] + pub(crate) fn add(&mut self, a: u64, b: u64, c: u64, d: u64) { + let ad = a & d; + let anticommute = (c & b) ^ ad; + self.hi ^= (self.lo ^ a ^ c ^ b ^ d ^ ad) & anticommute; + self.lo ^= anticommute; + } + + /// The accumulated phase, mod 4. + #[inline(always)] + pub(crate) fn total(self) -> u8 { + ((self.lo.count_ones() + 2 * self.hi.count_ones()) % 4) as u8 } - ((2 * sign_count + imag_count) % 4) as u8 } // ─── Column-wise row multiplication ─────────────────────────────────────── @@ -463,3 +521,35 @@ pub(crate) fn first_set(words: &[u64]) -> Option { .position(|&w| w != 0) .map(|i| i * super::BITS_PER_WORD + words[i].trailing_zeros() as usize) } + +#[cfg(test)] +mod tests { + use super::PhaseCounter; + + /// The Aaronson–Gottesman `g` term at one bit, as the replaced per-word + /// formula computed it. + fn g_term(a: bool, b: bool, c: bool, d: bool) -> u8 { + let sign = (a && b && c && !d) || (a && !b && !c && d) || (!a && b && c && d); + let imag = (a && !b && d) || (a && !c && d) || (!a && b && c) || (b && c && !d); + (2 * u8::from(sign) + u8::from(imag)) % 4 + } + + #[test] + fn phase_counter_adds_the_g_term_from_every_state() { + for bits in 0..16u8 { + let [a, b, c, d] = [0, 1, 2, 3].map(|k| bits >> k & 1 == 1); + for start in 0..4u8 { + let mut counter = PhaseCounter { + lo: u64::from(start & 1), + hi: u64::from(start >> 1), + }; + counter.add(a.into(), b.into(), c.into(), d.into()); + assert_eq!( + counter.total(), + (start + g_term(a, b, c, d)) % 4, + "{bits:04b} from {start}" + ); + } + } + } +} diff --git a/crates/ppvm-tableau-2/src/storage/inverse.rs b/crates/ppvm-tableau-2/src/storage/inverse.rs index 68d2bfbe3..eb31ff25c 100644 --- a/crates/ppvm-tableau-2/src/storage/inverse.rs +++ b/crates/ppvm-tableau-2/src/storage/inverse.rs @@ -236,7 +236,9 @@ impl TableauData { pub(crate) fn inv_pair_phase(&self, a: (InvRow, usize), b: (InvRow, usize)) -> u8 { let (ax, az) = self.inv_planes(a.0, a.1); let (bx, bz) = self.inv_planes(b.0, b.1); - let g = blocks::row_multiply_phase(ax, az, bx, bz); + // Padding words are zero and contribute no phase. + let live = self.live_words(); + let g = blocks::row_multiply_phase(&ax[..live], &az[..live], &bx[..live], &bz[..live]); (self.inv_sign(a.0, a.1) + self.inv_sign(b.0, b.1) + g) % 4 } diff --git a/crates/ppvm-tableau-2/src/storage/mod.rs b/crates/ppvm-tableau-2/src/storage/mod.rs index af7810aae..4ed93e17e 100644 --- a/crates/ppvm-tableau-2/src/storage/mod.rs +++ b/crates/ppvm-tableau-2/src/storage/mod.rs @@ -150,6 +150,8 @@ pub(crate) struct TableauData { n_qubits: usize, /// Words per major, `n.div_ceil(64)` rounded up to a whole block. stride: usize, + /// `n.div_ceil(64)`: the words of a major that can hold a set bit. + live: usize, orientation: Orientation, /// Signs of the inverse tableau's rows — a derived cache whose bits live in /// the quadrants above. See [`inverse`]; excluded from equality and hashing. @@ -166,6 +168,7 @@ impl TableauData { blocks: vec![Block([0; WORDS_PER_BLOCK]); words / WORDS_PER_BLOCK], n_qubits, stride, + live: n_qubits.div_ceil(BITS_PER_WORD), orientation: Orientation::ColumnMajor, inverse: InverseSigns::identity(stride), }; @@ -450,6 +453,37 @@ impl TableauData { // ─── Disjoint borrows for the gate kernels ──────────────────────────── + /// Words of a major that can hold a set bit, `n.div_ceil(64)`. The rest of + /// the stride is zero padding, and every gate kernel maps all-zero words to + /// all-zero words, so the gate borrows stop here. + #[inline] + fn live_words(&self) -> usize { + self.live + } + + /// Borrow `ranges` of the arena mutably at once, trimmed to + /// [`Self::live_words`]. The caller passes ranges that are disjoint by + /// layout (distinct majors, or a major and a phase plane), so the + /// overlap check `get_disjoint_mut` repeats on every gate is a debug check. + #[inline] + fn live_disjoint_mut( + &mut self, + ranges: [std::ops::Range; N], + ) -> [&mut [u64]; N] { + let live = self.live_words(); + let ranges = ranges.map(|r| r.start..r.start + live); + let len = self.blocks.len() * WORDS_PER_BLOCK; + debug_assert!(ranges.iter().all(|r| r.end <= len)); + debug_assert!(ranges.iter().enumerate().all(|(i, a)| { + ranges[i + 1..] + .iter() + .all(|b| a.end <= b.start || b.end <= a.start) + })); + // SAFETY: every range is in bounds and no two overlap (checked above in + // debug builds; guaranteed by the arena layout). + unsafe { self.words_mut().get_disjoint_unchecked_mut(ranges) } + } + /// The `(X, Z, phase-hi)` slices a one-qubit Clifford sweeps, for one half. /// /// Column-major only: `q` indexes a qubit column. @@ -460,15 +494,11 @@ impl TableauData { q: usize, ) -> (&mut [u64], &mut [u64], &mut [u64]) { debug_assert_eq!(self.orientation, Orientation::ColumnMajor); - let ranges = [ + let [x, z, ph] = self.live_disjoint_mut([ self.major_range(half, Plane::X, q), self.major_range(half, Plane::Z, q), self.phase_range(half, true), - ]; - let [x, z, ph] = self - .words_mut() - .get_disjoint_mut(ranges) - .expect("quadrant and phase-plane regions are disjoint by construction"); + ]); (x, z, ph) } @@ -478,20 +508,56 @@ impl TableauData { pub(crate) fn gate2_mut(&mut self, half: Half, a: usize, b: usize) -> Gate2Slices<'_> { debug_assert_eq!(self.orientation, Orientation::ColumnMajor); debug_assert_ne!(a, b, "two-qubit gate needs distinct qubits"); - let ranges = [ + let [xa, za, xb, zb, ph] = self.live_disjoint_mut([ self.major_range(half, Plane::X, a), self.major_range(half, Plane::Z, a), self.major_range(half, Plane::X, b), self.major_range(half, Plane::Z, b), self.phase_range(half, true), - ]; - let [xa, za, xb, zb, ph] = self - .words_mut() - .get_disjoint_mut(ranges) - .expect("distinct qubit columns and the phase plane are disjoint"); + ]); (xa, za, xb, zb, ph) } + /// `CNOT(c, t)` on both halves from one borrow of the live words, returning the + /// `g`-rule terms of the inverse-row products `ix_c·ix_t` and `iz_c·iz_t`, + /// read before any write. An inverse `X` row is the two halves' `Z` columns + /// and a `Z` row their `X` columns (see [`inverse`]), so these are the same + /// eight columns the forward update touches. + pub(crate) fn cnot_fused(&mut self, c: usize, t: usize) -> (u8, u8) { + debug_assert_eq!(self.orientation, Orientation::ColumnMajor); + debug_assert_ne!(c, t, "two-qubit gate needs distinct qubits"); + let [sxc, szc, sxt, szt, sph, dxc, dzc, dxt, dzt, dph] = self.live_disjoint_mut([ + self.major_range(Half::Stab, Plane::X, c), + self.major_range(Half::Stab, Plane::Z, c), + self.major_range(Half::Stab, Plane::X, t), + self.major_range(Half::Stab, Plane::Z, t), + self.phase_range(Half::Stab, true), + self.major_range(Half::Destab, Plane::X, c), + self.major_range(Half::Destab, Plane::Z, c), + self.major_range(Half::Destab, Plane::X, t), + self.major_range(Half::Destab, Plane::Z, t), + self.phase_range(Half::Destab, true), + ]); + // One live word (n <= 64): a single scalar pass saves the four loops' + // overhead. Otherwise separate loops, each simple enough to vectorize. + if sph.len() == 1 { + let (mut g_x, mut g_z) = ( + blocks::PhaseCounter::default(), + blocks::PhaseCounter::default(), + ); + g_x.add(szc[0], dzc[0], szt[0], dzt[0]); + g_z.add(sxc[0], dxc[0], sxt[0], dxt[0]); + blocks::cnot_word(sxc[0], &mut szc[0], &mut sxt[0], szt[0], &mut sph[0]); + blocks::cnot_word(dxc[0], &mut dzc[0], &mut dxt[0], dzt[0], &mut dph[0]); + return (g_x.total(), g_z.total()); + } + let g_x = blocks::row_multiply_phase(szc, dzc, szt, dzt); + let g_z = blocks::row_multiply_phase(sxc, dxc, sxt, dxt); + blocks::cnot(sxc, szc, sxt, szt, sph); + blocks::cnot(dxc, dzc, dxt, dzt, dph); + (g_x, g_z) + } + // ─── Logical bit access ─────────────────────────────────────────────── /// Read bit `i` of `words`. diff --git a/crates/ppvm-tableau-2/tests/behaviour.rs b/crates/ppvm-tableau-2/tests/behaviour.rs index 73075eb28..703c87eb8 100644 --- a/crates/ppvm-tableau-2/tests/behaviour.rs +++ b/crates/ppvm-tableau-2/tests/behaviour.rs @@ -858,3 +858,70 @@ fn wide_msd_shaped_circuit_runs_and_agrees_with_the_naive_form() { ); } } + +/// On a stabilizer state only random measurements draw, so `measure_batch` +/// draws in the same order as a per-qubit loop and must match it exactly, even +/// when a deterministic target precedes the random one it depends on. It must +/// also leave the same state behind. +#[test] +fn measure_batch_matches_a_per_qubit_loop_on_a_stabilizer_state() { + let mut base: Tab = GeneralizedTableau::new(6, 1e-10); + let mut setup_rng = rng(3); + base.h(0); + base.cnot(0, 1); + base.h(2); + base.x(4); + base.loss_channel(5, 1.0, &mut setup_rng); + let targets = [3, 1, 0, 4, 2, 5, 0]; + + for seed in 0..16 { + let mut a = base.fork(); + let mut b = base.fork(); + let (mut ar, mut br) = (rng(seed), rng(seed)); + let ra = a.measure_batch(&targets, &mut ar); + let rb: Vec> = targets.iter().map(|&q| b.measure(q, &mut br)).collect(); + + assert_eq!(ra, rb); + assert_eq!(ra[1], ra[2], "the Bell pair must agree"); + assert_eq!(ra[2], ra[6], "a repeated target must repeat its outcome"); + assert_eq!((ra[0], ra[3], ra[5]), (Some(false), Some(true), None)); + assert_eq!( + a.current_measurement_record(), + b.current_measurement_record() + ); + // The frames may differ (Stim's collapse picks another basis for the + // same state), so compare every Pauli expectation value instead. + for code in 0..4usize.pow(6) { + let w: String = (0..6) + .map(|q| b"IXYZ"[(code >> (2 * q)) & 3] as char) + .collect(); + let (ea, eb) = (a.expectation(&word(&w)), b.expectation(&word(&w))); + assert_close(ea, eb, 1e-12); + } + } +} + +/// With several amplitudes a deterministic-frame target also draws, so the +/// order differs from a per-qubit loop; the distribution must not. +#[test] +fn measure_batch_samples_the_per_qubit_distribution() { + let mut base: Tab = GeneralizedTableau::new(3, 1e-10); + base.h(0); + base.t(0); + base.h(0); // P(1) = sin²(π/8), from a frame where Z₀ is a stabilizer + base.h(1); + base.cnot(1, 2); + + let shots = 4000; + let (mut ones0, mut ones1) = (0, 0); + for seed in 0..shots { + let mut tab = base.fork(); + let r = tab.measure_batch(&[0, 2, 1], &mut rng(seed)); + assert_eq!(r[1], r[2], "the Bell pair must agree"); + ones0 += usize::from(r[0] == Some(true)); + ones1 += usize::from(r[2] == Some(true)); + } + let p1 = (std::f64::consts::PI / 8.0).sin().powi(2); + assert_close(ones0 as f64 / shots as f64, p1, 0.03); + assert_close(ones1 as f64 / shots as f64, 0.5, 0.04); +}