Improve performance of random_pauli! by 100x (QuantumSavory#454)

This commit improves the performance of the unbiased branch of `random_pauli!` by a factor of about 100. Each chunk in the data array is replaced with a random `UInt64`. --------- Co-authored-by: Stefan Krastanov <[email protected]>
amicciche · Dec 29, 2024 · 3297023 · 3297023
1 parent 6065fab
commit 3297023
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 7 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,8 +5,9 @@
 
 # News
 
-## v0.9.16 - dev
+## v0.9.16 - 2024-12-29
 
+- 100× faster unbiased `random_pauli`.
 - Enhancements to `GF(2)` Linear Algebra: unexported, experimental `gf2_row_echelon_with_pivots!`, `gf2_nullspace`, `gf2_rowspace_basis`.
 
 ## v0.9.15 - 2024-12-22

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]> and QuantumSavory community members"]
-version = "0.9.15"
+version = "0.9.16"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

diff --git a/src/randoms.jl b/src/randoms.jl
@@ -1,4 +1,4 @@
-using Random: randperm, AbstractRNG, GLOBAL_RNG
+using Random: randperm, AbstractRNG, GLOBAL_RNG, rand!
 using ILog2
 import Nemo
 
@@ -21,11 +21,24 @@ function random_pauli end
 """An in-place version of [`random_pauli`](@ref)"""
 function random_pauli! end
 
+
+# Modified from `Base` for `BitArray`
+@inline _msk_end(::Type{T}, l::Int) where {T<:Unsigned} = ~T(0) >>> _mod(T, -l)
+
+# A mask for the non-coding bits in the last chunk.
+@inline _msk_end(P::PauliOperator) = _msk_end(eltype(P.xz), length(P))
+
+# Unset the leftover bits in the data array that don't code the Pauli operator.
+@inline function _unset_noncoding_bits!(P::PauliOperator)
+    msk = _msk_end(P)
+    P.xz[end] &= msk
+    P.xz[end÷2] &= msk
+    nothing
+end
+
 function random_pauli!(rng::AbstractRNG, P::PauliOperator; nophase=true, realphase=true)
-    n = nqubits(P)
-    for i in 1:n
-        P[i] = rand(rng, (true, false)), rand(rng, (true,false))
-    end
+    rand!(rng, P.xz)
+    _unset_noncoding_bits!(P)
     P.phase[] = nophase ? 0x0 : (realphase ? rand(rng,(0x0,0x2)) : rand(rng,0x0:0x3))
     P
 end