diff --git a/.gitignore b/.gitignore
index fa59c73..e4f6948 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,3 +28,5 @@
 *.synctex.gz
 *.synctex.gz(busy)
 *.pdfsync
+test_data/
+__pycache__/
\ No newline at end of file
diff --git a/FFI.lean b/FFI.lean
new file mode 100644
index 0000000..75d1d9f
--- /dev/null
+++ b/FFI.lean
@@ -0,0 +1,7 @@
+module  -- shake: keep-all --deprecated_module: ignore
+
+public import FFI.FMA
+public import FFI.Float
+public import FFI.Float32
+public import FFI.Model.Float
+public import FFI.Model.Float32
diff --git a/FFI/FMA.lean b/FFI/FMA.lean
new file mode 100644
index 0000000..f663f99
--- /dev/null
+++ b/FFI/FMA.lean
@@ -0,0 +1,57 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+prelude
+public import Init.Data.Float.Model.Unpacked.Round
+
+-- This file is part of the logical model for floats which authors of float libraries
+-- need to rely on.
+@[expose] public section
+
+namespace Float.Model.UnpackedFloat
+
+/--
+Computes the fused multiply-add `x * y + z` of three floating point numbers and rounds the
+result according to the given specification. The product `x * y` is exact and never rounded
+on its own; only the final sum is rounded.
+-/
+def fma (spec : Format) : UnpackedFloat → UnpackedFloat → UnpackedFloat → UnpackedFloat
+  | .notANumber, _, _ => .notANumber
+  | _, .notANumber, _ => .notANumber
+  | _, _, .notANumber => .notANumber
+  | .zero _, .infinity _, _ => .notANumber
+  | .infinity _, .zero _, _ => .notANumber
+  | .infinity sign₁, .infinity sign₂, .infinity sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .infinity sign₃ else .notANumber
+  | .infinity sign₁, .finite sign₂ .., .infinity sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .infinity sign₃ else .notANumber
+  | .finite sign₁ .., .infinity sign₂, .infinity sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .infinity sign₃ else .notANumber
+  | .infinity sign₁, .infinity sign₂, _ => .infinity (sign₁ * sign₂)
+  | .infinity sign₁, .finite sign₂ .., _ => .infinity (sign₁ * sign₂)
+  | .finite sign₁ .., .infinity sign₂, _ => .infinity (sign₁ * sign₂)
+  | _, _, .infinity sign₃ => .infinity sign₃
+  | .zero sign₁, .zero sign₂, .zero sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .zero sign₃ else .zero .positive
+  | .zero sign₁, .finite sign₂ .., .zero sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .zero sign₃ else .zero .positive
+  | .finite sign₁ .., .zero sign₂, .zero sign₃ =>
+    if sign₁ * sign₂ == sign₃ then .zero sign₃ else .zero .positive
+  | .zero _, _, z => z
+  | _, .zero _, z => z
+  | .finite s₁ m₁ e₁ _, .finite s₂ m₂ e₂ _, .zero _ =>
+    roundWithAccuracy spec (s₁ * s₂) (m₁ * m₂) (e₁ + e₂) .exact
+  | .finite s₁ m₁ e₁ _, .finite s₂ m₂ e₂ _, .finite s₃ m₃ e₃ _ =>
+    let productMantissa := m₁ * m₂
+    let productExponent := e₁ + e₂
+    let smallerExponent := min productExponent e₃
+    let (productMantissa, _) := decreaseExponent productMantissa productExponent smallerExponent
+    let (m₃, _) := decreaseExponent m₃ e₃ smallerExponent
+    let mantissa := (s₁ * s₂).apply productMantissa + s₃.apply m₃
+    normalize spec mantissa smallerExponent .positive
+
+end Float.Model.UnpackedFloat
diff --git a/FFI/Float.lean b/FFI/Float.lean
new file mode 100644
index 0000000..292bb9f
--- /dev/null
+++ b/FFI/Float.lean
@@ -0,0 +1,31 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+prelude
+public import Init.Data.Float.Float
+public import FFI.Model.Float
+
+@[expose] public section
+
+namespace Float
+
+/--
+Computes the fused multiply-add `x * y + z` of three floating-point numbers. This operation is
+performed with a single rounding, which can be more accurate than performing the multiplication and
+addition separately.
+
+This function has a logical model in terms of `Float.Model`. It is implemented in compiled code by
+the C function `fma`.
+-/
+@[extern "fma"] def fma : Float → Float → Float → Float :=
+  fun x y z => .ofModel (x.toModel.fma y.toModel z.toModel)
+
+/-- `log (1 + x)`, the C99 `log1p`, accurate for small `x` where `log (1 + x)` loses the leading
+digits of the result to the rounding of `1 + x`. -/
+@[extern "log1p"] opaque log1p (x : Float) : Float
+
+end Float
diff --git a/FFI/Float32.lean b/FFI/Float32.lean
new file mode 100644
index 0000000..d51bb0e
--- /dev/null
+++ b/FFI/Float32.lean
@@ -0,0 +1,29 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+prelude
+public import Init.Data.Float.Float32
+public import FFI.Model.Float32
+
+@[expose] public section
+
+namespace Float32
+
+/--
+Computes the fused multiply-add `x * y + z` of three floating-point numbers. This operation is performed with a single rounding, which can be more accurate than performing the multiplication and addition separately.
+
+This function has a logical model in terms of `Float32.Model`. It is implemented in compiled code
+by the C function `fmaf`.
+-/
+@[extern "fmaf"] def fma : Float32 → Float32 → Float32 → Float32 :=
+  fun x y z => .ofModel (x.toModel.fma y.toModel z.toModel)
+
+/-- `log (1 + x)`, the C99 `log1p`, accurate for small `x` where `log (1 + x)` loses the leading
+digits of the result to the rounding of `1 + x`. -/
+@[extern "log1pf"] opaque log1p (x : Float32) : Float32
+
+end Float32
diff --git a/FFI/Model/Float.lean b/FFI/Model/Float.lean
new file mode 100644
index 0000000..4475657
--- /dev/null
+++ b/FFI/Model/Float.lean
@@ -0,0 +1,24 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+prelude
+public import Init.Data.Float.Model.Float
+public import FFI.FMA
+
+-- This file is part of the logical model for floats which authors of float libraries
+-- need to rely on.
+@[expose] public section
+
+namespace Float.Model
+
+/--
+Compute the fused multiply-add `a * b + c` of three `Float.Model`, with a single rounding.
+-/
+def fma (a b c : Float.Model) : Float.Model :=
+  pack (UnpackedFloat.fma Format.binary64 a.unpack b.unpack c.unpack)
+
+end Float.Model
diff --git a/FFI/Model/Float32.lean b/FFI/Model/Float32.lean
new file mode 100644
index 0000000..04b6842
--- /dev/null
+++ b/FFI/Model/Float32.lean
@@ -0,0 +1,26 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+prelude
+public import Init.Data.Float.Model.Float32
+public import FFI.FMA
+
+-- This file is part of the logical model for floats which authors of float libraries
+-- need to rely on.
+@[expose] public section
+
+namespace Float32.Model
+
+open Float.Model (Format UnpackedFloat)
+
+/--
+Compute the fused multiply-add `a * b + c` of three `Float32.Model`, with a single rounding.
+-/
+def fma (a b c : Float32.Model) : Float32.Model :=
+  pack (UnpackedFloat.fma Format.binary32 a.unpack b.unpack c.unpack)
+
+end Float32.Model
diff --git a/RandomDo.lean b/RandomDo.lean
index e04384d..b2d4043 100644
--- a/RandomDo.lean
+++ b/RandomDo.lean
@@ -6,8 +6,20 @@ public import RandomDo.Monad.ForInInstances
 public import RandomDo.Monad.Instances
 public import RandomDo.Monad.MeasurableSpace
 public import RandomDo.Monad.Notation
-public import RandomDo.Tactic.Deriving
-public import RandomDo.Tactic.Elab
-public import RandomDo.Tactic.ForInStep
-public import RandomDo.Tactic.IsMarkov
-public import RandomDo.Tactic.Lemmas
+public import RandomDo.NumLean.Binomial
+public import RandomDo.NumLean.Distributions
+public import RandomDo.NumLean.PCG64
+public import RandomDo.NumLean.SeedSequence
+public import RandomDo.NumLean.Ziggurat
+public import RandomDo.NumLean.ZigguratSampler
+public import RandomDo.Tactic.Computable.Counterparts
+public import RandomDo.Tactic.Computable.Defs
+public import RandomDo.Tactic.Computable.Deriving
+public import RandomDo.Tactic.Computable.Example
+public import RandomDo.Tactic.Computable.Polymorphic.Polymorphic
+public import RandomDo.Tactic.Computable.Polymorphic.Scalar
+public import RandomDo.Tactic.IsMarkov.Defs
+public import RandomDo.Tactic.IsMarkov.Deriving
+public import RandomDo.Tactic.IsMarkov.Elab
+public import RandomDo.Tactic.IsMarkov.ForInStep
+public import RandomDo.Tactic.IsMarkov.Lemmas
diff --git a/RandomDo/NumLean/Binomial.lean b/RandomDo/NumLean/Binomial.lean
new file mode 100644
index 0000000..0f1d7cc
--- /dev/null
+++ b/RandomDo/NumLean/Binomial.lean
@@ -0,0 +1,194 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import RandomDo.NumLean.PCG64
+
+/-!
+# The binomial distribution
+
+`binomial n p` draws exactly what numpy's `Generator.binomial` draws: the inversion of the
+cumulative distribution where the mean `n * p` is at most `30`, the BTPE algorithm of
+Kachitvichyanukul and Schmeiser beyond, and the mirror image of either when `p > 1 / 2`.
+
+`binomial`, which chooses between the two and mirrors them, is in
+`RandomDo.NumLean.Distributions`, with the other distributions.
+
+## Main definitions
+
+* `Inversion`, `inversionSetup`, `inversionDraw`: numpy's `random_binomial_inversion`.
+* `Btpe`, `btpeSetup`, `Btpe.accept`, `btpeDraw`: numpy's `random_binomial_btpe`.
+
+## References
+
+* V. Kachitvichyanukul and B. W. Schmeiser, *Binomial random variate generation*, Communications
+  of the ACM 31 (1988), 216-222.
+* numpy's `numpy/random/src/distributions/distributions.c`.
+-/
+
+@[expose] public section
+
+namespace NumLean
+
+/-! ## Inversion -/
+
+/-- The constants the inversion reads a draw against. -/
+structure Inversion where
+  /-- The number of trials. -/
+  n : Float
+  /-- The probability of a success, at most one half. -/
+  p : Float
+  /-- The probability of a failure, `1 - p`. -/
+  q : Float
+  /-- `q ^ n`, the probability that no trial succeeds. -/
+  qn : Float
+  /-- The number of successes past which the walk gives up and starts over. -/
+  bound : Float
+
+/-- The constants of `random_binomial_inversion`, for `p ≤ 0.5` and `n * p ≤ 30`. -/
+@[inline] def inversionSetup (n p : Float) : Inversion :=
+  let q := 1.0 - p
+  let np := n * p
+  let b := np + 10.0 * Float.sqrt (np * q + 1)
+  { n, p, q, qn := Float.exp (n * Float.log q), bound := if n < b then n else b }
+
+/-- Walk up the cumulative distribution from zero until it passes `u`, as the loop of
+`random_binomial_inversion`. Answers `-1` where the walk runs past `bound`, which the tail beyond
+it is too thin to reach and where numpy starts the draw over. -/
+partial def inversionWalk (s : Inversion) (x px u : Float) : Float :=
+  if u > px then
+    let x := x + 1
+    if x > s.bound then -1
+    else inversionWalk s x (((s.n - x + 1) * s.p * px) / (x * s.q)) (u - px)
+  else x
+
+/-- Sample by inverting the cumulative distribution, as numpy's `random_binomial_inversion`. -/
+partial def inversionDraw (s : Inversion) : RandPCG IO Float := do
+  let x := inversionWalk s 0 s.qn (← random)
+  if x < 0 then inversionDraw s else return x
+
+/-! ## BTPE -/
+
+/-- The constants BTPE reads a draw against. -/
+structure Btpe where
+  /-- The number of trials. -/
+  n : Float
+  /-- The probability of a success, at most one half. -/
+  r : Float
+  /-- The probability of a failure, `1 - r`. -/
+  q : Float
+  /-- The mode of the distribution. -/
+  m : Float
+  /-- The middle of the triangle, `m + 1 / 2`. -/
+  xm : Float
+  /-- The half-width of the triangle. -/
+  p1 : Float
+  /-- The left end of the parallelogram. -/
+  xl : Float
+  /-- The right end of the parallelogram. -/
+  xr : Float
+  /-- The height of the parallelogram, relative to the triangle. -/
+  c : Float
+  /-- The rate of the left exponential tail. -/
+  laml : Float
+  /-- The rate of the right exponential tail. -/
+  lamr : Float
+  /-- The area of the triangle and the parallelogram. -/
+  p2 : Float
+  /-- The area of the triangle, the parallelogram and the left tail. -/
+  p3 : Float
+  /-- The area of all four regions, which a draw is scaled by. -/
+  p4 : Float
+  /-- The variance `n * r * q`. -/
+  nrq : Float
+
+/-- The constants of `random_binomial_btpe`, for `p ≤ 0.5` and `n * p > 30`. -/
+@[inline] def btpeSetup (n p : Float) : Btpe :=
+  let r := if p < 1.0 - p then p else 1.0 - p
+  let q := 1.0 - r
+  let fm := n * r + r
+  let m := Float.floor fm
+  let p1 := Float.floor (2.195 * Float.sqrt (n * r * q) - 4.6 * q) + 0.5
+  let xm := m + 0.5
+  let xl := xm - p1
+  let xr := xm + p1
+  let c := 0.134 + 20.5 / (15.3 + m)
+  let al := (fm - xl) / (fm - xl * r)
+  let ar := (xr - fm) / (xr * q)
+  let laml := al * (1.0 + al / 2.0)
+  let lamr := ar * (1.0 + ar / 2.0)
+  let p2 := p1 * (1.0 + 2.0 * c)
+  let p3 := p2 + c / laml
+  { n, r, q, m, xm, p1, xl, xr, c, laml, lamr, p2, p3, p4 := p3 + c / lamr, nrq := n * r * q }
+
+/-- One term of the Stirling series bounding `log` of a factorial, as the last test of BTPE spells
+it out. `u2` is `u * u`. -/
+@[inline] def btpeStirling (u u2 : Float) : Float :=
+  (13680.0 - (462.0 - (132.0 - (99.0 - 140.0 / u2) / u2) / u2) / u2) / u / 166320.0
+
+/-- The ratios of the probabilities from the mode up to `y`, multiplied into `f` one at a time as
+the step 50 of `random_binomial_btpe` takes them. -/
+partial def btpeUp (a s f i y : Float) : Float :=
+  if i ≤ y then btpeUp a s (f * (a / i - s)) (i + 1) y else f
+
+/-- The ratios from `y` up to the mode, divided out of `f` one at a time. Dividing the running
+value and dividing by the product do not round alike, and BTPE reads the first. -/
+partial def btpeDown (a s f i m : Float) : Float :=
+  if i ≤ m then btpeDown a s (f / (a / i - s)) (i + 1) m else f
+
+/-- Whether BTPE accepts the candidate `y` drawn with `v`, as the steps 50 and 52 of
+`random_binomial_btpe`: by the explicit product of the ratios of the probabilities between the mode
+and `y` when the two are close, and by a squeeze then the Stirling bound otherwise. -/
+def Btpe.accept (b : Btpe) (y v : Float) : Bool := Id.run do
+  let k := Float.abs (y - b.m)
+  unless k > 20 && k < b.nrq / 2.0 - 1 do
+    let s := b.r / b.q
+    let a := s * (b.n + 1)
+    if b.m < y then return !(v > btpeUp a s 1.0 (b.m + 1) y)
+    if b.m > y then return !(v > btpeDown a s 1.0 (y + 1) b.m)
+    return !(v > 1.0)
+  let rho := (k / b.nrq) * ((k * (k / 3.0 + 0.625) + 0.16666666666666666) / b.nrq + 0.5)
+  let t := -k * k / (2 * b.nrq)
+  let a := Float.log v
+  if a < t - rho then return true
+  if a > t + rho then return false
+  let x1 := y + 1
+  let f1 := b.m + 1
+  let z := b.n + 1 - b.m
+  let w := b.n - y + 1
+  return !(a > b.xm * Float.log (f1 / x1) + (b.n - b.m + 0.5) * Float.log (z / w)
+    + (y - b.m) * Float.log (w * b.r / (x1 * b.q))
+    + btpeStirling f1 (f1 * f1) + btpeStirling z (z * z)
+    + btpeStirling x1 (x1 * x1) + btpeStirling w (w * w))
+
+/-- Draw a candidate from the triangle, the parallelogram or one of the two exponential tails, and
+start over until one is accepted, as the steps 10 to 60 of `random_binomial_btpe`. -/
+partial def btpeDraw (b : Btpe) : RandPCG IO Float := do
+  let u := (← random) * b.p4
+  let v ← random
+  if u ≤ b.p1 then
+    return Float.floor (b.xm - b.p1 * v + u)
+  else if u ≤ b.p2 then
+    let x := b.xl + (u - b.p1) / b.c
+    let v := v * b.c + 1.0 - Float.abs (b.m - x + 0.5) / b.p1
+    if v > 1.0 then btpeDraw b else
+      let y := Float.floor x
+      if b.accept y v then return y else btpeDraw b
+  else if u ≤ b.p3 then
+    let y := Float.floor (b.xl + Float.log v / b.laml)
+    -- `v` can be zero, and the floor of the resulting infinity is no candidate.
+    if y < 0 || v == 0.0 then btpeDraw b else
+      let v := v * (u - b.p2) * b.laml
+      if b.accept y v then return y else btpeDraw b
+  else
+    let y := Float.floor (b.xr - Float.log v / b.lamr)
+    if y > b.n || v == 0.0 then btpeDraw b else
+      let v := v * (u - b.p3) * b.lamr
+      if b.accept y v then return y else btpeDraw b
+
+end NumLean
+
+end
diff --git a/RandomDo/NumLean/Distributions.lean b/RandomDo/NumLean/Distributions.lean
new file mode 100644
index 0000000..22dd979
--- /dev/null
+++ b/RandomDo/NumLean/Distributions.lean
@@ -0,0 +1,163 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import RandomDo.NumLean.PCG64
+public meta import RandomDo.NumLean.PCG64
+public import FFI.Float
+public import RandomDo.NumLean.Binomial
+public import RandomDo.NumLean.Ziggurat
+public import RandomDo.NumLean.ZigguratSampler
+
+/-!
+# Sample from specific distributions using the PCG-64 generator.
+
+This file provides samplers for specific distributions using the PCG-64 generator. Each one draws
+exactly what its numpy counterpart draws, so that a `RandPCG` program and a numpy `Generator`
+seeded alike produce the same values.
+
+The draws taken straight from the generator, `randUInt64`, `randUInt32` and `random`, are in
+`RandomDo.NumLean.PCG64`; the ziggurat the normal and the exponential share is in
+`RandomDo.NumLean.ZigguratSampler` and its tables in `RandomDo.NumLean.Ziggurat`.
+
+## Main definitions
+* `randInt`: sample an integer in `[low, high)` or `[low, high]`, as numpy's `Generator.integers`.
+* `uniform`: sample a `Float` in `[low, high)`, as numpy's `Generator.uniform`.
+* `standardNormal` / `normal`: sample a normal deviate, as numpy's `Generator.normal`.
+* `standardExponential` / `exponential`: sample an exponential deviate, as numpy's
+  `Generator.exponential`.
+-/
+
+@[expose] public section
+
+namespace NumLean
+
+/-- Sample a `UInt32` uniformly in `[0, rng]` by Lemire's nearly divisionless rejection method, as
+numpy's `buffered_bounded_lemire_uint32`. -/
+@[inline] def randLemireUInt32 (rng : UInt32) : RandPCG IO UInt32 := do
+  let rngExcl := rng + 1
+  let mut m := (← randUInt32).toUInt64 * rngExcl.toUInt64
+  if m.toUInt32 < rngExcl then
+    let threshold := (0xFFFFFFFF - rng) % rngExcl
+    while m.toUInt32 < threshold do
+      m := (← randUInt32).toUInt64 * rngExcl.toUInt64
+  return (m >>> 32).toUInt32
+
+/-- Sample a `UInt64` uniformly in `[0, rng]` by Lemire's method, as numpy's
+`bounded_lemire_uint64`. -/
+@[inline] def randLemireUInt64 (rng : UInt64) : RandPCG IO UInt64 := do
+  let rngExcl := rng + 1
+  let mut x ← randUInt64
+  let mut leftover := x * rngExcl
+  if leftover < rngExcl then
+    let threshold := (0xFFFFFFFFFFFFFFFF - rng) % rngExcl
+    while leftover < threshold do
+      x ← randUInt64
+      leftover := x * rngExcl
+  return PCG64.mulHi x rngExcl
+
+/-- Sample a `UInt64` uniformly in `[0, rng]`, as numpy's `random_bounded_uint64` with Lemire's
+rejection. -/
+@[inline] def randBoundedUInt64 (rng : UInt64) : RandPCG IO UInt64 := do
+  if rng == 0 then return 0
+  else if rng == 0xFFFFFFFF then return (← randUInt32).toUInt64
+  else if rng < 0xFFFFFFFF then return (← randLemireUInt32 rng.toUInt32).toUInt64
+  else if rng == 0xFFFFFFFFFFFFFFFF then randUInt64
+  else randLemireUInt64 rng
+
+/-- Sample an integer uniformly in `[low, high)`, or in `[low, high]` when `endpoint` is set, as
+numpy's `Generator.integers`. -/
+@[inline] def randInt₀ (low high : Int) (endpoint : Bool := false) : RandPCG IO Int := do
+  let high := if endpoint then high else high - 1
+  if low < Int64.minValue.toInt then throw <| IO.userError "low is out of bounds for int64"
+  if high > Int64.maxValue.toInt then throw <| IO.userError "high is out of bounds for int64"
+  if low > high then throw <| IO.userError (if endpoint then "low > high" else "low >= high")
+  let x ← randBoundedUInt64 (high - low).toNat.toUInt64
+  return low + x.toNat
+
+/-- Sample an integer uniformly in `[0, high)`, or in `[0, high]` when `endpoint` is set. -/
+@[inline] def randInt (high : Int) (endpoint : Bool := false) : RandPCG IO Int :=
+    randInt₀ 0 high endpoint
+
+/-- Sample an integer uniformly in `[low, high)`, or in `[low, high]` when `endpoint` is set. -/
+@[inline] def randBoundedInt (low high : Int) (endpoint : Bool := false) : RandPCG IO Int :=
+    randInt₀ low high endpoint
+
+/-- Sample a `Float` uniformly in `[low, high)`. -/
+@[inline] def uniform (low high : Float) : RandPCG IO Float := do
+  if low > high then throw <| IO.userError "low > high"
+  let x ← random
+  return Float.fma x (high - low) low
+
+/-- Sample the tail of the standard normal beyond `Ziggurat.norR`, as the `idx == 0` branch of
+numpy's `random_standard_normal`: draw from an exponential tail until the pair of draws falls under
+the normal's, which is Marsaglia's method for the tail. `negate` carries the sign numpy reads off
+the integer already drawn, which the loop does not redraw. -/
+partial def normalTail (negate : Bool) : RandPCG IO Float := do
+  let xx := -Ziggurat.norInvR * Float.log1p (-(← random))
+  let yy := -Float.log1p (-(← random))
+  if yy + yy > xx * xx then
+    return if negate then -(Ziggurat.norR + xx) else Ziggurat.norR + xx
+  normalTail negate
+
+/-- Sample from the standard normal, as numpy's `random_standard_normal`. The 64-bit output gives
+the strip in its low byte, then the sign, then a 52-bit abscissa; the tail takes its sign from a
+further bit of that same abscissa, as numpy does. -/
+@[inline] def standardNormal : RandPCG IO Float :=
+  ziggurat Ziggurat.ki Ziggurat.wi Ziggurat.fi
+    (fun r =>
+      let idx := (r &&& 0xFF).toNat
+      let r := r >>> 8
+      (idx, (r >>> 1) &&& 0x000FFFFFFFFFFFFF, (r &&& 1) == 1))
+    (fun x => Float.exp ((-0.5) * x * x))
+    (fun rabs => normalTail (((rabs >>> 8) &&& 1) == 1))
+
+/-- Sample from the standard exponential, as numpy's `random_standard_exponential`. The 64-bit
+output is first shifted by three, then gives the strip in its low byte and a 53-bit abscissa; the
+tail is the exponential's own, memoryless, so one draw beyond `Ziggurat.expR` suffices. -/
+@[inline] def standardExponential : RandPCG IO Float :=
+  ziggurat Ziggurat.ke Ziggurat.we Ziggurat.fe
+    (fun r =>
+      let r := r >>> 3
+      ((r &&& 0xFF).toNat, r >>> 8, false))
+    (fun x => Float.exp (-x))
+    (fun _ => do return Ziggurat.expR - Float.log1p (-(← random)))
+
+/-- Draw random samples from a normal (Gaussian) distribution. -/
+@[inline] def normal (loc : Float := 0) (scale : Float := 1) : RandPCG IO Float := do
+  if scale < 0 then throw <| IO.userError "scale < 0"
+  return Float.fma scale (← standardNormal) loc
+
+/-- Draw random samples from a normal (Gaussian) distribution, with variance instead of standard
+deviation. -/
+@[inline] def normal' (loc : Float := 0) (var : Float := 1) : RandPCG IO Float := do
+  if var < 0 then throw <| IO.userError "var < 0"
+  normal loc (Float.sqrt var)
+
+/-- Draw samples from an exponential distribution. -/
+@[inline] def exponential (scale : Float := 1) : RandPCG IO Float := do
+  if scale < 0 then throw <| IO.userError "scale < 0"
+  return scale * (← standardExponential)
+
+/-- Draw samples from a binomial distribution. -/
+def binomial (n : Nat) (p : Float) : RandPCG IO Nat := do
+  -- The comparisons are the `Bool` ones: through `Decidable`, each costs more than a draw.
+  if p < 0.0 || 1.0 < p || p.isNaN then
+    throw <| IO.userError s!"expected 0 ≤ p ≤ 1, got {p}"
+  let n := n.toUInt64.toFloat
+  if n == 0 || p == 0.0 then return 0
+  if Float.le p 0.5 then
+    if Float.le (p * n) 30.0 then return (← inversionDraw (inversionSetup n p)).toUInt64.toNat
+    else return (← btpeDraw (btpeSetup n p)).toUInt64.toNat
+  else
+    let q := 1.0 - p
+    if Float.le (q * n) 30.0 then return (n - (← inversionDraw (inversionSetup n q))).toUInt64.toNat
+    else return (n - (← btpeDraw (btpeSetup n q))).toUInt64.toNat
+
+/-- Draw samples from a Bernoulli distribution. -/
+@[inline] def bernoulli (p : Float) := binomial 1 p
+
+end NumLean
diff --git a/RandomDo/NumLean/PCG64.lean b/RandomDo/NumLean/PCG64.lean
new file mode 100644
index 0000000..6a0a2e3
--- /dev/null
+++ b/RandomDo/NumLean/PCG64.lean
@@ -0,0 +1,251 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import Mathlib.Control.Random
+public import RandomDo.NumLean.SeedSequence
+public import Batteries.Lean.LawfulMonad
+
+/-!
+# The PCG-64 pseudo-random number generator
+
+This file provides `PCG64`, a drop-in replacement for `StdGen` implementing the
+`PCG-XSL-RR 128/64` variant of the permuted congruential generator family of O'Neill, i.e. the
+generator known as `pcg64` in the reference C implementation and as `numpy.random.PCG64` in numpy.
+
+The generator is a linear congruential generator on 128 bits,
+`state ← state * multiplier + increment`, whose state is scrambled down to 64 bits by the
+`XSL-RR` output permutation (xor the two halves together, then rotate by the top 6 bits of the
+state). The 128-bit arithmetic is emulated with pairs of `UInt64`, so every operation compiles to
+native machine arithmetic.
+
+## Main definitions
+
+* `PCG64`: the generator state, and its `RandomGen` instance
+* `PCG64.nextUInt64` / `PCG64.nextUInt32`: the 64-bit output, and numpy's 32-bit output, which
+  hands out the two halves of a 64-bit output in turn
+* `PCG64.seedWords` / `PCG64.seed` / `mkPCG64`: seeding, following the reference
+  `pcg_setseq_128_srandom_r`
+* `randUInt64` / `randUInt32` / `random`: the three outputs a `RandPCG` computation draws straight
+  from the generator.
+
+## References
+
+* M. E. O'Neill, *PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for
+  Random Number Generation*, 2014. <https://www.pcg-random.org/paper.html>
+* The reference implementation: <https://github.com/imneme/pcg-c>
+* numpy's vendored copy: `numpy/random/src/pcg64/pcg64.h`
+-/
+
+@[expose] public section
+
+namespace NumLean
+
+/-- The state of a PCG-64 generator: a 128-bit LCG state and a 128-bit increment (which selects the
+stream and must be odd), each stored as a pair of 64-bit words, together with the one-word buffer
+numpy's bit generator keeps for its 32-bit outputs, see `nextUInt32`. -/
+structure PCG64 where
+  /-- High 64 bits of the LCG state. -/
+  stateHi : UInt64
+  /-- Low 64 bits of the LCG state. -/
+  stateLo : UInt64
+  /-- High 64 bits of the increment. -/
+  incHi : UInt64
+  /-- Low 64 bits of the increment; it is odd for any generator built through the API below. -/
+  incLo : UInt64
+  /-- Whether `uinteger` holds the unused upper half of the last output drawn through `nextUInt32`.
+  Mirrors the `has_uint32` field of numpy's `pcg64_state`. -/
+  hasUInt32 : Bool := false
+  /-- The upper 32 bits of the last output drawn through `nextUInt32`, meaningful only while
+  `hasUInt32` is set. Mirrors the `uinteger` field of numpy's `pcg64_state`. -/
+  uinteger : UInt32 := 0
+  deriving Repr, DecidableEq
+
+namespace PCG64
+
+/-- High 64 bits of the PCG-64 multiplier `0x2360ED051FC65DA44385DF649FCCF645`. -/
+def multHi : UInt64 := 0x2360ED051FC65DA4
+
+/-- Low 64 bits of the PCG-64 multiplier `0x2360ED051FC65DA44385DF649FCCF645`. -/
+def multLo : UInt64 := 0x4385DF649FCCF645
+
+/-- The default stream, i.e. the reference default increment `0x5851F42D4C957F2D14057B7EF767814F`
+divided by two, since `seed` turns a stream `s` into the increment `2 * s + 1`. -/
+def defaultStream : Nat := 0x2C28FA16A64ABF968A02BDBF7BB3C0A7
+
+/-- The high 64 bits of the 128-bit product `a * b`, obtained from the four 32-bit limb products. -/
+@[inline] def mulHi (a b : UInt64) : UInt64 :=
+  let mask : UInt64 := 0xFFFFFFFF
+  let a₀ := a &&& mask
+  let a₁ := a >>> 32
+  let b₀ := b &&& mask
+  let b₁ := b >>> 32
+  let t := a₁ * b₀ + (a₀ * b₀) >>> 32
+  a₁ * b₁ + t >>> 32 + (a₀ * b₁ + (t &&& mask)) >>> 32
+
+/-- Rotate the 64-bit word `x` right by `r` bits. Only the low 6 bits of `r` are used. -/
+@[inline] def rotr (x r : UInt64) : UInt64 := (x >>> r) ||| (x <<< (64 - r))
+
+/-- One step of the underlying 128-bit LCG, `state ← state * multiplier + increment`. -/
+@[inline] def step (g : PCG64) : PCG64 :=
+  let lo := g.stateLo * multLo
+  let hi := mulHi g.stateLo multLo + g.stateLo * multHi + g.stateHi * multLo
+  let lo' := lo + g.incLo
+  -- unsigned addition wraps, so `lo' < lo` exactly when the low half carried
+  let carry : UInt64 := if lo' < lo then 1 else 0
+  { g with stateHi := hi + g.incHi + carry, stateLo := lo' }
+
+/-- The `XSL-RR` output permutation: fold the 128-bit state onto 64 bits by xoring its two halves,
+then rotate the result right by the 6 most significant bits of the state. -/
+@[inline] def output (g : PCG64) : UInt64 :=
+  rotr (g.stateHi ^^^ g.stateLo) (g.stateHi >>> 58)
+
+/-- The next 64-bit output, together with the advanced generator. As in the reference
+implementation, the state is stepped before the output permutation is applied. -/
+@[inline] def nextUInt64 (g : PCG64) : UInt64 × PCG64 :=
+  let g := g.step
+  (g.output, g)
+
+/-- The next 32-bit output, together with the advanced generator, as numpy's `pcg64_next32`: a
+64-bit output is drawn and its low half returned, its high half being kept in `uinteger` to be
+returned by the next call, which then does not step the generator. The buffer survives
+intermediate `nextUInt64` calls, as it does in numpy. -/
+@[inline] def nextUInt32 (g : PCG64) : UInt32 × PCG64 :=
+  if g.hasUInt32 then
+    (g.uinteger, { g with hasUInt32 := false })
+  else
+    let (x, g) := g.nextUInt64
+    (x.toUInt32, { g with hasUInt32 := true, uinteger := (x >>> 32).toUInt32 })
+
+/-- Seed a generator from four 64-bit words, following the reference
+`pcg_setseq_128_srandom_r`: `stateHi:stateLo` is the initial state and `seqHi:seqLo` selects the
+stream, whose increment is `2 * seq + 1`. -/
+def seedWords (stateHi stateLo seqHi seqLo : UInt64) : PCG64 :=
+  let g : PCG64 :=
+    { stateHi := 0, stateLo := 0,
+      incHi := (seqHi <<< 1) ||| (seqLo >>> 63), incLo := (seqLo <<< 1) ||| 1 }
+  let g := g.step
+  let lo := g.stateLo + stateLo
+  let carry : UInt64 := if lo < g.stateLo then 1 else 0
+  step { g with stateHi := g.stateHi + stateHi + carry, stateLo := lo }
+
+/-- Seed a generator from a `SeedSequence`, as numpy's `PCG64` constructor does: the first two
+words drawn from the mixer give the initial state and the next two the stream. -/
+def ofSeedSequence (s : SeedSequence) : PCG64 :=
+  let w := s.generateState 4
+  seedWords w[0]! w[1]! w[2]! w[3]!
+
+/-- Seed a generator from an initial state. -/
+def seed (n : Nat) : PCG64 := ofSeedSequence (SeedSequence.ofNat n)
+
+/-- The `SeedSequence` whose entropy is `n` outputs taken from `g`, split into 32-bit words, the
+advanced generator being returned alongside. Used to derive further generators from an existing
+one. Two outputs are enough by default: they fill the mixing pool exactly, and no amount of extra
+entropy would make it carry more than its `SeedSequence.poolSize` words. -/
+def toSeedSequence (g : PCG64) (n : Nat := 2) : SeedSequence × PCG64 :=
+  let (words, g) := Id.run do
+    let mut g := g
+    let mut words := Array.emptyWithCapacity (2 * n)
+    for _ in List.range n do
+      let (x, g') := g.nextUInt64
+      g := g'
+      words := (words.push x.toUInt32).push (x >>> 32).toUInt32
+    return (words, g)
+  (SeedSequence.ofWords words #[], g)
+
+/-- The range of values returned by `PCG64`, namely all of `[0, 2 ^ 64 - 1]`. -/
+def range : Nat × Nat := (0, UInt64.size - 1)
+
+/-- Derive `n` generators from one, together with the parent advanced past the outputs used as
+entropy. The children are the `SeedSequence.spawn` children of the entropy drawn from the parent,
+so they are obtained exactly as any other family of generators in this library. -/
+def spawn (g : PCG64) (n : Nat) : Array PCG64 × PCG64 :=
+  let (s, g) := g.toSeedSequence
+  ((s.spawn n).1.map ofSeedSequence, g)
+
+/-- Derive two generators from one: the first is the current one advanced by two steps, the second
+is the first `SeedSequence.spawn` child of those two outputs. Splitting is not part of the PCG
+specification and nothing here establishes that the two streams are independent; going through the
+mixer only removes the direct algebraic tie between the child's initial state and two consecutive
+states of the parent. -/
+def split (g : PCG64) : PCG64 × PCG64 :=
+  let (s, g) := g.toSeedSequence
+  (g, ofSeedSequence (s.spawn 1).1[0]!)
+
+/-- The first `n` outputs of `g`. -/
+def take (g : PCG64) (n : Nat) : Array UInt64 := Id.run do
+  let mut g := g
+  let mut out := Array.emptyWithCapacity n
+  for _ in [:n] do
+    let (x, g') := g.nextUInt64
+    g := g'
+    out := out.push x
+  return out
+
+end PCG64
+
+/-- Returns a PCG-64 generator seeded with `s`, on the default stream. The analogue of
+`mkStdGen`. -/
+def mkPCG64 (s : Nat := 0) : PCG64 := PCG64.seed s
+
+instance : Inhabited PCG64 := ⟨mkPCG64⟩
+
+instance : RandomGen PCG64 where
+  range _ := PCG64.range
+  next g := let (x, g) := g.nextUInt64; (x.toNat, g)
+  split := PCG64.split
+
+/-- A monad transformer to generate random objects using the generator type `PCG64`.
+`RandPCG m α` should be thought of a random value in `m α`. -/
+abbrev RandPCG := RandGT PCG64
+
+/-- Sample a `UInt64` from a PCG-64 generator. -/
+@[inline] def randUInt64 : RandPCG IO UInt64 := do
+  let (x, g) := (← get).down.nextUInt64
+  set (ULift.up g)
+  return x
+
+/-- Sample a `UInt32` from a PCG-64 generator, as numpy's `next_uint32` does for `PCG64`: the two
+halves of each 64-bit output are handed out in turn, see `PCG64.nextUInt32`. -/
+@[inline] def randUInt32 : RandPCG IO UInt32 := do
+  let (x, g) := (← get).down.nextUInt32
+  set (ULift.up g)
+  return x
+
+/-- Sample a `Float` in `[0, 1)` from a PCG-64 generator, as numpy's `next_double`: the top 53 bits
+of a 64-bit output, scaled by `2 ^ (-53)`. -/
+@[inline] def random : RandPCG IO Float := do
+  let x ← randUInt64
+  return (x >>> 11).toFloat * (Float.ofBits <| 0x3CA <<< (52 : UInt64))
+
+end NumLean
+
+namespace IO
+
+open NumLean
+
+/-- A global reference to a PCG-64 generator, seeded from the system's random source. -/
+initialize PCG64Ref : Ref PCG64 ←
+  let seed := UInt64.toNat (ByteArray.toUInt64LE! (← IO.getRandomBytes 8))
+  IO.mkRef (mkPCG64 seed)
+
+variable {m : Type* → Type*} {m₀ : Type → Type}
+variable [Monad m] [MonadLiftT (ST RealWorld) m₀] [ULiftable m₀ m]
+
+set_option autoImplicit true
+
+/-- Execute `RandPCG m α` using the global `PCG64Ref` as RNG. -/
+def runRandPCG (cmd : RandPCG m α) : m α := do
+  let PCG64 ← ULiftable.up (PCG64Ref.get : m₀ _)
+  let (res, new) ← StateT.run cmd PCG64
+  let _ ← ULiftable.up (PCG64Ref.set new.down : m₀ _)
+  pure res
+
+/-- Execute `RandPCG m α` using the global `PCG64Ref` as RNG and the given `seed`. -/
+def runRandPCGWith (seed : Nat) (cmd : RandPCG m α) : m α := do
+  pure <| (← cmd.run (ULift.up <| mkPCG64 seed)).1
+
+end IO
diff --git a/RandomDo/NumLean/SeedSequence.lean b/RandomDo/NumLean/SeedSequence.lean
new file mode 100644
index 0000000..b487cf1
--- /dev/null
+++ b/RandomDo/NumLean/SeedSequence.lean
@@ -0,0 +1,161 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+/-!
+# Numpy's `SeedSequence` entropy mixer
+
+Numpy instead runs the user's seed through `SeedSequence`, an entropy mixer whose job is to turn a
+small, structured seed (`0`, `1`, `2`, ...) into 128 well-distributed bits, so that consecutive
+seeds yield unrelated streams. The algorithm is a two-stage avalanche on 32-bit words: a pool of
+four words absorbs the entropy, every word is mixed into every other one, then the requested output
+words are drawn from the pool through a second hash.
+
+## Main definitions
+
+* `SeedSequence`: the mixer state, built by `SeedSequence.ofNat`;
+* `SeedSequence.generateState`: draw 64-bit words from the pool;
+* `SeedSequence.spawn`: derive independent child sequences, Numpy's principled alternative to
+  `RandomGen.split`;
+
+## References
+
+* Numpy's implementation: `numpy/random/bit_generator.pyx`.
+* The design it follows: M. E. O'Neill, *Developing a seed_seq Alternative*, 2015.
+  <https://www.pcg-random.org/posts/developing-a-seed_seq-alternative.html>
+-/
+
+@[expose] public section
+
+namespace NumLean
+
+/-- Numpy's entropy mixer. `entropy` is the user's seed and `spawnKey` identifies the sequence among
+the descendants of the original one; the `pool` is the mixed entropy both are absorbed into. -/
+structure SeedSequence where
+  /-- The seed, as 32-bit words, least significant first. -/
+  entropy : Array UInt32
+  /-- The path identifying this sequence among the descendants of the root one. -/
+  spawnKey : Array UInt32
+  /-- How many children have already been spawned from this sequence. -/
+  nChildrenSpawned : Nat
+  /-- The mixed entropy, of length `SeedSequence.poolSize`. -/
+  pool : Array UInt32
+  deriving Repr, Inhabited
+
+namespace SeedSequence
+
+/-- Number of 32-bit words held in the mixing pool. -/
+def poolSize : Nat := 4
+
+/-- Amount by which a hashed word is xor-shifted, half of a word's width. -/
+def xshift : UInt32 := 16
+
+/-- Initial hash constant of the entropy-absorbing stage. -/
+def initA : UInt32 := 0x43B0D7E5
+
+/-- Multiplier of the entropy-absorbing stage. -/
+def multA : UInt32 := 0x931E8875
+
+/-- Initial hash constant of the output stage. -/
+def initB : UInt32 := 0x8B51F9DD
+
+/-- Multiplier of the output stage. -/
+def multB : UInt32 := 0x58F38DED
+
+/-- Left multiplier of the pool mixing function. -/
+def mixMultL : UInt32 := 0xCA01F9DD
+
+/-- Right multiplier of the pool mixing function. -/
+def mixMultR : UInt32 := 0x4973F715
+
+/-- Hash a word against the running hash constant, returning the hashed word and the advanced
+constant. -/
+@[inline] def hashmix (value hashConst : UInt32) : UInt32 × UInt32 :=
+  let value := value ^^^ hashConst
+  let hashConst := hashConst * multA
+  let value := value * hashConst
+  (value ^^^ (value >>> xshift), hashConst)
+
+/-- Combine two pool words. -/
+@[inline] def mix (x y : UInt32) : UInt32 :=
+  let r := mixMultL * x - mixMultR * y
+  r ^^^ (r >>> xshift)
+
+/-- Auxiliary function for `toWords`, accumulating the 32-bit words of `n`. -/
+def toWordsAux (n : Nat) (acc : Array UInt32) : Array UInt32 :=
+  if h : n = 0 then acc
+  else toWordsAux (n / 4294967296) (acc.push (UInt32.ofNat (n % 4294967296)))
+termination_by n
+decreasing_by exact Nat.div_lt_self (Nat.pos_of_ne_zero h) (by omega)
+
+/-- Split a natural number into 32-bit words, least significant first. Zero maps to `#[0]`, as in
+Numpy's `_int_to_uint32_array`. -/
+def toWords (n : Nat) : Array UInt32 := if n = 0 then #[0] else toWordsAux n #[]
+
+/-- Absorb an entropy array into a pool of `poolSize` words: seed the pool, then mix every word
+into every other one so that late bits influence early ones, and finally fold in whatever entropy
+did not fit in the pool. -/
+def mixEntropy (entropy : Array UInt32) : Array UInt32 := Id.run do
+  let mut pool : Array UInt32 := Array.replicate poolSize 0
+  let mut hashConst := initA
+  for i in [:poolSize] do
+    let (v, hc) := hashmix (entropy[i]?.getD 0) hashConst
+    pool := pool.set! i v
+    hashConst := hc
+  for src in [:poolSize] do
+    for dst in [:poolSize] do
+      if src ≠ dst then
+        let (h, hc) := hashmix pool[src]! hashConst
+        pool := pool.set! dst (mix pool[dst]! h)
+        hashConst := hc
+  for src in [poolSize:entropy.size] do
+    for dst in [:poolSize] do
+      let (h, hc) := hashmix entropy[src]! hashConst
+      pool := pool.set! dst (mix pool[dst]! h)
+      hashConst := hc
+  return pool
+
+/-- Build a sequence from an entropy array and a spawn key. When a spawn key is present, the
+entropy is padded with zeros up to the pool size, so that a child cannot collide with a root
+sequence whose seed happens to look like the concatenation of the two. -/
+def ofWords (entropy spawnKey : Array UInt32) : SeedSequence :=
+  let entropy' :=
+    if spawnKey.isEmpty || poolSize ≤ entropy.size then entropy
+    else entropy ++ Array.replicate (poolSize - entropy.size) 0
+  { entropy, spawnKey, nChildrenSpawned := 0, pool := mixEntropy (entropy' ++ spawnKey) }
+
+/-- The sequence Numpy builds from an integer seed, i.e. `numpy.random.SeedSequence(n)`. -/
+def ofNat (n : Nat) : SeedSequence := ofWords (toWords n) #[]
+
+/-- Draw `n` 64-bit words from the pool. Each is assembled from two consecutive 32-bit outputs,
+least significant first, as Numpy does when asked for `dtype=np.uint64`. -/
+def generateState (s : SeedSequence) (n : Nat) : Array UInt64 := Id.run do
+  let mut words : Array UInt32 := Array.emptyWithCapacity (2 * n)
+  let mut hashConst := initB
+  for i in [:2 * n] do
+    let v := s.pool[i % poolSize]! ^^^ hashConst
+    hashConst := hashConst * multB
+    let v := v * hashConst
+    words := words.push (v ^^^ (v >>> xshift))
+  let mut out := Array.emptyWithCapacity n
+  for i in [:n] do
+    out := out.push (words[2 * i]!.toUInt64 ||| (words[2 * i + 1]!.toUInt64 <<< 32))
+  return out
+
+/-- Derive `n` child sequences, together with the parent updated to remember how many children it
+has produced. Each child absorbs a distinct spawn key into the mixer, so it is a pure function of
+the root entropy and of its key path: it can be rebuilt directly, without replaying any draw, and
+does not depend on how much the parent has been used. Their statistical independence rests on the
+avalanche quality of the mixer, not on a proof. -/
+def spawn (s : SeedSequence) (n : Nat) : Array SeedSequence × SeedSequence := Id.run do
+  let mut children := Array.emptyWithCapacity n
+  for i in [:n] do
+    children := children.push (ofWords s.entropy (s.spawnKey ++ toWords (s.nChildrenSpawned + i)))
+  return (children, { s with nChildrenSpawned := s.nChildrenSpawned + n })
+
+end SeedSequence
+
+end NumLean
diff --git a/RandomDo/NumLean/Ziggurat.lean b/RandomDo/NumLean/Ziggurat.lean
new file mode 100644
index 0000000..89dce60
--- /dev/null
+++ b/RandomDo/NumLean/Ziggurat.lean
@@ -0,0 +1,527 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+/-!
+# The ziggurat tables
+
+The tables numpy's `random_standard_normal` and `random_standard_exponential` read, copied verbatim
+from its `numpy/random/src/distributions/ziggurat_constants.h`. Lean parses each of these literals
+to the same `Float` the C compiler produces, so the tables are the same numbers numpy uses.
+
+A ziggurat covers the density by 256 horizontal strips of equal area. `w` scales a drawn integer to
+an abscissa in its strip, `k` is the largest such integer that falls in the rectangular part of the
+strip, where the draw is accepted at once, and `f` is the density at the strip's edge, used by the
+rejection test on the strips that stick out of the curve. The `i` tables serve the normal and the
+`e` tables the exponential.
+
+The single-precision tables are not copied: no sampler here needs them yet.
+-/
+
+@[expose] public section
+
+namespace NumLean.Ziggurat
+
+/-- The abscissa of the normal's base strip, beyond which its tail is sampled separately. -/
+def norR : Float := 3.6541528853610087963519472518
+
+/-- `1 / norR`. -/
+def norInvR : Float := 0.27366123732975827203338247596
+
+/-- The abscissa of the exponential's base strip, beyond which its tail is sampled separately. -/
+def expR : Float := 7.6971174701310497140446280481
+
+/-- For each strip of the normal, the largest 52-bit draw landing in its rectangular part. -/
+def ki : Array UInt64 := #[
+  0x000EF33D8025EF6A, 0x0000000000000000, 0x000C08BE98FBC6A8, 0x000DA354FABD8142,
+  0x000E51F67EC1EEEA, 0x000EB255E9D3F77E, 0x000EEF4B817ECAB9, 0x000F19470AFA44AA,
+  0x000F37ED61FFCB18, 0x000F4F469561255C, 0x000F61A5E41BA396, 0x000F707A755396A4,
+  0x000F7CB2EC28449A, 0x000F86F10C6357D3, 0x000F8FA6578325DE, 0x000F9724C74DD0DA,
+  0x000F9DA907DBF509, 0x000FA360F581FA74, 0x000FA86FDE5B4BF8, 0x000FACF160D354DC,
+  0x000FB0FB6718B90F, 0x000FB49F8D5374C6, 0x000FB7EC2366FE77, 0x000FBAECE9A1E50E,
+  0x000FBDAB9D040BED, 0x000FC03060FF6C57, 0x000FC2821037A248, 0x000FC4A67AE25BD1,
+  0x000FC6A2977AEE31, 0x000FC87AA92896A4, 0x000FCA325E4BDE85, 0x000FCBCCE902231A,
+  0x000FCD4D12F839C4, 0x000FCEB54D8FEC99, 0x000FD007BF1DC930, 0x000FD1464DD6C4E6,
+  0x000FD272A8E2F450, 0x000FD38E4FF0C91E, 0x000FD49A9990B478, 0x000FD598B8920F53,
+  0x000FD689C08E99EC, 0x000FD76EA9C8E832, 0x000FD848547B08E8, 0x000FD9178BAD2C8C,
+  0x000FD9DD07A7ADD2, 0x000FDA9970105E8C, 0x000FDB4D5DC02E20, 0x000FDBF95C5BFCD0,
+  0x000FDC9DEBB99A7D, 0x000FDD3B8118729D, 0x000FDDD288342F90, 0x000FDE6364369F64,
+  0x000FDEEE708D514E, 0x000FDF7401A6B42E, 0x000FDFF46599ED40, 0x000FE06FE4BC24F2,
+  0x000FE0E6C225A258, 0x000FE1593C28B84C, 0x000FE1C78CBC3F99, 0x000FE231E9DB1CAA,
+  0x000FE29885DA1B91, 0x000FE2FB8FB54186, 0x000FE35B33558D4A, 0x000FE3B799D0002A,
+  0x000FE410E99EAD7F, 0x000FE46746D47734, 0x000FE4BAD34C095C, 0x000FE50BAED29524,
+  0x000FE559F74EBC78, 0x000FE5A5C8E41212, 0x000FE5EF3E138689, 0x000FE6366FD91078,
+  0x000FE67B75C6D578, 0x000FE6BE661E11AA, 0x000FE6FF55E5F4F2, 0x000FE73E5900A702,
+  0x000FE77B823E9E39, 0x000FE7B6E37070A2, 0x000FE7F08D774243, 0x000FE8289053F08C,
+  0x000FE85EFB35173A, 0x000FE893DC840864, 0x000FE8C741F0CEBC, 0x000FE8F9387D4EF6,
+  0x000FE929CC879B1D, 0x000FE95909D388EA, 0x000FE986FB939AA2, 0x000FE9B3AC714866,
+  0x000FE9DF2694B6D5, 0x000FEA0973ABE67C, 0x000FEA329CF166A4, 0x000FEA5AAB32952C,
+  0x000FEA81A6D5741A, 0x000FEAA797DE1CF0, 0x000FEACC85F3D920, 0x000FEAF07865E63C,
+  0x000FEB13762FEC13, 0x000FEB3585FE2A4A, 0x000FEB56AE3162B4, 0x000FEB76F4E284FA,
+  0x000FEB965FE62014, 0x000FEBB4F4CF9D7C, 0x000FEBD2B8F449D0, 0x000FEBEFB16E2E3E,
+  0x000FEC0BE31EBDE8, 0x000FEC2752B15A15, 0x000FEC42049DAFD3, 0x000FEC5BFD29F196,
+  0x000FEC75406CEEF4, 0x000FEC8DD2500CB4, 0x000FECA5B6911F12, 0x000FECBCF0C427FE,
+  0x000FECD38454FB15, 0x000FECE97488C8B3, 0x000FECFEC47F91B7, 0x000FED1377358528,
+  0x000FED278F844903, 0x000FED3B10242F4C, 0x000FED4DFBAD586E, 0x000FED605498C3DD,
+  0x000FED721D414FE8, 0x000FED8357E4A982, 0x000FED9406A42CC8, 0x000FEDA42B85B704,
+  0x000FEDB3C8746AB4, 0x000FEDC2DF416652, 0x000FEDD171A46E52, 0x000FEDDF813C8AD3,
+  0x000FEDED0F909980, 0x000FEDFA1E0FD414, 0x000FEE06AE124BC4, 0x000FEE12C0D95A06,
+  0x000FEE1E579006E0, 0x000FEE29734B6524, 0x000FEE34150AE4BC, 0x000FEE3E3DB89B3C,
+  0x000FEE47EE2982F4, 0x000FEE51271DB086, 0x000FEE59E9407F41, 0x000FEE623528B42E,
+  0x000FEE6A0B5897F1, 0x000FEE716C3E077A, 0x000FEE7858327B82, 0x000FEE7ECF7B06BA,
+  0x000FEE84D2484AB2, 0x000FEE8A60B66343, 0x000FEE8F7ACCC851, 0x000FEE94207E25DA,
+  0x000FEE9851A829EA, 0x000FEE9C0E13485C, 0x000FEE9F557273F4, 0x000FEEA22762CCAE,
+  0x000FEEA4836B42AC, 0x000FEEA668FC2D71, 0x000FEEA7D76ED6FA, 0x000FEEA8CE04FA0A,
+  0x000FEEA94BE8333B, 0x000FEEA950296410, 0x000FEEA8D9C0075E, 0x000FEEA7E7897654,
+  0x000FEEA678481D24, 0x000FEEA48AA29E83, 0x000FEEA21D22E4DA, 0x000FEE9F2E352024,
+  0x000FEE9BBC26AF2E, 0x000FEE97C524F2E4, 0x000FEE93473C0A3A, 0x000FEE8E40557516,
+  0x000FEE88AE369C7A, 0x000FEE828E7F3DFD, 0x000FEE7BDEA7B888, 0x000FEE749BFF37FF,
+  0x000FEE6CC3A9BD5E, 0x000FEE64529E007E, 0x000FEE5B45A32888, 0x000FEE51994E57B6,
+  0x000FEE474A0006CF, 0x000FEE3C53E12C50, 0x000FEE30B2E02AD8, 0x000FEE2462AD8205,
+  0x000FEE175EB83C5A, 0x000FEE09A22A1447, 0x000FEDFB27E349CC, 0x000FEDEBEA76216C,
+  0x000FEDDBE422047E, 0x000FEDCB0ECE39D3, 0x000FEDB964042CF4, 0x000FEDA6DCE938C9,
+  0x000FED937237E98D, 0x000FED7F1C38A836, 0x000FED69D2B9C02B, 0x000FED538D06AE00,
+  0x000FED3C41DEA422, 0x000FED23E76A2FD8, 0x000FED0A732FE644, 0x000FECEFDA07FE34,
+  0x000FECD4100EB7B8, 0x000FECB708956EB4, 0x000FEC98B61230C1, 0x000FEC790A0DA978,
+  0x000FEC57F50F31FE, 0x000FEC356686C962, 0x000FEC114CB4B335, 0x000FEBEB948E6FD0,
+  0x000FEBC429A0B692, 0x000FEB9AF5EE0CDC, 0x000FEB6FE1C98542, 0x000FEB42D3AD1F9E,
+  0x000FEB13B00B2D4B, 0x000FEAE2591A02E9, 0x000FEAAEAE992257, 0x000FEA788D8EE326,
+  0x000FEA3FCFFD73E5, 0x000FEA044C8DD9F6, 0x000FE9C5D62F563B, 0x000FE9843BA947A4,
+  0x000FE93F471D4728, 0x000FE8F6BD76C5D6, 0x000FE8AA5DC4E8E6, 0x000FE859E07AB1EA,
+  0x000FE804F690A940, 0x000FE7AB488233C0, 0x000FE74C751F6AA5, 0x000FE6E8102AA202,
+  0x000FE67DA0B6ABD8, 0x000FE60C9F38307E, 0x000FE5947338F742, 0x000FE51470977280,
+  0x000FE48BD436F458, 0x000FE3F9BFFD1E37, 0x000FE35D35EEB19C, 0x000FE2B5122FE4FE,
+  0x000FE20003995557, 0x000FE13C82788314, 0x000FE068C4EE67B0, 0x000FDF82B02B71AA,
+  0x000FDE87C57EFEAA, 0x000FDD7509C63BFD, 0x000FDC46E529BF13, 0x000FDAF8F82E0282,
+  0x000FD985E1B2BA75, 0x000FD7E6EF48CF04, 0x000FD613ADBD650B, 0x000FD40149E2F012,
+  0x000FD1A1A7B4C7AC, 0x000FCEE204761F9E, 0x000FCBA8D85E11B2, 0x000FC7D26ECD2D22,
+  0x000FC32B2F1E22ED, 0x000FBD6581C0B83A, 0x000FB606C4005434, 0x000FAC40582A2874,
+  0x000F9E971E014598, 0x000F89FA48A41DFC, 0x000F66C5F7F0302C, 0x000F1A5A4B331C4A]
+
+/-- For each strip of the normal, the scale taking a 52-bit draw to an abscissa. -/
+def wi : Array Float := #[
+  8.68362706080130616677e-16, 4.77933017572773682428e-17, 6.35435241740526230246e-17,
+  7.45487048124769627714e-17, 8.32936681579309972857e-17, 9.06806040505948228243e-17,
+  9.71486007656776183958e-17, 1.02947503142410192108e-16, 1.08234302884476839838e-16,
+  1.13114701961090307945e-16, 1.17663594570229211411e-16, 1.21936172787143633280e-16,
+  1.25974399146370927864e-16, 1.29810998862640315416e-16, 1.33472037368241227547e-16,
+  1.36978648425712032797e-16, 1.40348230012423820659e-16, 1.43595294520569430270e-16,
+  1.46732087423644219083e-16, 1.49769046683910367425e-16, 1.52715150035961979750e-16,
+  1.55578181694607639484e-16, 1.58364940092908853989e-16, 1.61081401752749279325e-16,
+  1.63732852039698532012e-16, 1.66323990584208352778e-16, 1.68859017086765964015e-16,
+  1.71341701765596607184e-16, 1.73775443658648593310e-16, 1.76163319230009959832e-16,
+  1.78508123169767272927e-16, 1.80812402857991522674e-16, 1.83078487648267501776e-16,
+  1.85308513886180189386e-16, 1.87504446393738816849e-16, 1.89668097007747596212e-16,
+  1.91801140648386198029e-16, 1.93905129306251037069e-16, 1.95981504266288244037e-16,
+  1.98031606831281739736e-16, 2.00056687762733300198e-16, 2.02057915620716538808e-16,
+  2.04036384154802118313e-16, 2.05993118874037063144e-16, 2.07929082904140197311e-16,
+  2.09845182223703516690e-16, 2.11742270357603418769e-16, 2.13621152594498681022e-16,
+  2.15482589785814580926e-16, 2.17327301775643674990e-16, 2.19155970504272708519e-16,
+  2.20969242822353175995e-16, 2.22767733047895534948e-16, 2.24552025294143552381e-16,
+  2.26322675592856786566e-16, 2.28080213834501706782e-16, 2.29825145544246839061e-16,
+  2.31557953510408037008e-16, 2.33279099280043561128e-16, 2.34989024534709550938e-16,
+  2.36688152357916037468e-16, 2.38376888404542434981e-16, 2.40055621981350627349e-16,
+  2.41724727046750252175e-16, 2.43384563137110286400e-16, 2.45035476226149539878e-16,
+  2.46677799523270498158e-16, 2.48311854216108767769e-16, 2.49937950162045242375e-16,
+  2.51556386532965786439e-16, 2.53167452417135826983e-16, 2.54771427381694417303e-16,
+  2.56368581998939683749e-16, 2.57959178339286723500e-16, 2.59543470433517070146e-16,
+  2.61121704706701939097e-16, 2.62694120385972564623e-16, 2.64260949884118951286e-16,
+  2.65822419160830680292e-16, 2.67378748063236329361e-16, 2.68930150647261591777e-16,
+  2.70476835481199518794e-16, 2.72019005932773206655e-16, 2.73556860440867908686e-16,
+  2.75090592773016664571e-16, 2.76620392269639032183e-16, 2.78146444075954410103e-16,
+  2.79668929362423005309e-16, 2.81188025534502074329e-16, 2.82703906432447923059e-16,
+  2.84216742521840606520e-16, 2.85726701075460149289e-16, 2.87233946347097994381e-16,
+  2.88738639737848191815e-16, 2.90240939955384233230e-16, 2.91741003166694553259e-16,
+  2.93238983144718163965e-16, 2.94735031409293489611e-16, 2.96229297362806647792e-16,
+  2.97721928420902891115e-16, 2.99213070138601307081e-16, 3.00702866332133102993e-16,
+  3.02191459196806151971e-16, 3.03678989421180184427e-16, 3.05165596297821922381e-16,
+  3.06651417830895451744e-16, 3.08136590840829717032e-16, 3.09621251066292253306e-16,
+  3.11105533263689296831e-16, 3.12589571304399892784e-16, 3.14073498269944617203e-16,
+  3.15557446545280064031e-16, 3.17041547910402852545e-16, 3.18525933630440648871e-16,
+  3.20010734544401137886e-16, 3.21496081152744704901e-16, 3.22982103703941557538e-16,
+  3.24468932280169778077e-16, 3.25956696882307838340e-16, 3.27445527514370671802e-16,
+  3.28935554267536967851e-16, 3.30426907403912838589e-16, 3.31919717440175233652e-16,
+  3.33414115231237245918e-16, 3.34910232054077845412e-16, 3.36408199691876507948e-16,
+  3.37908150518594979994e-16, 3.39410217584148914282e-16, 3.40914534700312603713e-16,
+  3.42421236527501816058e-16, 3.43930458662583133920e-16, 3.45442337727858401604e-16,
+  3.46957011461378353333e-16, 3.48474618808741370700e-16, 3.49995300016538099813e-16,
+  3.51519196727607440975e-16, 3.53046452078274009054e-16, 3.54577210797743572160e-16,
+  3.56111619309838843415e-16, 3.57649825837265051035e-16, 3.59191980508602994994e-16,
+  3.60738235468235137839e-16, 3.62288744989419151904e-16, 3.63843665590734438546e-16,
+  3.65403156156136995766e-16, 3.66967378058870090021e-16, 3.68536495289491401456e-16,
+  3.70110674588289834952e-16, 3.71690085582382297792e-16, 3.73274900927794352614e-16,
+  3.74865296456848868882e-16, 3.76461451331202869131e-16, 3.78063548200896037651e-16,
+  3.79671773369794425924e-16, 3.81286316967837738238e-16, 3.82907373130524317507e-16,
+  3.84535140186095955858e-16, 3.86169820850914927119e-16, 3.87811622433558721164e-16,
+  3.89460757048192620674e-16, 3.91117441837820542060e-16, 3.92781899208054153270e-16,
+  3.94454357072087711446e-16, 3.96135049107613542983e-16, 3.97824215026468259474e-16,
+  3.99522100857856502444e-16, 4.01228959246062907451e-16, 4.02945049763632792393e-16,
+  4.04670639241074995115e-16, 4.06406002114225038723e-16, 4.08151420790493873480e-16,
+  4.09907186035326643447e-16, 4.11673597380302570170e-16, 4.13450963554423599878e-16,
+  4.15239602940268833891e-16, 4.17039844056831587498e-16, 4.18852026071011229572e-16,
+  4.20676499339901510978e-16, 4.22513625986204937320e-16, 4.24363780509307796137e-16,
+  4.26227350434779809917e-16, 4.28104737005311666397e-16, 4.29996355916383230161e-16,
+  4.31902638100262944617e-16, 4.33824030562279080411e-16, 4.35760997273684900553e-16,
+  4.37714020125858747008e-16, 4.39683599951052137423e-16, 4.41670257615420348435e-16,
+  4.43674535190656726604e-16, 4.45696997211204306674e-16, 4.47738232024753387312e-16,
+  4.49798853244554968009e-16, 4.51879501313005876278e-16, 4.53980845187003400947e-16,
+  4.56103584156742206384e-16, 4.58248449810956667052e-16, 4.60416208163115281428e-16,
+  4.62607661954784567754e-16, 4.64823653154320737780e-16, 4.67065065671263059081e-16,
+  4.69332828309332890697e-16, 4.71627917983835129766e-16, 4.73951363232586715165e-16,
+  4.76304248053313737663e-16, 4.78687716104872284247e-16, 4.81102975314741720538e-16,
+  4.83551302941152515162e-16, 4.86034051145081195402e-16, 4.88552653135360343280e-16,
+  4.91108629959526955862e-16, 4.93703598024033454728e-16, 4.96339277440398725619e-16,
+  4.99017501309182245754e-16, 5.01740226071808946011e-16, 5.04509543081872748637e-16,
+  5.07327691573354207058e-16, 5.10197073234156184149e-16, 5.13120268630678373200e-16,
+  5.16100055774322824569e-16, 5.19139431175769859873e-16, 5.22241633800023428760e-16,
+  5.25410172417759732697e-16, 5.28648856950494511482e-16, 5.31961834533840037535e-16,
+  5.35353631181649688145e-16, 5.38829200133405320160e-16, 5.42393978220171234073e-16,
+  5.46053951907478041166e-16, 5.49815735089281410703e-16, 5.53686661246787600374e-16,
+  5.57674893292657647836e-16, 5.61789555355541665830e-16, 5.66040892008242216739e-16,
+  5.70440462129138908417e-16, 5.75001376891989523684e-16, 5.79738594572459365014e-16,
+  5.84669289345547900201e-16, 5.89813317647789942685e-16, 5.95193814964144415532e-16,
+  6.00837969627190832234e-16, 6.06778040933344851394e-16, 6.13052720872528159123e-16,
+  6.19708989458162555387e-16, 6.26804696330128439415e-16, 6.34412240712750598627e-16,
+  6.42623965954805540945e-16, 6.51560331734499356881e-16, 6.61382788509766415145e-16,
+  6.72315046250558662913e-16, 6.84680341756425875856e-16, 6.98971833638761995415e-16,
+  7.15999493483066421560e-16, 7.37242430179879890722e-16, 7.65893637080557275482e-16,
+  8.11384933765648418565e-16]
+
+/-- For each strip of the normal, the value of the density at its edge. -/
+def fi : Array Float := #[
+  1.00000000000000000000e+00, 9.77101701267671596263e-01, 9.59879091800106665211e-01,
+  9.45198953442299649730e-01, 9.32060075959230460718e-01, 9.19991505039347012840e-01,
+  9.08726440052130879366e-01, 8.98095921898343418910e-01, 8.87984660755833377088e-01,
+  8.78309655808917399966e-01, 8.69008688036857046555e-01, 8.60033621196331532488e-01,
+  8.51346258458677951353e-01, 8.42915653112204177333e-01, 8.34716292986883434679e-01,
+  8.26726833946221373317e-01, 8.18929191603702366642e-01, 8.11307874312656274185e-01,
+  8.03849483170964274059e-01, 7.96542330422958966274e-01, 7.89376143566024590648e-01,
+  7.82341832654802504798e-01, 7.75431304981187174974e-01, 7.68637315798486264740e-01,
+  7.61953346836795386565e-01, 7.55373506507096115214e-01, 7.48892447219156820459e-01,
+  7.42505296340151055290e-01, 7.36207598126862650112e-01, 7.29995264561476231435e-01,
+  7.23864533468630222401e-01, 7.17811932630721960535e-01, 7.11834248878248421200e-01,
+  7.05928501332754310127e-01, 7.00091918136511615067e-01, 6.94321916126116711609e-01,
+  6.88616083004671808432e-01, 6.82972161644994857355e-01, 6.77388036218773526009e-01,
+  6.71861719897082099173e-01, 6.66391343908750100056e-01, 6.60975147776663107813e-01,
+  6.55611470579697264149e-01, 6.50298743110816701574e-01, 6.45035480820822293424e-01,
+  6.39820277453056585060e-01, 6.34651799287623608059e-01, 6.29528779924836690007e-01,
+  6.24450015547026504592e-01, 6.19414360605834324325e-01, 6.14420723888913888899e-01,
+  6.09468064925773433949e-01, 6.04555390697467776029e-01, 5.99681752619125263415e-01,
+  5.94846243767987448159e-01, 5.90047996332826008015e-01, 5.85286179263371453274e-01,
+  5.80559996100790898232e-01, 5.75868682972353718164e-01, 5.71211506735253227163e-01,
+  5.66587763256164445025e-01, 5.61996775814524340831e-01, 5.57437893618765945014e-01,
+  5.52910490425832290562e-01, 5.48413963255265812791e-01, 5.43947731190026262382e-01,
+  5.39511234256952132426e-01, 5.35103932380457614215e-01, 5.30725304403662057062e-01,
+  5.26374847171684479008e-01, 5.22052074672321841931e-01, 5.17756517229756352272e-01,
+  5.13487720747326958914e-01, 5.09245245995747941592e-01, 5.05028667943468123624e-01,
+  5.00837575126148681903e-01, 4.96671569052489714213e-01, 4.92530263643868537748e-01,
+  4.88413284705458028423e-01, 4.84320269426683325253e-01, 4.80250865909046753544e-01,
+  4.76204732719505863248e-01, 4.72181538467730199660e-01, 4.68180961405693596422e-01,
+  4.64202689048174355069e-01, 4.60246417812842867345e-01, 4.56311852678716434184e-01,
+  4.52398706861848520777e-01, 4.48506701507203064949e-01, 4.44635565395739396077e-01,
+  4.40785034665803987508e-01, 4.36954852547985550526e-01, 4.33144769112652261445e-01,
+  4.29354541029441427735e-01, 4.25583931338021970170e-01, 4.21832709229495894654e-01,
+  4.18100649837848226120e-01, 4.14387534040891125642e-01, 4.10693148270188157500e-01,
+  4.07017284329473372217e-01, 4.03359739221114510510e-01, 3.99720314980197222177e-01,
+  3.96098818515832451492e-01, 3.92495061459315619512e-01, 3.88908860018788715696e-01,
+  3.85340034840077283462e-01, 3.81788410873393657674e-01, 3.78253817245619183840e-01,
+  3.74736087137891138443e-01, 3.71235057668239498696e-01, 3.67750569779032587814e-01,
+  3.64282468129004055601e-01, 3.60830600989648031529e-01, 3.57394820145780500731e-01,
+  3.53974980800076777232e-01, 3.50570941481406106455e-01, 3.47182563956793643900e-01,
+  3.43809713146850715049e-01, 3.40452257044521866547e-01, 3.37110066637006045021e-01,
+  3.33783015830718454708e-01, 3.30470981379163586400e-01, 3.27173842813601400970e-01,
+  3.23891482376391093290e-01, 3.20623784956905355514e-01, 3.17370638029913609834e-01,
+  3.14131931596337177215e-01, 3.10907558126286509559e-01, 3.07697412504292056035e-01,
+  3.04501391976649993243e-01, 3.01319396100803049698e-01, 2.98151326696685481377e-01,
+  2.94997087799961810184e-01, 2.91856585617095209972e-01, 2.88729728482182923521e-01,
+  2.85616426815501756042e-01, 2.82516593083707578948e-01, 2.79430141761637940157e-01,
+  2.76356989295668320494e-01, 2.73297054068577072172e-01, 2.70250256365875463072e-01,
+  2.67216518343561471038e-01, 2.64195763997261190426e-01, 2.61187919132721213522e-01,
+  2.58192911337619235290e-01, 2.55210669954661961700e-01, 2.52241126055942177508e-01,
+  2.49284212418528522415e-01, 2.46339863501263828249e-01, 2.43408015422750312329e-01,
+  2.40488605940500588254e-01, 2.37581574431238090606e-01, 2.34686861872330010392e-01,
+  2.31804410824338724684e-01, 2.28934165414680340644e-01, 2.26076071322380278694e-01,
+  2.23230075763917484855e-01, 2.20396127480151998723e-01, 2.17574176724331130872e-01,
+  2.14764175251173583536e-01, 2.11966076307030182324e-01, 2.09179834621125076977e-01,
+  2.06405406397880797353e-01, 2.03642749310334908452e-01, 2.00891822494656591136e-01,
+  1.98152586545775138971e-01, 1.95425003514134304483e-01, 1.92709036903589175926e-01,
+  1.90004651670464985713e-01, 1.87311814223800304768e-01, 1.84630492426799269756e-01,
+  1.81960655599522513892e-01, 1.79302274522847582272e-01, 1.76655321443734858455e-01,
+  1.74019770081838553999e-01, 1.71395595637505754327e-01, 1.68782774801211288285e-01,
+  1.66181285764481906364e-01, 1.63591108232365584074e-01, 1.61012223437511009516e-01,
+  1.58444614155924284882e-01, 1.55888264724479197465e-01, 1.53343161060262855866e-01,
+  1.50809290681845675763e-01, 1.48286642732574552861e-01, 1.45775208005994028060e-01,
+  1.43274978973513461566e-01, 1.40785949814444699690e-01, 1.38308116448550733057e-01,
+  1.35841476571253755301e-01, 1.33386029691669155683e-01, 1.30941777173644358090e-01,
+  1.28508722279999570981e-01, 1.26086870220185887081e-01, 1.23676228201596571932e-01,
+  1.21276805484790306533e-01, 1.18888613442910059947e-01, 1.16511665625610869035e-01,
+  1.14145977827838487895e-01, 1.11791568163838089811e-01, 1.09448457146811797824e-01,
+  1.07116667774683801961e-01, 1.04796225622487068629e-01, 1.02487158941935246892e-01,
+  1.00189498768810017482e-01, 9.79032790388624646338e-02, 9.56285367130089991594e-02,
+  9.33653119126910124859e-02, 9.11136480663737591268e-02, 8.88735920682758862021e-02,
+  8.66451944505580717859e-02, 8.44285095703534715916e-02, 8.22235958132029043366e-02,
+  8.00305158146630696292e-02, 7.78493367020961224423e-02, 7.56801303589271778804e-02,
+  7.35229737139813238622e-02, 7.13779490588904025339e-02, 6.92451443970067553879e-02,
+  6.71246538277884968737e-02, 6.50165779712428976156e-02, 6.29210244377581412456e-02,
+  6.08381083495398780614e-02, 5.87679529209337372930e-02, 5.67106901062029017391e-02,
+  5.46664613248889208474e-02, 5.26354182767921896513e-02, 5.06177238609477817000e-02,
+  4.86135532158685421122e-02, 4.66230949019303814174e-02, 4.46465522512944634759e-02,
+  4.26841449164744590750e-02, 4.07361106559409394401e-02, 3.88027074045261474722e-02,
+  3.68842156885673053135e-02, 3.49809414617161251737e-02, 3.30932194585785779961e-02,
+  3.12214171919203004046e-02, 2.93659397581333588001e-02, 2.75272356696031131329e-02,
+  2.57058040085489103443e-02, 2.39022033057958785407e-02, 2.21170627073088502113e-02,
+  2.03510962300445102935e-02, 1.86051212757246224594e-02, 1.68800831525431419000e-02,
+  1.51770883079353092332e-02, 1.34974506017398673818e-02, 1.18427578579078790488e-02,
+  1.02149714397014590439e-02, 8.61658276939872638800e-03, 7.05087547137322242369e-03,
+  5.52240329925099155545e-03, 4.03797259336302356153e-03, 2.60907274610215926189e-03,
+  1.26028593049859797236e-03]
+
+/-- For each strip of the exponential, the largest 53-bit draw landing in its rectangular part. -/
+def ke : Array UInt64 := #[
+  0x001C5214272497C6, 0x0000000000000000, 0x00137D5BD79C317E, 0x00186EF58E3F3C10,
+  0x001A9BB7320EB0AE, 0x001BD127F719447C, 0x001C951D0F88651A, 0x001D1BFE2D5C3972,
+  0x001D7E5BD56B18B2, 0x001DC934DD172C70, 0x001E0409DFAC9DC8, 0x001E337B71D47836,
+  0x001E5A8B177CB7A2, 0x001E7B42096F046C, 0x001E970DAF08AE3E, 0x001EAEF5B14EF09E,
+  0x001EC3BD07B46556, 0x001ED5F6F08799CE, 0x001EE614AE6E5688, 0x001EF46ECA361CD0,
+  0x001F014B76DDD4A4, 0x001F0CE313A796B6, 0x001F176369F1F77A, 0x001F20F20C452570,
+  0x001F29AE1951A874, 0x001F31B18FB95532, 0x001F39125157C106, 0x001F3FE2EB6E694C,
+  0x001F463332D788FA, 0x001F4C10BF1D3A0E, 0x001F51874C5C3322, 0x001F56A109C3ECC0,
+  0x001F5B66D9099996, 0x001F5FE08210D08C, 0x001F6414DD445772, 0x001F6809F6859678,
+  0x001F6BC52A2B02E6, 0x001F6F4B3D32E4F4, 0x001F72A07190F13A, 0x001F75C8974D09D6,
+  0x001F78C71B045CC0, 0x001F7B9F12413FF4, 0x001F7E5346079F8A, 0x001F80E63BE21138,
+  0x001F835A3DAD9162, 0x001F85B16056B912, 0x001F87ED89B24262, 0x001F8A10759374FA,
+  0x001F8C1BBA3D39AC, 0x001F8E10CC45D04A, 0x001F8FF102013E16, 0x001F91BD968358E0,
+  0x001F9377AC47AFD8, 0x001F95204F8B64DA, 0x001F96B878633892, 0x001F98410C968892,
+  0x001F99BAE146BA80, 0x001F9B26BC697F00, 0x001F9C85561B717A, 0x001F9DD759CFD802,
+  0x001F9F1D6761A1CE, 0x001FA058140936C0, 0x001FA187EB3A3338, 0x001FA2AD6F6BC4FC,
+  0x001FA3C91ACE0682, 0x001FA4DB5FEE6AA2, 0x001FA5E4AA4D097C, 0x001FA6E55EE46782,
+  0x001FA7DDDCA51EC4, 0x001FA8CE7CE6A874, 0x001FA9B793CE5FEE, 0x001FAA9970ADB858,
+  0x001FAB745E588232, 0x001FAC48A3740584, 0x001FAD1682BF9FE8, 0x001FADDE3B5782C0,
+  0x001FAEA008F21D6C, 0x001FAF5C2418B07E, 0x001FB012C25B7A12, 0x001FB0C41681DFF4,
+  0x001FB17050B6F1FA, 0x001FB2179EB2963A, 0x001FB2BA2BDFA84A, 0x001FB358217F4E18,
+  0x001FB3F1A6C9BE0C, 0x001FB486E10CACD6, 0x001FB517F3C793FC, 0x001FB5A500C5FDAA,
+  0x001FB62E2837FE58, 0x001FB6B388C9010A, 0x001FB7353FB50798, 0x001FB7B368DC7DA8,
+  0x001FB82E1ED6BA08, 0x001FB8A57B0347F6, 0x001FB919959A0F74, 0x001FB98A85BA7204,
+  0x001FB9F861796F26, 0x001FBA633DEEE286, 0x001FBACB2F41EC16, 0x001FBB3048B49144,
+  0x001FBB929CAEA4E2, 0x001FBBF23CC8029E, 0x001FBC4F39D22994, 0x001FBCA9A3E140D4,
+  0x001FBD018A548F9E, 0x001FBD56FBDE729C, 0x001FBDAA068BD66A, 0x001FBDFAB7CB3F40,
+  0x001FBE491C7364DE, 0x001FBE9540C9695E, 0x001FBEDF3086B128, 0x001FBF26F6DE6174,
+  0x001FBF6C9E828AE2, 0x001FBFB031A904C4, 0x001FBFF1BA0FFDB0, 0x001FC03141024588,
+  0x001FC06ECF5B54B2, 0x001FC0AA6D8B1426, 0x001FC0E42399698A, 0x001FC11BF9298A64,
+  0x001FC151F57D1942, 0x001FC1861F770F4A, 0x001FC1B87D9E74B4, 0x001FC1E91620EA42,
+  0x001FC217EED505DE, 0x001FC2450D3C83FE, 0x001FC27076864FC2, 0x001FC29A2F90630E,
+  0x001FC2C23CE98046, 0x001FC2E8A2D2C6B4, 0x001FC30D654122EC, 0x001FC33087DE9C0E,
+  0x001FC3520E0B7EC6, 0x001FC371FADF66F8, 0x001FC390512A2886, 0x001FC3AD137497FA,
+  0x001FC3C844013348, 0x001FC3E1E4CCAB40, 0x001FC3F9F78E4DA8, 0x001FC4107DB85060,
+  0x001FC4257877FD68, 0x001FC438E8B5BFC6, 0x001FC44ACF15112A, 0x001FC45B2BF447E8,
+  0x001FC469FF6C4504, 0x001FC477495001B2, 0x001FC483092BFBB8, 0x001FC48D3E457FF6,
+  0x001FC495E799D21A, 0x001FC49D03DD30B0, 0x001FC4A29179B432, 0x001FC4A68E8E07FC,
+  0x001FC4A8F8EBFB8C, 0x001FC4A9CE16EA9E, 0x001FC4A90B41FA34, 0x001FC4A6AD4E28A0,
+  0x001FC4A2B0C82E74, 0x001FC49D11E62DE2, 0x001FC495CC852DF4, 0x001FC48CDC265EC0,
+  0x001FC4823BEC237A, 0x001FC475E696DEE6, 0x001FC467D6817E82, 0x001FC458059DC036,
+  0x001FC4466D702E20, 0x001FC433070BCB98, 0x001FC41DCB0D6E0E, 0x001FC406B196BBF6,
+  0x001FC3EDB248CB62, 0x001FC3D2C43E593C, 0x001FC3B5DE0591B4, 0x001FC396F599614C,
+  0x001FC376005A4592, 0x001FC352F3069370, 0x001FC32DC1B22818, 0x001FC3065FBD7888,
+  0x001FC2DCBFCBF262, 0x001FC2B0D3B99F9E, 0x001FC2828C8FFCF0, 0x001FC251DA79F164,
+  0x001FC21EACB6D39E, 0x001FC1E8F18C6756, 0x001FC1B09637BB3C, 0x001FC17586DCCD10,
+  0x001FC137AE74D6B6, 0x001FC0F6F6BB2414, 0x001FC0B348184DA4, 0x001FC06C898BAFF0,
+  0x001FC022A092F364, 0x001FBFD5710F72B8, 0x001FBF84DD29488E, 0x001FBF30C52FC60A,
+  0x001FBED907770CC6, 0x001FBE7D80327DDA, 0x001FBE1E094BA614, 0x001FBDBA7A354408,
+  0x001FBD52A7B9F826, 0x001FBCE663C6201A, 0x001FBC757D2C4DE4, 0x001FBBFFBF63B7AA,
+  0x001FBB84F23FE6A2, 0x001FBB04D9A0D18C, 0x001FBA7F351A70AC, 0x001FB9F3BF92B618,
+  0x001FB9622ED4ABFC, 0x001FB8CA33174A16, 0x001FB82B76765B54, 0x001FB7859C5B895C,
+  0x001FB6D840D55594, 0x001FB622F7D96942, 0x001FB5654C6F37E0, 0x001FB49EBFBF69D2,
+  0x001FB3CEC803E746, 0x001FB2F4CF539C3E, 0x001FB21032442852, 0x001FB1203E5A9604,
+  0x001FB0243042E1C2, 0x001FAF1B31C479A6, 0x001FAE045767E104, 0x001FACDE9DBF2D72,
+  0x001FABA8E640060A, 0x001FAA61F399FF28, 0x001FA908656F66A2, 0x001FA79AB3508D3C,
+  0x001FA61726D1F214, 0x001FA47BD48BEA00, 0x001FA2C693C5C094, 0x001FA0F4F47DF314,
+  0x001F9F04336BBE0A, 0x001F9CF12B79F9BC, 0x001F9AB84415ABC4, 0x001F98555B782FB8,
+  0x001F95C3ABD03F78, 0x001F92FDA9CEF1F2, 0x001F8FFCDA9AE41C, 0x001F8CB99E7385F8,
+  0x001F892AEC479606, 0x001F8545F904DB8E, 0x001F80FDC336039A, 0x001F7C427839E926,
+  0x001F7700A3582ACC, 0x001F71200F1A241C, 0x001F6A8234B7352A, 0x001F630000A8E266,
+  0x001F5A66904FE3C4, 0x001F50724ECE1172, 0x001F44C7665C6FDA, 0x001F36E5A38A59A2,
+  0x001F26143450340A, 0x001F113E047B0414, 0x001EF6AEFA57CBE6, 0x001ED38CA188151E,
+  0x001EA2A61E122DB0, 0x001E5961C78B267C, 0x001DDDF62BAC0BB0, 0x001CDB4DD9E4E8C0]
+
+/-- For each strip of the exponential, the scale taking a 53-bit draw to an abscissa. -/
+def we : Array Float := #[
+  9.655740063209182975e-16, 7.089014243955414331e-18, 1.163941249669122378e-17,
+  1.524391512353216015e-17, 1.833284885723743916e-17, 2.108965109464486630e-17,
+  2.361128077843138196e-17, 2.595595772310893952e-17, 2.816173554197752338e-17,
+  3.025504130321382330e-17, 3.225508254836375280e-17, 3.417632340185027033e-17,
+  3.602996978734452488e-17, 3.782490776869649048e-17, 3.956832198097553231e-17,
+  4.126611778175946428e-17, 4.292321808442525631e-17, 4.454377743282371417e-17,
+  4.613133981483185932e-17, 4.768895725264635940e-17, 4.921928043727962847e-17,
+  5.072462904503147014e-17, 5.220704702792671737e-17, 5.366834661718192181e-17,
+  5.511014372835094717e-17, 5.653388673239667134e-17, 5.794088004852766616e-17,
+  5.933230365208943081e-17, 6.070922932847179572e-17, 6.207263431163193485e-17,
+  6.342341280303076511e-17, 6.476238575956142121e-17, 6.609030925769405241e-17,
+  6.740788167872722244e-17, 6.871574991183812442e-17, 7.001451473403929616e-17,
+  7.130473549660643409e-17, 7.258693422414648352e-17, 7.386159921381791997e-17,
+  7.512918820723728089e-17, 7.639013119550825792e-17, 7.764483290797848102e-17,
+  7.889367502729790548e-17, 8.013701816675454434e-17, 8.137520364041762206e-17,
+  8.260855505210038174e-17, 8.383737972539139383e-17, 8.506196999385323132e-17,
+  8.628260436784112996e-17, 8.749954859216182511e-17, 8.871305660690252281e-17,
+  8.992337142215357066e-17, 9.113072591597909173e-17, 9.233534356381788123e-17,
+  9.353743910649128938e-17, 9.473721916312949566e-17, 9.593488279457997317e-17,
+  9.713062202221521206e-17, 9.832462230649511362e-17, 9.951706298915071878e-17,
+  1.007081177024294931e-16, 1.018979547484694078e-16, 1.030867374515421954e-16,
+  1.042746244856188556e-16, 1.054617701794576406e-16, 1.066483248011914702e-16,
+  1.078344348241948498e-16, 1.090202431758350473e-16, 1.102058894705578110e-16,
+  1.113915102286197502e-16, 1.125772390816567488e-16, 1.137632069661684705e-16,
+  1.149495423059009298e-16, 1.161363711840218308e-16, 1.173238175059045788e-16,
+  1.185120031532669434e-16, 1.197010481303465158e-16, 1.208910707027385520e-16,
+  1.220821875294706151e-16, 1.232745137888415193e-16, 1.244681632985112523e-16,
+  1.256632486302898513e-16, 1.268598812200397542e-16, 1.280581714730749379e-16,
+  1.292582288654119552e-16, 1.304601620412028847e-16, 1.316640789066572582e-16,
+  1.328700867207380889e-16, 1.340782921828999433e-16, 1.352888015181175458e-16,
+  1.365017205594397770e-16, 1.377171548282880964e-16, 1.389352096127063919e-16,
+  1.401559900437571538e-16, 1.413796011702485188e-16, 1.426061480319665444e-16,
+  1.438357357315790180e-16, 1.450684695053687684e-16, 1.463044547929475721e-16,
+  1.475437973060951633e-16, 1.487866030968626066e-16, 1.500329786250736949e-16,
+  1.512830308253539427e-16, 1.525368671738125550e-16, 1.537945957544996933e-16,
+  1.550563253257577148e-16, 1.563221653865837505e-16, 1.575922262431176140e-16,
+  1.588666190753684151e-16, 1.601454560042916733e-16, 1.614288501593278662e-16,
+  1.627169157465130500e-16, 1.640097681172717950e-16, 1.653075238380036909e-16,
+  1.666103007605742067e-16, 1.679182180938228863e-16, 1.692313964762022267e-16,
+  1.705499580496629830e-16, 1.718740265349031656e-16, 1.732037273081008369e-16,
+  1.745391874792533975e-16, 1.758805359722491379e-16, 1.772279036068006489e-16,
+  1.785814231823732619e-16, 1.799412295642463721e-16, 1.813074597718501559e-16,
+  1.826802530695252266e-16, 1.840597510598587828e-16, 1.854460977797569461e-16,
+  1.868394397994192684e-16, 1.882399263243892051e-16, 1.896477093008616722e-16,
+  1.910629435244376536e-16, 1.924857867525243818e-16, 1.939163998205899420e-16,
+  1.953549467624909132e-16, 1.968015949351037382e-16, 1.982565151475019047e-16,
+  1.997198817949342081e-16, 2.011918729978734671e-16, 2.026726707464198289e-16,
+  2.041624610503588774e-16, 2.056614340951917875e-16, 2.071697844044737034e-16,
+  2.086877110088159721e-16, 2.102154176219292789e-16, 2.117531128241075913e-16,
+  2.133010102535779087e-16, 2.148593288061663316e-16, 2.164282928437604723e-16,
+  2.180081324120784027e-16, 2.195990834682870728e-16, 2.212013881190495942e-16,
+  2.228152948696180545e-16, 2.244410588846308588e-16, 2.260789422613173739e-16,
+  2.277292143158621037e-16, 2.293921518837311354e-16, 2.310680396348213318e-16,
+  2.327571704043534613e-16, 2.344598455404957859e-16, 2.361763752697773994e-16,
+  2.379070790814276700e-16, 2.396522861318623520e-16, 2.414123356706293277e-16,
+  2.431875774892255956e-16, 2.449783723943070217e-16, 2.467850927069288738e-16,
+  2.486081227895851719e-16, 2.504478596029557040e-16, 2.523047132944217013e-16,
+  2.541791078205812227e-16, 2.560714816061770759e-16, 2.579822882420530896e-16,
+  2.599119972249746917e-16, 2.618610947423924219e-16, 2.638300845054942823e-16,
+  2.658194886341845120e-16, 2.678298485979525166e-16, 2.698617262169488933e-16,
+  2.719157047279818500e-16, 2.739923899205814823e-16, 2.760924113487617126e-16,
+  2.782164236246436081e-16, 2.803651078006983464e-16, 2.825391728480253184e-16,
+  2.847393572388174091e-16, 2.869664306419817679e-16, 2.892211957417995598e-16,
+  2.915044901905293183e-16, 2.938171887070028633e-16, 2.961602053345465687e-16,
+  2.985344958730045276e-16, 3.009410605012618141e-16, 3.033809466085003416e-16,
+  3.058552518544860874e-16, 3.083651274815310004e-16, 3.109117819034266344e-16,
+  3.134964845996663118e-16, 3.161205703467105734e-16, 3.187854438219713117e-16,
+  3.214925846206797361e-16, 3.242435527309451638e-16, 3.270399945182240440e-16,
+  3.298836492772283149e-16, 3.327763564171671408e-16, 3.357200633553244075e-16,
+  3.387168342045505162e-16, 3.417688593525636996e-16, 3.448784660453423890e-16,
+  3.480481301037442286e-16, 3.512804889222979418e-16, 3.545783559224791863e-16,
+  3.579447366604276541e-16, 3.613828468219060593e-16, 3.648961323764542545e-16,
+  3.684882922095621322e-16, 3.721633036080207290e-16, 3.759254510416256532e-16,
+  3.797793587668874387e-16, 3.837300278789213687e-16, 3.877828785607895292e-16,
+  3.919437984311428867e-16, 3.962191980786774996e-16, 4.006160751056541688e-16,
+  4.051420882956573177e-16, 4.098056438903062509e-16, 4.146159964290904582e-16,
+  4.195833672073398926e-16, 4.247190841824385048e-16, 4.300357481667470702e-16,
+  4.355474314693952008e-16, 4.412699169036069903e-16, 4.472209874259932285e-16,
+  4.534207798565834480e-16, 4.598922204905932469e-16, 4.666615664711475780e-16,
+  4.737590853262492027e-16, 4.812199172829237933e-16, 4.890851827392209900e-16,
+  4.974034236191939753e-16, 5.062325072144159699e-16, 5.156421828878082953e-16,
+  5.257175802022274839e-16, 5.365640977112021618e-16, 5.483144034258703912e-16,
+  5.611387454675159622e-16, 5.752606481503331688e-16, 5.909817641652102998e-16,
+  6.087231416180907671e-16, 6.290979034877557049e-16, 6.530492053564040799e-16,
+  6.821393079028928626e-16, 7.192444966089361564e-16, 7.706095350032096755e-16,
+  8.545517038584027421e-16]
+
+/-- For each strip of the exponential, the value of the density at its edge. -/
+def fe : Array Float := #[
+  1.000000000000000000e+00, 9.381436808621747003e-01, 9.004699299257464817e-01,
+  8.717043323812035949e-01, 8.477855006239896074e-01, 8.269932966430503241e-01,
+  8.084216515230083777e-01, 7.915276369724956185e-01, 7.759568520401155522e-01,
+  7.614633888498962833e-01, 7.478686219851951034e-01, 7.350380924314234843e-01,
+  7.228676595935720206e-01, 7.112747608050760117e-01, 7.001926550827881623e-01,
+  6.895664961170779872e-01, 6.793505722647653622e-01, 6.695063167319247333e-01,
+  6.600008410789997004e-01, 6.508058334145710999e-01, 6.418967164272660897e-01,
+  6.332519942143660652e-01, 6.248527387036659775e-01, 6.166821809152076561e-01,
+  6.087253820796220127e-01, 6.009689663652322267e-01, 5.934009016917334289e-01,
+  5.860103184772680329e-01, 5.787873586028450257e-01, 5.717230486648258170e-01,
+  5.648091929124001709e-01, 5.580382822625874484e-01, 5.514034165406412891e-01,
+  5.448982376724396115e-01, 5.385168720028619127e-01, 5.322538802630433219e-01,
+  5.261042139836197284e-01, 5.200631773682335979e-01, 5.141263938147485613e-01,
+  5.082897764106428795e-01, 5.025495018413477233e-01, 4.969019872415495476e-01,
+  4.913438695940325340e-01, 4.858719873418849144e-01, 4.804833639304542103e-01,
+  4.751751930373773747e-01, 4.699448252839599771e-01, 4.647897562504261781e-01,
+  4.597076156421376902e-01, 4.546961574746155033e-01, 4.497532511627549967e-01,
+  4.448768734145485126e-01, 4.400651008423538957e-01, 4.353161032156365740e-01,
+  4.306281372884588343e-01, 4.259995411430343437e-01, 4.214287289976165751e-01,
+  4.169141864330028757e-01, 4.124544659971611793e-01, 4.080481831520323954e-01,
+  4.036940125305302773e-01, 3.993906844752310725e-01, 3.951369818332901573e-01,
+  3.909317369847971069e-01, 3.867738290841376547e-01, 3.826621814960098344e-01,
+  3.785957594095807899e-01, 3.745735676159021588e-01, 3.705946484351460013e-01,
+  3.666580797815141568e-01, 3.627629733548177748e-01, 3.589084729487497794e-01,
+  3.550937528667874599e-01, 3.513180164374833381e-01, 3.475804946216369817e-01,
+  3.438804447045024082e-01, 3.402171490667800224e-01, 3.365899140286776059e-01,
+  3.329980687618089852e-01, 3.294409642641363267e-01, 3.259179723935561879e-01,
+  3.224284849560891675e-01, 3.189719128449572394e-01, 3.155476852271289490e-01,
+  3.121552487741795501e-01, 3.087940669345601852e-01, 3.054636192445902565e-01,
+  3.021634006756935276e-01, 2.988929210155817917e-01, 2.956517042812611962e-01,
+  2.924392881618925744e-01, 2.892552234896777485e-01, 2.860990737370768255e-01,
+  2.829704145387807457e-01, 2.798688332369729248e-01, 2.767939284485173568e-01,
+  2.737453096528029706e-01, 2.707225967990600224e-01, 2.677254199320447947e-01,
+  2.647534188350622042e-01, 2.618062426893629779e-01, 2.588835497490162285e-01,
+  2.559850070304153791e-01, 2.531102900156294577e-01, 2.502590823688622956e-01,
+  2.474310756653276266e-01, 2.446259691318921070e-01, 2.418434693988772144e-01,
+  2.390832902624491774e-01, 2.363451524570596429e-01, 2.336287834374333461e-01,
+  2.309339171696274118e-01, 2.282602939307167011e-01, 2.256076601166840667e-01,
+  2.229757680581201940e-01, 2.203643758433594946e-01, 2.177732471487005272e-01,
+  2.152021510753786837e-01, 2.126508619929782795e-01, 2.101191593889882581e-01,
+  2.076068277242220372e-01, 2.051136562938377095e-01, 2.026394390937090173e-01,
+  2.001839746919112650e-01, 1.977470661050988732e-01, 1.953285206795632167e-01,
+  1.929281499767713515e-01, 1.905457696631953912e-01, 1.881811994042543179e-01,
+  1.858342627621971110e-01, 1.835047870977674633e-01, 1.811926034754962889e-01,
+  1.788975465724783054e-01, 1.766194545904948843e-01, 1.743581691713534942e-01,
+  1.721135353153200598e-01, 1.698854013025276610e-01, 1.676736186172501919e-01,
+  1.654780418749360049e-01, 1.632985287519018169e-01, 1.611349399175920349e-01,
+  1.589871389693142123e-01, 1.568549923693652315e-01, 1.547383693844680830e-01,
+  1.526371420274428570e-01, 1.505511850010398944e-01, 1.484803756438667910e-01,
+  1.464245938783449441e-01, 1.443837221606347754e-01, 1.423576454324722018e-01,
+  1.403462510748624548e-01, 1.383494288635802039e-01, 1.363670709264288572e-01,
+  1.343990717022136294e-01, 1.324453279013875218e-01, 1.305057384683307731e-01,
+  1.285802045452281717e-01, 1.266686294375106714e-01, 1.247709185808309612e-01,
+  1.228869795095451356e-01, 1.210167218266748335e-01, 1.191600571753276827e-01,
+  1.173168992115555670e-01, 1.154871635786335338e-01, 1.136707678827443141e-01,
+  1.118676316700562973e-01, 1.100776764051853845e-01, 1.083008254510337970e-01,
+  1.065370040500016602e-01, 1.047861393065701724e-01, 1.030481601712577161e-01,
+  1.013229974259536315e-01, 9.961058367063713170e-02, 9.791085331149219917e-02,
+  9.622374255043279756e-02, 9.454918937605585882e-02, 9.288713355604354127e-02,
+  9.123751663104015530e-02, 8.960028191003285847e-02, 8.797537446727021759e-02,
+  8.636274114075691288e-02, 8.476233053236811865e-02, 8.317409300963238272e-02,
+  8.159798070923741931e-02, 8.003394754231990538e-02, 7.848194920160642130e-02,
+  7.694194317048050347e-02, 7.541388873405840965e-02, 7.389774699236474620e-02,
+  7.239348087570873780e-02, 7.090105516237182881e-02, 6.942043649872875477e-02,
+  6.795159342193660135e-02, 6.649449638533977414e-02, 6.504911778675374900e-02,
+  6.361543199980733421e-02, 6.219341540854099459e-02, 6.078304644547963265e-02,
+  5.938430563342026597e-02, 5.799717563120065922e-02, 5.662164128374287675e-02,
+  5.525768967669703741e-02, 5.390531019604608703e-02, 5.256449459307169225e-02,
+  5.123523705512628146e-02, 4.991753428270637172e-02, 4.861138557337949667e-02,
+  4.731679291318154762e-02, 4.603376107617516977e-02, 4.476229773294328196e-02,
+  4.350241356888818328e-02, 4.225412241331623353e-02, 4.101744138041481941e-02,
+  3.979239102337412542e-02, 3.857899550307485742e-02, 3.737728277295936097e-02,
+  3.618728478193142251e-02, 3.500903769739741045e-02, 3.384258215087432992e-02,
+  3.268796350895953468e-02, 3.154523217289360859e-02, 3.041444391046660423e-02,
+  2.929566022463739317e-02, 2.818894876397863569e-02, 2.709438378095579969e-02,
+  2.601204664513421735e-02, 2.494202641973178314e-02, 2.388442051155817078e-02,
+  2.283933540638524023e-02, 2.180688750428358066e-02, 2.078720407257811723e-02,
+  1.978042433800974303e-02, 1.878670074469603046e-02, 1.780620041091136169e-02,
+  1.683910682603994777e-02, 1.588562183997316302e-02, 1.494596801169114850e-02,
+  1.402039140318193759e-02, 1.310916493125499106e-02, 1.221259242625538123e-02,
+  1.133101359783459695e-02, 1.046481018102997894e-02, 9.614413642502209895e-03,
+  8.780314985808975251e-03, 7.963077438017040002e-03, 7.163353183634983863e-03,
+  6.381905937319179087e-03, 5.619642207205483020e-03, 4.877655983542392333e-03,
+  4.157295120833795314e-03, 3.460264777836904049e-03, 2.788798793574076128e-03,
+  2.145967743718906265e-03, 1.536299780301572356e-03, 9.672692823271745359e-04,
+  4.541343538414967652e-04]
+
+end NumLean.Ziggurat
diff --git a/RandomDo/NumLean/ZigguratSampler.lean b/RandomDo/NumLean/ZigguratSampler.lean
new file mode 100644
index 0000000..ba56ae3
--- /dev/null
+++ b/RandomDo/NumLean/ZigguratSampler.lean
@@ -0,0 +1,87 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import RandomDo.NumLean.PCG64
+public import FFI.Float
+
+/-!
+# The ziggurat sampler
+
+The sampler shape numpy's `random_standard_normal` and `random_standard_exponential` share, the
+ziggurat of Marsaglia and Tsang. The density is covered by 256 strips of equal area; a draw picks a
+strip and a point in it, and is accepted at once when the point falls in the strip's rectangular
+part, which is where about 99% of the draws end.
+
+The tables the two samplers read live in `RandomDo.NumLean.Ziggurat`, and the samplers themselves
+in `RandomDo.NumLean.Distributions`; this file holds only the shape they share, which is generic in
+the tables `k`, `w` and `f` and in the three functions that say how a 64-bit output is cut up, what
+the density is, and how the base strip's tail is sampled.
+
+## Main definitions
+
+* `ziggurat`: the sampler, whose fast path is inlined into its caller
+* `zigguratSlow`: the branches the remaining ~1% of the draws take
+* `zigguratWedge`: the rejection test on a strip that sticks out of the curve
+
+## References
+
+* G. Marsaglia and W. W. Tsang, *The Ziggurat Method for Generating Random Variables*, Journal of
+  Statistical Software, 2000.
+* numpy's samplers: `numpy/random/src/distributions/distributions.c`
+-/
+
+@[expose] public section
+
+namespace NumLean
+
+/-- The rejection test on a strip that sticks out of the curve: the point drawn at height `u`
+between the density at the strip's two edges lies under the curve. -/
+@[inline] def zigguratWedge (f : Array Float) (idx : Nat) (u density : Float) : Bool :=
+  Float.fma (f[idx - 1]! - f[idx]!) u f[idx]! < density
+
+/-- The rare branches of `ziggurat`, reached by about 1% of the draws: the strip either is the base
+one, whose unbounded part `tail` samples from the integer drawn, or sticks out of the curve, and
+then the point is tested against `density` and the whole draw is started over on rejection.
+
+`idx`, `ri` and `x` are the strip, the integer and the abscissa `ziggurat` has already drawn and
+found not to land in the rectangular part of its strip. Each redraw retries the fast path here
+rather than returning to `ziggurat`, so the two together run exactly the loop of numpy's samplers.
+
+This is kept apart from `ziggurat` so that the fast path can be inlined into its caller: a
+recursive function cannot be, and behind a call boundary every draw would have to box the generator
+state and its result, which costs several times the draw itself. -/
+@[specialize] partial def zigguratSlow (k : Array UInt64) (w f : Array Float)
+    (split : UInt64 → Nat × UInt64 × Bool) (density : Float → Float)
+    (tail : UInt64 → RandPCG IO Float) (idx : Nat) (ri : UInt64) (x : Float) :
+    RandPCG IO Float := do
+  if idx == 0 then tail ri
+  else if zigguratWedge f idx (← random) (density x) then return x
+  else
+    let (idx, ri, negate) := split (← randUInt64)
+    let x := ri.toFloat * w[idx]!
+    let x := if negate then -x else x
+    if ri < k[idx]! then return x
+    else zigguratSlow k w f split density tail idx ri x
+
+/-- The sampler shape shared by numpy's `random_standard_normal` and
+`random_standard_exponential`, the ziggurat of Marsaglia and Tsang: the density is covered by 256
+strips of equal area, and one 64-bit output supplies at once the strip `idx` and the integer `ri`
+that `w` scales to an abscissa, `split` saying how those bits are laid out and whether the deviate
+comes out negated.
+
+The draw is returned as it stands when `ri` falls below `k[idx]`, which is where about 99% of the
+draws end; `zigguratSlow` takes over the remaining ones. -/
+@[inline] def ziggurat (k : Array UInt64) (w f : Array Float)
+    (split : UInt64 → Nat × UInt64 × Bool) (density : Float → Float)
+    (tail : UInt64 → RandPCG IO Float) : RandPCG IO Float := do
+  let (idx, ri, negate) := split (← randUInt64)
+  let x := ri.toFloat * w[idx]!
+  let x := if negate then -x else x
+  if ri < k[idx]! then return x
+  else zigguratSlow k w f split density tail idx ri x
+
+end NumLean
diff --git a/RandomDo/Tactic/Computable/Counterparts.lean b/RandomDo/Tactic/Computable/Counterparts.lean
new file mode 100644
index 0000000..ee61542
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Counterparts.lean
@@ -0,0 +1,47 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import RandomDo.Tactic.Computable.Defs
+public import RandomDo.NumLean.Distributions
+public import Mathlib.Probability.Distributions.Gaussian.Real
+public import Mathlib.Probability.Distributions.Bernoulli
+
+/-!
+# Computable counterparts of the pieces an `rdo` program is made of
+
+The `@[computable]` attribute translates a program by replacing each piece it is made of by the
+counterpart recorded here through `@[computable_as]`. There is one entry per piece the programs of
+`RandomDo.Tactic.Computable.Example` are built from: the two types their values live in, and the
+one distribution they draw from.
+-/
+
+@[expose] public section
+
+/-! ## Types -/
+
+attribute [computable_as Float] Real
+attribute [computable_as Float] NNReal
+
+/-! ## Distributions -/
+
+attribute [computable_as NumLean.normal'] ProbabilityTheory.gaussianReal
+
+/-- Draw from a Bernoulli distribution with probability `p` of returning `x` and `1 - p` of
+returning `y`. -/
+def bernoulliChoice (α : Type) [_h : MeasurableSpace α] (x y : α) (p : Float) :
+    NumLean.RandPCG IO α := do
+  return if (← NumLean.bernoulli p) == 1 then x else y
+
+attribute [computable_as bernoulliChoice] ProbabilityTheory.bernoulliMeasure
+
+/-! ## Classical functions -/
+
+attribute [computable_as Float.sqrt] Real.sqrt
+
+attribute [computable_as Float.log] Real.log
+
+attribute [computable_as Float.exp] Real.exp
diff --git a/RandomDo/Tactic/Computable/Defs.lean b/RandomDo/Tactic/Computable/Defs.lean
new file mode 100644
index 0000000..66b5c81
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Defs.lean
@@ -0,0 +1,45 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public meta import Batteries.Lean.NameMapAttribute
+public meta import Lean.ReservedNameAction
+
+/-!
+# The `@[computable_as]` attribute
+
+An `rdo` program is written over the Giry monad: it draws from measures on a measurable space, which no machine samples. Turning it into a program that runs, which the `@[computable]` attribute of `RandomDo.Tactic.Computable.Deriving` does, asks for a counterpart of each piece the program is
+built from, e.g, `NumLean.normal'` for `gaussianReal`. This file holds the attribute
+recording them; `RandomDo.Tactic.Computable.Counterparts` holds the counterparts themselves.
+
+`@[computable_as f]` on a declaration `d` reads: `f` is what `d` becomes in a translated program.
+-/
+
+public meta section
+
+open Lean
+
+namespace RDo.Tactic
+
+/-- The counterparts recorded by `@[computable_as]`, keyed by the declaration they translate. -/
+initialize computableAsExt : NameMapExtension Name ←
+  registerNameMapAttribute {
+    name := `computable_as
+    descr := "record the computable counterpart of this declaration"
+    add := fun _ stx ↦ do
+      let f := stx[1][0]
+      unless f.isIdent do
+        throwError "`computable_as` takes the name of one declaration"
+      realizeGlobalConstNoOverload f
+  }
+
+/-- The computable counterpart of `declName`, when `@[computable_as]` recorded one. -/
+def computableAs? (declName : Name) : CoreM (Option Name) :=
+  return computableAsExt.find? (← getEnv) declName
+
+end RDo.Tactic
+
+end
diff --git a/RandomDo/Tactic/Computable/Deriving.lean b/RandomDo/Tactic/Computable/Deriving.lean
new file mode 100644
index 0000000..b1b1ce7
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Deriving.lean
@@ -0,0 +1,179 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public meta import RandomDo.Tactic.Computable.Counterparts
+public import RandomDo.Monad.MeasurableSpace
+public meta import Lean.Elab.Tactic.Basic
+
+/-!
+# The `@[computable]` attribute
+
+`@[computable]` reads an `rdo` program and adds the program that samples from it:
+
+```
+@[computable]
+noncomputable def shifted : Measure ℝ := rdo
+  let x ← gaussianReal 0 1
+  return x + 1
+```
+
+adds `shiftedComputable : RandPCG IO Float`, which draws from `NumLean.normal' 0 1` and adds one.
+
+The Giry monad and its two operations become `RandPCG IO`, `pure` and `bind`, the `for` loop of
+`rdo` becomes the `for` loop of that monad, and a `let` stays a `let`. Anything else is rebuilt
+from the counterpart `@[computable_as]` records for its head, with its arguments translated in turn
+and its instances synthesized anew. A term translates into a term of the translation of its
+type; where the rebuilt one does not, its head is a definition nothing is known about, and its body
+is read in its place. `@[computable]` records the program it writes, so a program drawing from
+another translates into one calling that other's translation.
+
+Two things extend the attribute: an `@[computable_as]` entry, and an alternative of `translate` for
+a construct of `rdo` it has not been taught.
+
+The counterparts are imported `meta` as well as publicly, and so reach every module the attribute
+does. A translated program is written to be run, and a `run_cmd` or an `#eval` runs it in the very
+module that writes it: the interpreter then asks for the code of what it calls, which a plain
+`import` does not carry.
+
+`set_option trace.computable true` prints what each piece became.
+-/
+
+public meta section
+
+open Lean Meta NumLean
+
+namespace RDo.Tactic
+
+initialize registerTraceClass `computable
+
+/-- The monad the translated programs live in. -/
+def computableMonad : MetaM Expr := mkAppM ``RandPCG #[mkConst ``IO]
+
+mutual
+
+/-- Translate a piece of an `rdo` program. `σ` maps the variables the program binds to the ones the
+translated program binds in their place. -/
+partial def translate (σ : FVarSubst) (e : Expr) : MetaM Expr :=
+  withTraceNode `computable (fun
+      | .ok e' => return m!"{e} ↦ {e'}"
+      | .error _ => return m!"{e}: not translated") do
+    match_expr e with
+    | MeasurableSpacePure.mPure _ _ _ _ a =>
+      mkAppOptM ``Pure.pure #[← computableMonad, none, none, ← translate σ a]
+    | MeasurableSpaceBind.mBind _ _ _ _ _ _ p k =>
+      mkAppOptM ``Bind.bind #[← computableMonad, none, none, none,
+        ← translate σ p, ← translate σ k]
+    | MeasurableSpaceForIn.forIn _ _ _ _ _ _ xs init body =>
+      mkAppOptM ``ForIn.forIn #[← computableMonad, none, none, none, none,
+        ← translate σ xs, ← translate σ init, ← translate σ body]
+    | MeasureTheory.Measure α _ => return mkApp (← computableMonad) (← translate σ α)
+    /- A subtype is its carrier, and one of its values is the value it carries: the constraint and
+    the proof of it are what a computable counterpart does not have. -/
+    | Subtype α _ => translate σ α
+    | Subtype.mk _ _ v _ => translate σ v
+    | _ => match e with
+      | .fvar x => return σ.get x
+      | .sort .. | .lit .. => return e
+      | .mdata _ b => translate σ b
+      | .lam .. => lambdaBoundedTelescope e 1 fun xs body ↦ do
+        let x := xs[0]!.fvarId!
+        withLocalDeclD (← x.getUserName) (← translate σ (← x.getType)) fun y ↦ do
+          mkLambdaFVars #[y] (← translate (σ.insert x y) body)
+      | .forallE .. => forallBoundedTelescope e (some 1) fun xs body ↦ do
+        let x := xs[0]!.fvarId!
+        withLocalDeclD (← x.getUserName) (← translate σ (← x.getType)) fun y ↦ do
+          mkForallFVars #[y] (← translate (σ.insert x y) body)
+      | .letE n t v b _ => withLetDecl n t v fun x ↦ do
+        withLetDecl n (← translate σ t) (← translate σ v) fun y ↦ do
+          mkLetFVars #[y] (← translate (σ.insert x.fvarId! y) (b.instantiate1 x))
+      | _ => translateApp σ e
+
+/-- Rebuild an application from the counterpart of its head; where nothing known about that head
+fits, look through it and read its body in its place. -/
+partial def translateApp (σ : FVarSubst) (e : Expr) : MetaM Expr := do
+  if e.getAppFn.isFVar then
+    return mkAppN (← translate σ e.getAppFn) (← e.getAppArgs.mapM (translate σ))
+  if e.getAppFn.isLambda then return ← translate σ e.headBeta
+  let .const declName _ := e.getAppFn | throwError "`computable`: cannot translate{indentExpr e}"
+  let counterpart? ← computableAs? declName
+  let head := counterpart?.getD declName
+  /- We save the state of metavariables, so that a failed rebuild does not leave them in a
+  half-built state. -/
+  let s ← saveState
+  try
+    /- A term translates into a term of the translation of its type. A head nothing is known about
+    rebuilds into itself, of the type it had, and might fail. We check that the rebuilt term has
+    the translation of the type, and if not, we look through it. -/
+    let f ← rebuild σ head e
+    -- Useless for a type, whose type is a sort either way.
+    if (← inferType f).isSort then return f
+    let expected ← translate σ (← inferType e)
+    unless ← isDefEq (← inferType f) expected do
+      throwError "`computable`: {f} is of type{indentExpr (← inferType f)}\n\
+        where the translation asks for{indentExpr expected}"
+    return f
+  catch ex =>
+    restoreState s
+    /- The rebuild failed, we try to look through the head, and read its body in its place. -/
+    if counterpart?.isNone then
+      if let some e' ← unfoldDefinition? e then
+        trace[computable] "nothing known about {declName}, looking through it"
+        return ← translate σ e'
+    throw ex
+
+/-- Apply `head` to the arguments of `e`, translated in turn, and check that the result has the
+translation of the type of `e`. -/
+partial def rebuild (σ : FVarSubst) (head : Name) (e : Expr) : MetaM Expr := do
+  let mut f ← mkConstWithFreshMVarLevels head
+  for arg in e.getAppArgs do
+    let .forallE _ t _ bi ← whnf (← inferType f)
+      | throwError "`computable`: {f} does not take the argument{indentExpr arg}"
+    let arg ←
+      if bi.isInstImplicit then do
+        -- An instance is not translated: it is asked for anew, at the translated types.
+        let inst ← synthInstance t
+        trace[computable] "instance: {inst}"
+        pure inst
+      else
+        translate σ arg
+    unless ← isDefEq (← inferType arg) t do
+      throwError "`computable`: {arg} does not fit the argument of {f}, of type{indentExpr t}"
+    f := mkApp f arg
+  return f
+
+end
+
+/-- Translate the `rdo` program `declName` and add the translation to the environment, under the
+name `declName` followed by `Computable`. -/
+def addComputableDecl (declName : Name) : MetaM Unit := do
+  let info ← getConstInfo declName
+  let some value := info.value?
+    | throwError "`computable` can only be derived for a definition, but {declName} has no value"
+  let value ← instantiateMVars (← translate {} value)
+  let type ← instantiateMVars (← inferType value)
+  let translated := declName.appendAfter "Computable"
+  addAndCompile <| .defnDecl <| ← mkDefinitionValInferringUnsafe translated info.levelParams type
+    value (.regular (getMaxHeight (← getEnv) value + 1))
+  addDocStringCore translated s!"The computable program that samples from `{declName}` \
+    (automatically generated by the `@[computable]` attribute)."
+  computableAsExt.add declName translated
+  trace[computable] "wrote {translated}:{indentExpr type}"
+
+@[inherit_doc addComputableDecl]
+initialize registerBuiltinAttribute {
+  name := `computable
+  descr := "translate this `rdo` program into the program that samples from it"
+  applicationTime := .afterCompilation
+  add := fun declName _stx kind ↦ do
+    unless kind == AttributeKind.global do
+      throwError "`computable` must be a global attribute"
+    (addComputableDecl declName).run'
+}
+
+end RDo.Tactic
+
+end
diff --git a/RandomDo/Tactic/Computable/Example.lean b/RandomDo/Tactic/Computable/Example.lean
new file mode 100644
index 0000000..13f7e54
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Example.lean
@@ -0,0 +1,40 @@
+module
+/- /-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import RandomDo.Monad.Instances
+public import RandomDo.Monad.Notation
+public import RandomDo.NumLean.Distributions
+public import RandomDo.Tactic.Computable.Deriving
+public import Mathlib
+meta import RandomDo.NumLean.Distributions
+meta import Batteries.Data.Float.Basic
+
+/-!
+# `@[computable]` on an `rdo` program
+
+`test` denotes a distribution: a Gaussian draw, shifted by one. The attribute reads it and writes
+the program that samples from it, drawing from `NumLean.normal` instead and shifting the draw the
+same way.
+-/
+
+@[expose] public section
+
+open MeasureTheory ProbabilityTheory NumLean
+
+/-- Test -/
+@[computable]
+noncomputable def test (m : ℝ) : Measure ℝ := rdo
+  let x ← gaussianReal m 1
+  return x + 1
+
+/- And it runs: seeded alike, it draws what numpy 2.3.4 draws from `default_rng(42).normal() + 1`,
+down to the last bit. -/
+run_cmd do IO.runRandPCG do
+  let x ← (IO.runRandPCG <| test_computable 10 : IO Float)
+  Lean.logInfo m!"`test_computable` drew {x.toStringFull}"
+ -/
diff --git a/RandomDo/Tactic/Computable/Polymorphic/Polymorphic.lean b/RandomDo/Tactic/Computable/Polymorphic/Polymorphic.lean
new file mode 100644
index 0000000..963ccb1
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Polymorphic/Polymorphic.lean
@@ -0,0 +1,98 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import Mathlib.Probability.Distributions.Gaussian.Real
+public import Mathlib.Probability.Distributions.Bernoulli
+public import Mathlib.MeasureTheory.Function.SpecialFunctions.Basic
+public import RandomDo.Monad.Instances
+public import RandomDo.NumLean.Distributions
+
+/-!
+# Polymorphic `rdo` programs
+
+A program written over an arbitrary `MeasurableSpaceMonad` `m`, drawing through the classes of this
+file, is read at `m := Measure` to prove things about it and run at `m := RandM` to sample from it.
+Each class has an instance of each kind: the distribution of Mathlib on `ℝ`, and the sampler of
+`NumLean` on `Float`. The scalar classes `HasExp`, `HasLog` and `HasSqrt` do the same for the
+functions a program computes with.
+-/
+
+@[expose] public section
+
+open MeasureTheory ProbabilityTheory NumLean
+
+universe v
+
+/-- A typeclass for monads that can draw from a Gaussian distribution. -/
+class HasGaussian (m : (α : Type) → [MeasurableSpace α] → Type v)
+   (α β : Type*) (R : Type) [MeasurableSpace R] where
+  /-- Draw a sample from a Gaussian distribution with mean `μ` and variance `v`. -/
+  gaussian : α → β → m R
+
+noncomputable instance : HasGaussian Measure ℝ NNReal ℝ where
+  gaussian μ v := gaussianReal μ v
+
+/-- The monad that samples, seen as a `MeasurableSpaceMonad`. -/
+abbrev RandM := Monad.toMeasurableSpaceMonad (RandPCG IO)
+
+instance instMeasurableSpaceFloat : MeasurableSpace Float := ⊤
+
+instance : HasGaussian RandM Float Float Float where
+  gaussian μ v := normal' μ v
+
+/-- A typeclass for monads that can draw from a Bernoulli distribution. -/
+class HasBernoulli (m : (α : Type) → [MeasurableSpace α] → Type v) (P : Type) where
+  /-- Draw `true` with probability `p`, and `false` otherwise. -/
+  bernoulli : P → m Bool
+
+/-- A probability outside `[0, 1]` is clamped to it. -/
+noncomputable instance : HasBernoulli Measure ℝ where
+  bernoulli p := bernoulliMeasure true false (Set.projIcc 0 1 zero_le_one p)
+
+/-- A probability outside `[0, 1]` is clamped to it, as in the instance on `Measure`. -/
+instance : HasBernoulli RandM Float where
+  bernoulli p := return (← NumLean.bernoulli p) == 1
+
+/-- A typeclass for scalars with an exponential. -/
+class HasExp (R : Type) where
+  /-- The exponential. -/
+  exp : R → R
+
+noncomputable instance : HasExp ℝ := ⟨Real.exp⟩
+
+@[fun_prop]
+lemma HasExp.measurable_exp_real : Measurable (HasExp.exp : ℝ → ℝ) := Real.measurable_exp
+
+instance : HasExp Float := ⟨Float.exp⟩
+
+/-- A typeclass for scalars with a logarithm. -/
+class HasLog (R : Type) where
+  /-- The logarithm. -/
+  log : R → R
+
+noncomputable instance : HasLog ℝ := ⟨Real.log⟩
+
+@[fun_prop]
+lemma HasLog.measurable_log_real : Measurable (HasLog.log : ℝ → ℝ) := Real.measurable_log
+
+instance : HasLog Float := ⟨Float.log⟩
+
+/-- A typeclass for scalars with a square root. -/
+class HasSqrt (R : Type) where
+  /-- The square root. -/
+  sqrt : R → R
+
+noncomputable instance : HasSqrt ℝ := ⟨Real.sqrt⟩
+
+@[fun_prop]
+lemma HasSqrt.measurable_sqrt_real : Measurable (HasSqrt.sqrt : ℝ → ℝ) :=
+  Real.continuous_sqrt.measurable
+
+instance : HasSqrt Float := ⟨Float.sqrt⟩
+
+/-- A natural number as a `Float`, so that programs polymorphic in the scalars can cast counts. -/
+instance instNatCastFloat : NatCast Float := ⟨Nat.toFloat⟩
diff --git a/RandomDo/Tactic/Computable/Polymorphic/Scalar.lean b/RandomDo/Tactic/Computable/Polymorphic/Scalar.lean
new file mode 100644
index 0000000..30cfe3b
--- /dev/null
+++ b/RandomDo/Tactic/Computable/Polymorphic/Scalar.lean
@@ -0,0 +1,22 @@
+/-
+Copyright (c) 2026 Gaëtan Serré. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré
+-/
+module
+
+public import Mathlib.MeasureTheory.MeasurableSpace.Defs
+
+/-!
+# Polymorphic scalar classes
+
+A typeclass for measurable spaces with the operations and numerals of a field, but none of its
+axioms, so that it applies to both `ℝ` and `Float`.
+-/
+
+@[expose] public section
+
+/-- `Scalar R` bundles a measurable space structure on `R` with the operations and numerals of a
+field, but none of its axioms, so that it applies to both `ℝ` and `Float`. -/
+class abbrev Scalar (R : Type*) :=
+  MeasurableSpace R, Add R, Mul R, Div R, Sub R, Zero R, One R, NatCast R, OfScientific R
diff --git a/RandomDo/Tactic/IsMarkov.lean b/RandomDo/Tactic/IsMarkov/Defs.lean
similarity index 100%
rename from RandomDo/Tactic/IsMarkov.lean
rename to RandomDo/Tactic/IsMarkov/Defs.lean
diff --git a/RandomDo/Tactic/Deriving.lean b/RandomDo/Tactic/IsMarkov/Deriving.lean
similarity index 97%
rename from RandomDo/Tactic/Deriving.lean
rename to RandomDo/Tactic/IsMarkov/Deriving.lean
index 6129c23..2c90e5e 100644
--- a/RandomDo/Tactic/Deriving.lean
+++ b/RandomDo/Tactic/IsMarkov/Deriving.lean
@@ -5,7 +5,7 @@ Authors: Rémy Degenne
 -/
 module
 
-public import RandomDo.Tactic.Elab
+public import RandomDo.Tactic.IsMarkov.Elab
 
 /-!
 # The `@[is_markov]` attribute
@@ -63,7 +63,7 @@ def addIsMarkovInstance (declName : Name) : TermElabM Unit := do
                       value := ← instantiateMVars proof })
   Meta.addInstance instName .global 1000
 
-/-- The `@[is_markov]` attribute. -/
+@[inherit_doc isMarkovStatement]
 initialize registerBuiltinAttribute {
   name := `is_markov
   descr := "prove that this `rdo` program is a Markov kernel, and register it as an instance"
diff --git a/RandomDo/Tactic/Elab.lean b/RandomDo/Tactic/IsMarkov/Elab.lean
similarity index 96%
rename from RandomDo/Tactic/Elab.lean
rename to RandomDo/Tactic/IsMarkov/Elab.lean
index 4f7c8ce..fa013d8 100644
--- a/RandomDo/Tactic/Elab.lean
+++ b/RandomDo/Tactic/IsMarkov/Elab.lean
@@ -5,7 +5,7 @@ Authors: Gaëtan Serré
 -/
 module
 
-public import RandomDo.Tactic.Lemmas
+public import RandomDo.Tactic.IsMarkov.Lemmas
 public meta import Lean.Elab.Tactic.Basic
 
 /-!
@@ -136,6 +136,9 @@ def closeLeaf (g : MVarId) : MetaM (List MVarId) := do
   if let some gs ← observing? (g.applyConst ``IsMarkov.gaussianReal) then
     trace[is_markov] "`gaussianReal` leaf: handing back the measurability of its parameters"
     return gs
+  if let some gs ← observing? (g.applyConst ``IsMarkov.bernoulliMeasure) then
+    trace[is_markov] "`bernoulliMeasure` leaf: handing back the measurability of its parameters"
+    return gs
   return [g]
 
 /-- The constant heading the body of `κ`, when it is a definition the tactic could look through. -/
@@ -172,7 +175,9 @@ def abstractLoopVars (vars : Array FVarId) (g : MVarId) : MetaM MVarId := do
 
 /-- Turn a goal `IsMarkov κ` into the list of goals the user is left with. -/
 partial def isMarkovCore (g : MVarId) (fuel : Nat) : MetaM (List MVarId) := g.withContext do
-  let target ← instantiateMVars (← g.getType)
+  /- The annotations a goal may carry, e.g. the one a tactic `have` leaves on the goal after it,
+  would hide the head of the statement. -/
+  let target := (← instantiateMVars (← g.getType)).cleanupAnnotations
   -- `IsMarkov` takes five arguments: `γ`, `α`, their `MeasurableSpace` instances, and `κ`.
   unless target.isAppOfArity ``IsMarkov 5 do
     trace[is_markov] "not an `IsMarkov` goal, handed back: {target}"
@@ -251,7 +256,8 @@ partial def isMarkovCore (g : MVarId) (fuel : Nat) : MetaM (List MVarId) := g.wi
         let mut goals := []
         let mut side := []
         for g' in gs do
-          if (← instantiateMVars (← g'.getType)).isAppOfArity ``IsMarkov 5 then
+          let t := (← instantiateMVars (← g'.getType)).cleanupAnnotations
+          if t.isAppOfArity ``IsMarkov 5 then
             goals := goals ++ (← isMarkovCore g' fuel)
           else
             side := side ++ [g']
@@ -280,7 +286,7 @@ partial def isMarkovCore (g : MVarId) (fuel : Nat) : MetaM (List MVarId) := g.wi
       if ← g.isAssigned then return leftover
       /- The goal was not closed, so we try to unfold names in the head of the program until we
       reach a known shape. If that fails, we leave the goal to the user. -/
-      match ← unfoldToKnownShape (← instantiateMVars (← g.getType)) fuel with
+      match ← unfoldToKnownShape (← instantiateMVars (← g.getType)).cleanupAnnotations fuel with
       | some target =>
         trace[is_markov] "unfolded the head definition to: {target.appArg!}"
         isMarkovCore (← g.change target) (fuel - 1)
@@ -316,7 +322,7 @@ lemma _root_.isProbabilityMeasure_of_isMarkov {α : Type*} [MeasurableSpace α]
 /-- Bring a goal of the form `IsProbabilityMeasure μ` into the form `IsMarkov fun _ : Unit ↦ μ`, so
 that `isMarkovCore` can be applied. -/
 def toIsMarkovGoal (g : MVarId) : MetaM MVarId := do
-  let target ← instantiateMVars (← g.getType)
+  let target := (← instantiateMVars (← g.getType)).cleanupAnnotations
   unless target.isAppOfArity ``IsProbabilityMeasure 3 do
     return g
   match ← g.applyConst ``isProbabilityMeasure_of_isMarkov with
diff --git a/RandomDo/Tactic/ForInStep.lean b/RandomDo/Tactic/IsMarkov/ForInStep.lean
similarity index 98%
rename from RandomDo/Tactic/ForInStep.lean
rename to RandomDo/Tactic/IsMarkov/ForInStep.lean
index e1cce94..6351f45 100644
--- a/RandomDo/Tactic/ForInStep.lean
+++ b/RandomDo/Tactic/IsMarkov/ForInStep.lean
@@ -5,7 +5,7 @@ Authors: Gaëtan Serré
 -/
 module
 
-public import RandomDo.Tactic.IsMarkov
+public import RandomDo.Tactic.IsMarkov.Defs
 public import RandomDo.Monad.MeasurableSpace
 
 /-!
diff --git a/RandomDo/Tactic/Lemmas.lean b/RandomDo/Tactic/IsMarkov/Lemmas.lean
similarity index 92%
rename from RandomDo/Tactic/Lemmas.lean
rename to RandomDo/Tactic/IsMarkov/Lemmas.lean
index 1284d31..38cc0d0 100644
--- a/RandomDo/Tactic/Lemmas.lean
+++ b/RandomDo/Tactic/IsMarkov/Lemmas.lean
@@ -8,10 +8,11 @@ module
 public import RandomDo.Monad.Instances
 public import RandomDo.Monad.ForInInstances
 public import RandomDo.Measurable
-public import RandomDo.Tactic.ForInStep
+public import RandomDo.Tactic.IsMarkov.ForInStep
 public import Mathlib.MeasureTheory.Measure.ProbabilityMeasure
 public import Mathlib.Data.List.OfFn
 public import Mathlib.Probability.Distributions.Gaussian.Real
+public import Mathlib.Probability.Distributions.Bernoulli
 
 /-!
 # Markov property of `rdo` programs
@@ -31,6 +32,8 @@ complex program to the Markov property/measurability of its underlying mathemati
   the bound variable is Markovian in the parameter.
 * `gaussianReal`: A Gaussian distribution whose mean and variance depend measurably on the parameter
   is Markovian in the parameter.
+* `bernoulliMeasure`: A Bernoulli distribution whose two outcomes and probability depend measurably
+  on the parameter is Markovian in the parameter.
 * `comp`: Composing a Markov kernel `κ` with a measurable function `g` is Markovian in the
   parameter.
 * `ite`: A conditional `rdo` program that chooses between two Markov kernels `κ` and `η` based on a
@@ -43,6 +46,7 @@ complex program to the Markov property/measurability of its underlying mathemati
 * `forInList_comp`, `forInArray_comp`, `forInVector_comp`: The same three, for a loop over a
   collection the program takes as an argument. The body is then asked to be Markovian jointly in the
   parameter and in the element, which the fixed collections do not need.
+* `forIn_nil`, `forIn_cons`: A `for` loop over a list, unrolled one element at a time.
 * `breakRunK`: The case analysis a program performs after a loop that returns early, on the `Option`
   slot holding the returned value, is Markovian as soon as both of its branches are.
 -/
@@ -51,6 +55,7 @@ complex program to the Markov property/measurability of its underlying mathemati
 
 open MeasureTheory ProbabilityTheory Function
 open MeasurableSpacePure
+open scoped ENNReal
 
 namespace IsMarkov
 
@@ -91,6 +96,18 @@ lemma gaussianReal {m : γ → ℝ} {v : γ → NNReal} (hm : Measurable m) (hv
     IsMarkov fun c ↦ ProbabilityTheory.gaussianReal (m c) (v c) :=
   ⟨ProbabilityTheory.measurable_gaussianReal.comp (hm.prodMk hv), fun _ ↦ inferInstance⟩
 
+lemma bernoulliMeasure {x y : γ → α} {p : γ → unitInterval} (hx : Measurable x)
+    (hy : Measurable y) (hp : Measurable p) :
+    IsMarkov fun c ↦ ProbabilityTheory.bernoulliMeasure (x c) (y c) (p c) := by
+  refine ⟨Measure.measurable_of_measurable_coe _ fun s hs ↦ ?_, fun _ ↦ inferInstance⟩
+  simp only [bernoulliMeasure_def, Measure.add_apply, Measure.smul_apply,
+    Measure.dirac_apply' _ hs, ENNReal.smul_def, smul_eq_mul]
+  have hp' : Measurable fun c ↦ ((unitInterval.toNNReal (p c) : ℝ≥0∞)) := by fun_prop
+  have hq' : Measurable fun c ↦ ((unitInterval.toNNReal (unitInterval.symm (p c)) : ℝ≥0∞)) := by
+    fun_prop
+  exact (hp'.mul ((measurable_one.indicator hs).comp hx)).add
+    (hq'.mul ((measurable_one.indicator hs).comp hy))
+
 lemma comp {κ : γ → Measure α} (hκ : IsMarkov κ) {g : σ → γ} (hg : Measurable g) :
     IsMarkov fun c ↦ κ (g c) := ⟨hκ.measurable.comp hg, fun _ ↦ hκ.isProbabilityMeasure _⟩
 
@@ -228,11 +245,14 @@ private lemma forIn_eq_listLoop (l : List ι) (b : σ) (g : ι → σ → Measur
     MeasurableSpaceForIn.forIn (m := Measure) l b g = listLoop g l b :=
   loop_eq_listLoop g l b l _ (fun _ _ _ ↦ rfl) ⟨[], rfl⟩
 
-private lemma forIn_nil (b : σ) (g : ι → σ → Measure (ForInStep σ)) :
+/-- A `for` loop over the empty list returns its initial state. -/
+lemma forIn_nil (b : σ) (g : ι → σ → Measure (ForInStep σ)) :
     MeasurableSpaceForIn.forIn (m := Measure) ([] : List ι) b g = mPure b :=
   forIn_eq_listLoop _ _ _
 
-private lemma forIn_cons (a : ι) (l : List ι) (b : σ) (g : ι → σ → Measure (ForInStep σ)) :
+/-- A `for` loop over `a :: l` runs its body on `a`, then stops or carries on with the loop over
+`l`. -/
+lemma forIn_cons (a : ι) (l : List ι) (b : σ) (g : ι → σ → Measure (ForInStep σ)) :
     MeasurableSpaceForIn.forIn (m := Measure) (a :: l) b g
       = g a b >>=ₘ fun step ↦ ForInStep.casesOn (motive := fun _ ↦ Measure σ) step mPure
           fun b' ↦ MeasurableSpaceForIn.forIn (m := Measure) l b' g := by
diff --git a/Test.lean b/Test.lean
index 7862fd2..bba557b 100644
--- a/Test.lean
+++ b/Test.lean
@@ -2,9 +2,12 @@ module  -- shake: keep-all --deprecated_module: ignore
 
 public import Test.Bind
 public import Test.Common
+public import Test.Computable
 public import Test.Control
 public import Test.Gaps
 public import Test.Instances
 public import Test.IsMarkov
 public import Test.Loops
 public import Test.MonadLaws
+public import Test.Polymorphic
+public import Test.Sample
diff --git a/Test/Computable.lean b/Test/Computable.lean
new file mode 100644
index 0000000..0b7752b
--- /dev/null
+++ b/Test/Computable.lean
@@ -0,0 +1,62 @@
+module
+
+public import Test.IsMarkov
+public import Test.Bind
+public meta import RandomDo
+import Batteries.Data.Float.Basic
+
+set_option linter.style.header false
+
+set_option trace.computable true
+
+namespace Test.Computable
+
+open Test.IsMarkov NumLean Lean.Elab.Command MeasureTheory ProbabilityTheory
+
+def logComputable {α : Type} [Lean.ToMessageData α] (prog : RandPCG IO α) : CommandElabM Unit := do
+  let x ← (IO.runRandPCG prog : IO α)
+  let y ← (IO.runRandPCGWith 42 prog : IO α)
+  Lean.logInfo m!"x = {x}"
+  Lean.logInfo m!"y (seed 42) = {y}"
+
+attribute [computable] sumTwo
+
+run_cmd logComputable sumTwoComputable
+
+@[computable]
+noncomputable
+def unfoldSumTwo : Measure ℝ := rdo
+  let y ← sumTwo
+  let x ← gaussianReal 0 1
+  return x + y
+
+attribute [computable] centred
+
+run_cmd logComputable (centredComputable 20)
+
+attribute [computable] branchOn
+
+run_cmd do logComputable (branchOnComputable 20)
+
+run_cmd logComputable (branchOnComputable (-1))
+
+attribute [computable] fairCoin
+
+run_cmd logComputable (fairCoinComputable)
+
+attribute [computable] Bind.twoCoins
+
+run_cmd logComputable (Bind.twoCoinsComputable)
+
+@[computable]
+noncomputable
+def ex1 : Measure ℝ := rdo
+  let mut x := 0
+  for _ in List.range 1000 rdo
+    let y ← gaussianReal 0 1
+    x := x + y
+  return x
+
+run_cmd logComputable (ex1Computable)
+
+end Test.Computable
diff --git a/Test/Polymorphic.lean b/Test/Polymorphic.lean
new file mode 100644
index 0000000..566bf7d
--- /dev/null
+++ b/Test/Polymorphic.lean
@@ -0,0 +1,119 @@
+module
+
+public import Test.IsMarkov
+public meta import RandomDo
+import Batteries.Data.Float.Basic
+
+set_option linter.style.header false
+
+/-!
+# Polymorphic `rdo` programs
+
+The programs of `Test.Computable`, written once over an arbitrary `MeasurableSpaceMonad` `m` and
+drawing through `HasGaussian` and `HasBernoulli`. Read at `m := Measure`, each one is a probability
+measure, checked by `is_markov`, and is the program of `Test.IsMarkov` when there is one. Run at
+`m := RandM`, it samples.
+-/
+
+@[expose] public section
+
+namespace Test.Polymorphic
+
+open Test.IsMarkov NumLean Lean.Elab.Command MeasureTheory ProbabilityTheory
+
+universe v
+
+variable {m : (α : Type) → [MeasurableSpace α] → Type v} [MeasurableSpaceMonad m]
+  {R V : Type} [Scalar R] [Scalar V]
+
+def logPolymorphic {α : Type} [MeasurableSpace α] [Lean.ToMessageData α] (prog : RandM α) :
+    CommandElabM Unit := do
+  let x ← (IO.runRandPCG prog : IO α)
+  let y ← (IO.runRandPCGWith 42 prog : IO α)
+  Lean.logInfo m!"x = {x}"
+  Lean.logInfo m!"y (seed 42) = {y}"
+
+def sumTwo [HasGaussian m R V R] : m R := rdo
+  let x ← HasGaussian.gaussian (m := m) (0 : R) (1 : V)
+  let y ← HasGaussian.gaussian (m := m) (0 : R) (1 : V)
+  return x + y
+
+example : sumTwo (m := Measure) (R := ℝ) (V := NNReal) = Test.IsMarkov.sumTwo := rfl
+
+run_cmd logPolymorphic (sumTwo (m := RandM) (R := Float) (V := Float))
+
+def unfoldSumTwo [HasGaussian m R V R] : m R := rdo
+  let y ← sumTwo (m := m) (V := V)
+  let x ← HasGaussian.gaussian (m := m) (0 : R) (1 : V)
+  return x + y
+
+example : IsProbabilityMeasure (unfoldSumTwo (m := Measure) (R := ℝ) (V := NNReal)) := by
+  is_markov
+
+run_cmd logPolymorphic (unfoldSumTwo (m := RandM) (R := Float) (V := Float))
+
+def centred [HasGaussian m R V R] (c : R) : m R := rdo
+  let x ← HasGaussian.gaussian (m := m) c (1 : V)
+  return x
+
+example : centred (m := Measure) (R := ℝ) (V := NNReal) = Test.IsMarkov.centred := rfl
+
+run_cmd logPolymorphic (centred (m := RandM) (R := Float) (V := Float) 20)
+
+def branchOn [LT R] [DecidableLT R] [HasGaussian m R V R] (c : R) : m R := rdo
+  if 0 < c then
+    let x ← HasGaussian.gaussian (m := m) c (1 : V)
+    return x
+  else
+    let x ← HasGaussian.gaussian (m := m) (0 : R) (1 : V)
+    return x
+
+example : branchOn (m := Measure) (R := ℝ) (V := NNReal) = Test.IsMarkov.branchOn := rfl
+
+run_cmd logPolymorphic (branchOn (m := RandM) (R := Float) (V := Float) 20)
+
+run_cmd logPolymorphic (branchOn (m := RandM) (R := Float) (V := Float) (-1))
+
+def coin [HasBernoulli m R] (p : R) : m Bool := rdo
+  let b ← HasBernoulli.bernoulli (m := m) p
+  return b
+
+example : IsProbabilityMeasure (coin (m := Measure) (1 / 2 : ℝ)) := by is_markov
+
+run_cmd logPolymorphic (coin (m := RandM) (0.5 : Float))
+
+/--
+error: expected 0 ≤ p ≤ 1, got 2.000000
+-/
+#guard_msgs in
+run_cmd logPolymorphic (coin (m := RandM) (2 : Float))
+
+/--
+error: expected 0 ≤ p ≤ 1, got -1.000000
+-/
+#guard_msgs in
+run_cmd logPolymorphic (coin (m := RandM) (-1 : Float))
+
+def twoCoins [HasBernoulli m R] (p : R) : m Bool := rdo
+  let x ← coin (m := m) p
+  let y ← coin (m := m) p
+  return x && y
+
+example : IsProbabilityMeasure (twoCoins (m := Measure) (1 / 2 : ℝ)) := by is_markov
+
+run_cmd logPolymorphic (twoCoins (m := RandM) (0.5 : Float))
+
+def ex1 [HasGaussian m R V R] : m R := rdo
+  let mut x : R := 0
+  for _ in List.range 1000 rdo
+    let y ← HasGaussian.gaussian (m := m) (0 : R) (1 : V)
+    x := x + y
+  return x
+
+example : IsProbabilityMeasure (ex1 (m := Measure) (R := ℝ) (V := NNReal)) := by is_markov
+
+run_cmd logPolymorphic (ex1 (m := RandM) (R := Float) (V := Float))
+
+end Test.Polymorphic
+
+end
diff --git a/Test/Sample.lean b/Test/Sample.lean
new file mode 100644
index 0000000..1d1a133
--- /dev/null
+++ b/Test/Sample.lean
@@ -0,0 +1,31 @@
+module
+
+public import RandomDo
+
+set_option linter.style.header false
+
+open NumLean IO Lean
+
+run_cmd do
+  let x ← rand 0 1000
+  logInfo m!"Random number between 0 and 1000: {x}"
+
+run_cmd do
+  let x ← (runRandPCG <| randInt 1000 : IO Int)
+  logInfo m!"Random number between 0 and 1000 (PCG): {x}"
+
+run_cmd do
+  let x ← (runRandPCG <| normal 10 2 : IO Float)
+  logInfo m!"{x}"
+
+run_cmd do
+  let x ← (runRandPCG <| exponential 10 : IO Float)
+  logInfo m!"{x}"
+
+run_cmd do
+  let x ← (runRandPCG <| binomial 10 0.5 : IO Nat)
+  logInfo m!"{x}"
+
+run_cmd do
+  let x ← (runRandPCG <| bernoulli 0.5 : IO Nat)
+  logInfo m!"{x}"
diff --git a/lakefile.toml b/lakefile.toml
index f51bc95..1dd6684 100644
--- a/lakefile.toml
+++ b/lakefile.toml
@@ -14,8 +14,18 @@ name = "LeanMachineLearning"
 git = "https://github.com/LeanMachineLearning/LML"
 rev = "main"
 
+[[lean_lib]]
+name = "FFI"
+precompileModules = true
+
 [[lean_lib]]
 name = "RandomDo"
 
 [[lean_lib]]
 name = "Test"
+
+# Used to run the tests in `scripts`
+[[lean_exe]]
+name = "dump"
+root = "Dump"
+srcDir = "test_data"
diff --git a/scripts/check_all.sh b/scripts/check_all.sh
new file mode 100755
index 0000000..bdf06af
--- /dev/null
+++ b/scripts/check_all.sh
@@ -0,0 +1,11 @@
+#!/usr/bin/env bash
+#
+
+cd "$(dirname "$0")/.." || exit 1
+
+status=0
+for script in scripts/check_*.py; do
+  echo "== $script"
+  python3 "$script" || status=1
+done
+exit $status
diff --git a/scripts/check_binomial.py b/scripts/check_binomial.py
new file mode 100644
index 0000000..a26caf1
--- /dev/null
+++ b/scripts/check_binomial.py
@@ -0,0 +1,27 @@
+import numpy as np
+from common import compare
+
+print("Checking binomial distribution...")
+
+# One pair per branch of numpy's `random_binomial`: inversion and BTPE, each on both sides of
+# `p = 1/2`, on both sides of the `n * p = 30` threshold, and the three degenerate cases.
+PARAMS = [(10, 0.3), (60, 0.5), (100, 0.31), (1000, 0.5), (1000, 0.9), (5, 0.99),
+          (0, 0.5), (10, 0.0), (10, 1.0)]
+
+LEAN = """import RandomDo
+def params : List (Nat × Float) :=
+  [%s]
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range (@N@ / %d) do
+          for (n, p) in params do
+            h.putStrLn (toString (← NumLean.binomial n p))
+""" % (",\n   ".join(f"({n}, {p})" for n, p in PARAMS), len(PARAMS))
+
+def distrib(seed, N):
+    g = np.random.default_rng(seed)
+    return [g.binomial(n, p) for _ in range(N // len(PARAMS)) for n, p in PARAMS]
+
+compare(LEAN, distrib)
diff --git a/scripts/check_exponential.py b/scripts/check_exponential.py
new file mode 100644
index 0000000..77f14e0
--- /dev/null
+++ b/scripts/check_exponential.py
@@ -0,0 +1,17 @@
+import numpy as np
+from common import compare
+
+print("Checking exponential distribution...")
+
+LEAN = """import RandomDo
+import Batteries.Data.Float.Basic
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range @N@ do h.putStrLn (← NumLean.exponential).toStringFull
+"""
+
+distrib = lambda seed, N: np.random.default_rng(seed).exponential(size=N)
+
+compare(LEAN, distrib)
\ No newline at end of file
diff --git a/scripts/check_float.py b/scripts/check_float.py
new file mode 100755
index 0000000..c6eb55b
--- /dev/null
+++ b/scripts/check_float.py
@@ -0,0 +1,17 @@
+import numpy as np
+from common import compare
+
+print("Checking uniform distribution over [0, 1)...")
+
+LEAN = """import RandomDo
+import Batteries.Data.Float.Basic
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range @N@ do h.putStrLn (← NumLean.random).toStringFull
+"""
+
+distrib = lambda seed, N: np.random.default_rng(seed).random(N)
+
+compare(LEAN, distrib)
\ No newline at end of file
diff --git a/scripts/check_int.py b/scripts/check_int.py
new file mode 100755
index 0000000..8cb7dc6
--- /dev/null
+++ b/scripts/check_int.py
@@ -0,0 +1,17 @@
+import numpy as np
+from common import compare
+
+print("Checking uniform distribution over integers...")
+
+LEAN = """import RandomDo
+import Batteries.Data.Float.Basic
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range @N@ do h.putStrLn <| toString (← NumLean.randInt 1000000)
+"""
+
+distrib = lambda seed, N: np.random.default_rng(seed).integers(1000000, size=N)
+
+compare(LEAN, distrib)
\ No newline at end of file
diff --git a/scripts/check_normal.py b/scripts/check_normal.py
new file mode 100644
index 0000000..366f44e
--- /dev/null
+++ b/scripts/check_normal.py
@@ -0,0 +1,17 @@
+import numpy as np
+from common import compare
+
+print("Checking normal distribution...")
+
+LEAN = """import RandomDo
+import Batteries.Data.Float.Basic
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range @N@ do h.putStrLn (← NumLean.normal).toStringFull
+"""
+
+distrib = lambda seed, N: np.random.default_rng(seed).normal(size=N)
+
+compare(LEAN, distrib)
\ No newline at end of file
diff --git a/scripts/check_uniform.py b/scripts/check_uniform.py
new file mode 100755
index 0000000..e979b1c
--- /dev/null
+++ b/scripts/check_uniform.py
@@ -0,0 +1,17 @@
+import numpy as np
+from common import compare
+
+print("Checking uniform distribution...")
+
+LEAN = """import RandomDo
+import Batteries.Data.Float.Basic
+def main (args : List String) : IO Unit := do
+  for s in args do
+    IO.FS.withFile (System.FilePath.mk s!"@DIR@/pcg64-{s}.txt") .write fun h ↦
+      IO.runRandPCGWith s.toNat! do
+        for _ in List.range @N@ do h.putStrLn (← NumLean.uniform (-1000000) 1000000).toStringFull
+"""
+
+distrib = lambda seed, N: np.random.default_rng(seed).uniform(-1000000, 1000000, size=N)
+
+compare(LEAN, distrib)
\ No newline at end of file
diff --git a/scripts/common.py b/scripts/common.py
new file mode 100644
index 0000000..5e90412
--- /dev/null
+++ b/scripts/common.py
@@ -0,0 +1,30 @@
+import os, shutil, subprocess, sys
+from decimal import Decimal
+import numpy as np
+
+N, SEEDS, DIR = 1_000_000, np.random.randint(1_000_000_000, size=5), "test_data"
+
+def check_distrib(path, xs):
+    with open(path) as f:
+        for i, (line, x) in enumerate(zip(f, xs)):
+            if line.rstrip("\n") != format(Decimal(float(x)), "f"):
+                return i
+    return None
+
+def compare(lean_code, distrib): 
+  shutil.rmtree(DIR, ignore_errors=True)
+  os.makedirs(DIR)
+  open(f"{DIR}/Dump.lean", "w").write(lean_code.replace("@DIR@", DIR).replace("@N@", str(N)))
+  subprocess.run(["lake", "exe", "dump", *map(str, SEEDS)], check=True)
+
+  ok = True
+  for seed in SEEDS:
+      xs = distrib(seed, N)
+      bad_line = check_distrib(f"{DIR}/pcg64-{seed}.txt", xs)
+      if bad_line is None:
+          print(f"seed {seed} {N} identical draws")
+      else:
+          print(f"seed {seed} DIVERGENCE at line {bad_line}")
+      ok = ok and bad_line is None
+
+  sys.exit(0 if ok else 1)