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..1b94f05 100644 --- a/RandomDo.lean +++ b/RandomDo.lean @@ -6,6 +6,11 @@ public import RandomDo.Monad.ForInInstances public import RandomDo.Monad.Instances public import RandomDo.Monad.MeasurableSpace public import RandomDo.Monad.Notation +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.Deriving public import RandomDo.Tactic.Elab public import RandomDo.Tactic.ForInStep diff --git a/RandomDo/NumLean/Distributions.lean b/RandomDo/NumLean/Distributions.lean new file mode 100644 index 0000000..73640d4 --- /dev/null +++ b/RandomDo/NumLean/Distributions.lean @@ -0,0 +1,138 @@ +/- +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.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 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) + +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. +* The reference implementation: +* 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. + +-/ + +@[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/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_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)