diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6d362c6a0..8b17c8dec 100644
--- 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
index 0b84b0ac0..a5cc637fa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org> 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
index c87932fd4..efac990e4 100644
--- 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