diff --git a/libmath/desk/lib/fixed.hoon b/libmath/desk/lib/fixed.hoon index dd313f2..c9949e6 100644 --- a/libmath/desk/lib/fixed.hoon +++ b/libmath/desk/lib/fixed.hoon @@ -189,6 +189,23 @@ (~(sun rs %n) (abs:si s)) (~(mul rs %n) .-1 (~(sun rs %n) (abs:si s))) (~(div rs %n) fs (~(sun rs %n) (bex b.prec))) +:: +from-rd: [@rd prec] -> @ +:: +:: Quantize a double-precision IEEE float to fixed-point, rounding the +:: scaled value f * 2^b to the nearest integer (ties to even). +:: Out-of-range results wrap (s-to-twoc). Mirrors +from-rs at double +:: precision -- added for /lib/i754rand's distribution-to-fixed-point +:: composition (librand's rand-spec.md section 12.1: "sample at @rd, +:: quantize round-to-nearest"). +:: Examples +:: > (from-rd .~1.5 [8 8]) +:: 0x180 +:: Source +++ from-rd + |= [f=@rd =prec] + ^- @ + =/ i=@s (need (~(toi rd %n) (~(mul rd %n) f (~(sun rd %n) (bex b.prec))))) + (s-to-twoc:(ng prec) i) :: +neg: [@ prec] -> @ :: :: Two's-complement negation of a fixed-point number at its own width. diff --git a/libmath/desk/lib/math.hoon b/libmath/desk/lib/math.hoon index 11eedbe..703658f 100644 --- a/libmath/desk/lib/math.hoon +++ b/libmath/desk/lib/math.hoon @@ -70,9 +70,9 @@ :: Returns the value 1/sqrt(2) (OEIS A010503). :: Examples :: > invsqt2 - :: .70710677 + :: .0.70710677 :: Source - ++ invsqt2 .70710677 + ++ invsqt2 .0.70710677 :: +log2: @rs :: :: Returns the value log(2) (OEIS A002162). diff --git a/libmath/desk/tests/lib/fixed.hoon b/libmath/desk/tests/lib/fixed.hoon index 95e2a3d..1f46d9f 100644 --- a/libmath/desk/tests/lib/fixed.hoon +++ b/libmath/desk/tests/lib/fixed.hoon @@ -106,4 +106,13 @@ %+ expect-eq !>(`@rs`.-1.5) !>((to-rs:fixed n1-5 q88)) :: -1.5 %+ expect-eq !>(`@`0x180) !>((from-rs:fixed (to-rs:fixed 0x180 q88) q88)) :: round-trip == +:: +:: +from-rd: mirrors +from-rs at double precision (added for librand's +:: /lib/i754rand-to-fixed-point composition, rand-spec.md section 12.1). +:: No +to-rd counterpart exists (not asked for), so no round-trip check. +++ test-rd-bridge ^- tang + ;: weld + %+ expect-eq !>(`@`0x180) !>((from-rd:fixed .~1.5 q88)) :: 1.5 + %+ expect-eq !>(`@`0x1.fdc0) !>((from-rd:fixed .~-2.25 q88)) :: -2.25 + == -- diff --git a/libmath/desk/tests/lib/math-constants.hoon b/libmath/desk/tests/lib/math-constants.hoon new file mode 100644 index 0000000..d698bd3 --- /dev/null +++ b/libmath/desk/tests/lib/math-constants.hoon @@ -0,0 +1,56 @@ +:::: /tests/lib/math-constants -- tau/pi/phi/sqt2/invsqt2 at all four +:::: precisions (@rs/@rd/@rh/@rq). +:: +:: Regression for a real bug caught while building librand's /lib/complexrand +:: adapter: /lib/math's @rs +invsqt2 was written `.70710677` (missing the +:: leading `.0.` before the fractional digits), which Hoon parses as the +:: INTEGER 70,710,677.0, not 0.70710677 -- every other precision (rd/rh/rq) +:: had the correct `int.frac` form, so this only affected @rs. Nothing +:: exercised +invsqt2:rs:math until then. Expected bits below are ship- +:: verified (dojo `@ux` casts), not hand-derived. +:: +/+ *test, math +|% +++ test-tau ^- tang + ;: weld + %+ expect-eq !>(`@`0x40c9.0fdb) !>(tau:rs:math) + %+ expect-eq !>(`@`0x4019.21fb.5444.2d18) !>(tau:rd:math) + %+ expect-eq !>(`@`0x4648) !>(tau:rh:math) + %+ expect-eq !>(`@`0x4001.921f.b544.42d1.8469.898c.c517.01b8) !>(tau:rq:math) + == +:: +++ test-pi ^- tang + ;: weld + %+ expect-eq !>(`@`0x4049.0fdb) !>(pi:rs:math) + %+ expect-eq !>(`@`0x4009.21fb.5444.2d18) !>(pi:rd:math) + %+ expect-eq !>(`@`0x4248) !>(pi:rh:math) + %+ expect-eq !>(`@`0x4000.921f.b544.42d1.8469.898c.c517.01b8) !>(pi:rq:math) + == +:: +++ test-phi ^- tang + ;: weld + %+ expect-eq !>(`@`0x3fcf.1bbd) !>(phi:rs:math) + %+ expect-eq !>(`@`0x3ff9.e377.9b97.f4a8) !>(phi:rd:math) + %+ expect-eq !>(`@`0x3e79) !>(phi:rh:math) + %+ expect-eq !>(`@`0x3fff.9e37.79b9.7f4a.7c15.f39c.c060.5cee) !>(phi:rq:math) + == +:: +++ test-sqt2 ^- tang + ;: weld + %+ expect-eq !>(`@`0x3fb5.04f3) !>(sqt2:rs:math) + %+ expect-eq !>(`@`0x3ff6.a09e.667f.3bcd) !>(sqt2:rd:math) + %+ expect-eq !>(`@`0x3da8) !>(sqt2:rh:math) + %+ expect-eq !>(`@`0x3fff.6a09.e667.f3bc.c908.b2fb.1366.ea95) !>(sqt2:rq:math) + == +:: +:: +invsqt2:rs is the arm the bug above lived in -- checked against its +:: IEEE-754 value directly (0.7071067690849304 as a float32), not just +:: bit-equality with a value that could itself be wrong. +++ test-invsqt2 ^- tang + ;: weld + %+ expect-eq !>(`@`0x3f35.04f3) !>(invsqt2:rs:math) + %+ expect-eq !>(`@`0x3fe6.a09e.667f.3bcd) !>(invsqt2:rd:math) + %+ expect-eq !>(`@`0x39a8) !>(invsqt2:rh:math) + %+ expect-eq !>(`@`0x3ffe.6a09.e667.f3bc.c908.b2fb.1366.ea95) !>(invsqt2:rq:math) + == +-- diff --git a/librand/NEXT-STEPS.md b/librand/NEXT-STEPS.md new file mode 100644 index 0000000..c9d1b73 --- /dev/null +++ b/librand/NEXT-STEPS.md @@ -0,0 +1,347 @@ +# `/lib/rand` — next steps + +Status as of 2026-07-13: **all nine milestones of `rand-spec.md` section 13 +are done.** This file is a development log (decisions, footguns, bugs +found and fixed) kept in full rather than rewritten, per this project's own +convention of not silently erasing history — see "Deferred to v2" at the +bottom for what's actually left to do. Milestones 1-6 (`rand-spec.md` +section 13) landed first: +`++split-mix`, `++philox`, `++seed`, `+step`, `+fork`, `++gen`, `++uni`, +`++pcg`, `++dist` (both `++rd` and `++rs`), `++sample` -- KAT-verified +against Vigna's reference SplitMix64, the Random123 `kat_vectors` file, an +independent Python IEEE-754 encoder for `++uni`'s float arms, pcg-c's own +seed=42/seq=54 demo convention for `++pcg`, NumPy's `random_poisson_ptrs` +for `++dist`'s PTRS branch, and on-ship-computed regression values plus +moment tests elsewhere. + +## Architecture: split into `/lib/rand` (plumbing) + `/lib/i754rand` (porcelain) + +After milestone 6 landed, the single `rand.hoon` file was split in two, +matching how `/lib/twoc`/`fixed`/`complex`/`unum` are already kept separate +from `/lib/math` in this codebase: + +- **`/sur/rand`** -- `+$phil`/`+$sm64`/`+$pcg64`/`+$rng`, extracted so any + file needing just the type can use a lightweight `/-` import. +- **`/lib/rand`** (the plumbing) -- `++philox`, `++split-mix`, `++pcg`, + `++seed`, `+step`, `+fork`, `++gen`, `++uni` (integer-only: `+bits`/ + `+below`/`+between`), `++sample` (shuffle/permutation/choice/reservoir). + Every planned adapter library (`i754rand`, and the non-float ones below) + imports this file. Tests: `-test %/tests/lib/rand ~` (44 tests). +- **`/lib/i754rand`** (the IEEE-754 porcelain) -- `++uni`'s float arms + (`+rs`/`+rd`/`+rh`/`+rq`/`+rs-oo`/`+rd-oo`), `++alias` (Vose's method, + moved here from `++sample`), and `++dist`. Tests: + `-test %/tests/lib/i754rand ~` (47 tests). + +`++alias` moved out of `++sample` specifically because `+draw` needs an +actual uniform `@rd` float draw to decide accept-vs-redirect -- it isn't +float-free the way the rest of `++sample` is, so it belongs with the +porcelain, not the plumbing. `+categorical` (in `++dist`) references it as +a same-file sibling (`alias-table`/`draw:alias`, unqualified). + +Import mechanics worth remembering if this gets refactored again: `/+ rand` +alone does NOT make the unwrapped `/sur/rand` types resolve as `rng:rand` +externally unless `/lib/rand` itself does `=+ rand` after its own `/- rand` +import (confirmed empirically -- `rng:rand` works from any consuming file +once `rand.hoon` does the unwrap). Adding a second `=+ rand` inside a +*consumer* file (e.g. `i754rand.hoon`) to get `rng` bare, instead of +`rng:rand`, breaks qualified sibling lookups like `bits:uni:rand` in ways +that were not obvious from the error message (`-find.uni`) -- so +`i754rand.hoon` does NOT unwrap; every bare `rng` in its arm signatures is +written `rng:rand` explicitly instead. + +**Spec correction found during milestone 4**: `rand-spec.md` section 3.3 +originally said PCG64 outputs the current state's permutation and THEN +advances. That's backwards -- pcg-c's actual reference +(`pcg_setseq_128_xsl_rr_64_random_r`, cross-checked against NumPy's vendored +copy) advances the state FIRST and outputs from the new state. Both the +spec and `++pcg` now follow the verified reference order; see the struck- +through note left in the spec for the record. + +**Two real algorithm bugs caught during milestone 6's on-ship testing** +(both are exactly why this library tests against ship-computed values, +not hand-derivation): `++alias`'s `+build-loop` had the small/large +worklist assignment reversed (an index whose weight said it belonged on +one worklist was pushed onto the other); `++dist`'s `+binomial` compared +the uniform draw against the individual pmf term instead of the +accumulated CDF, which either returned 0 too often or (for larger n) ran +the recurrence past n and crashed with subtract-underflow. Both are now +documented in the arms' own doccords so a future reader sees why the code +looks the way it does. + +Note on `++dist`: `+chi2`/`+student-t` aren't in milestone 5's literal arm +list but are one-line compositions of `+gamma`/`+normal` with no new +machinery, so they were added alongside rather than left in limbo until an +unscheduled milestone. `+categorical` has no `++rs` mirror: alias-table's +`.prob` is fixed at `@rd`, so it's precision-invariant. + +Remaining milestones, in dependency order (see `rand-spec.md` section 13 for +full detail): + +7. Non-float adapters, as four SEPARATE libraries (not nested in `/lib/rand` + or `/lib/i754rand`), each importing `/lib/rand` for the engine/uni/sample + primitives: + - DONE: `/lib/twocrand` -- `+twoc-full`, `+twoc-between`. Needs only + `/lib/rand` + `/lib/twoc`; no floats at all. + - DONE: `/lib/fixedrand` -- `+fixed`, `+fixed-unit`, `+fixed-between` + (delegates to `twocrand`'s `+twoc-between`). Needs `/lib/rand` + + `twocrand` + `/lib/fixed` (needed a `+from-rd` added, mirroring its + existing `+from-rs` -- done). Footgun hit and documented: `+fixed` + (the arm, matching rand-spec.md's public API) collides with `fixed` + (the imported library face), since a same-named battery arm shadows + an imported face throughout its own core -- same class of bug + `/lib/fixed` itself hit renaming `+twoc` to `+neg`. Fixed here by + aliasing the import (`/+ fx=fixed`) instead of renaming the arm, + since the arm name is spec-mandated and the import name isn't. + - DONE: `/lib/complexrand` -- `+cuniform`, `+normal-parts`, `+cnormal`, + `+on-circle`, `+in-disk`. Every arm here draws floats directly, so + this one DOES depend on `/lib/i754rand` (for `+rd:uni`/`+rd-oo:uni`), + plus `/lib/complex` and `/lib/math` (cos/sin). One arm set per + component-width door (`+cd` first, `+cs` mirror second, matching + `/lib/complex`'s own ship order). Caught a real, previously- + uncovered bug while building `+cnormal:cs`: `/lib/math`'s `@rs` + `+invsqt2` was written `.70710677` (missing the leading `.0.`), + which Hoon parses as the integer 70,710,677.0, not 0.70710677 -- + every other precision (`@rd`/`@rh`/`@rq`) had the correct form. + Fixed, with a new regression suite (`tests/lib/math-constants.hoon` + in `libmath`) covering `tau`/`pi`/`phi`/`sqt2`/`invsqt2` at all four + precisions, since nothing previously exercised any of them. + - DONE: `/lib/unumrand` -- `+posit-lattice`, `+posit-unit`. Needs only + `/lib/rand` + `/lib/unum`; no floats (`+posit-unit`'s construction is + explicitly float-free per rand-spec.md section 12.4). Distributions + ("sample at @rd, convert") need NO new /lib/unum plumbing, unlike + fixedrand's `+from-rd` -- `/lib/unum` already ships `+from-rh/rs/rd/rq` + at every width door. + - Range subtlety (worked through with the user before implementation): + `+posit-lattice` (uniform over bit patterns, minus NaR) and + `+posit-unit` (uniform over VALUES on [0,1), exact) are fundamentally + different because posits are tapered -- consecutive bit patterns are + NOT evenly spaced in value. `+posit-lattice` ships at all FIVE width + doors (rpb/rph/rps/rpd/rpq); it has no bit-count subtlety (simple + NaR-rejection, exact at any width). + - `+posit-unit`'s bit-count k is the dangerous part: draw k raw bits u, + encode the dyadic u*2^-k via /lib/unum's existing +bit (RNE, + saturating -- confirmed no rounding-mode parameter needed). For this + to be EXACTLY the round-to-nearest image of continuous uniform (not + approximately), k must exceed the finest rounding-cell width anywhere + in [0,1) -- which, because of the taper, is NEAR ZERO and shrinks fast + with width. Hand-traced /lib/unum's own `+sea` decode of posit8's + pattern `1` (minpos) and confirmed minpos = 2^-4(n-2) exactly (2^-24 + for n=8), so the rounding boundary nearest zero sits at 2^-(4(n-2)+1). + The spec's k values (32/64/128 for posit8/16/32) are exactly k=4n, + which gives a CONSTANT 7-bit safety margin at any width + (4n - (4(n-2)+1) = 7 always) -- naive choices like k=n or k=n+8 would + be silently wrong (they'd truncate resolution near zero and bias the + distribution with no obvious symptom). u=0 must produce `[%z ~]` + explicitly per spec (not a degenerate `[%p ...]` with a=0). Scoped to + `[0,1)` only -- a general bounded-range posit uniform is out of scope. + - Decided with the user: `+posit-unit` stays scoped to posit8/16/32 + (matching /lib/unum's existing "unverified until oracle sweep extends" + caveat for posit64/128 -- even though the k=4n margin argument would + work mathematically at any width, since +bit's encoding logic isn't a + convergence-based transcendental, staying scoped avoids a false sense + of rigor where the rest of the library hasn't independently verified + those widths either). + - Decided with the user on oracle rigor: chi-square posit8 (256 patterns, + exactly enumerable) AND posit16 (65,536 patterns, still cheap to + enumerate exactly) against the mpmath oracle; posit32 (2^32 patterns, + exact enumeration infeasible) gets ship-verified value-regression + tests only, relying on the same proven k=4n margin argument rather + than independent re-verification. + - Plan revised during implementation: originally planned as a real Hoon + `-test` regression (hardcoded expected statistic/threshold). Dropped + in favor of `librand/tools/posit_unit_check.py` staying a standalone + oracle (no mpmath needed -- everything here is already exactly dyadic, + so plain `Fraction` arithmetic is exact throughout): embedding a + 256- or 65536-row expected-probability table as Hoon literals isn't + practical, and unlike the moment tests (which check a couple of + summary statistics), this needs the FULL per-pattern distribution to + mean anything. `posit_unit_check.py table ` prints the exact + table; `posit_unit_check.py chi2 --bins N` + chi-squares ship-drawn raw patterns (one per line) against it, + quantile-binned so N stays modest. The oracle script's OWN + correctness is validated in-repo (a simulated correct-distribution + sample passes with a sane p-value; a deliberately-wrong one is + rejected with p~0 -- proving the test has real power, not just + rubber-stamping). Actually run once at posit8 with N=100,000 ship- + drawn draws (fixed seed 12345), tallied on-ship into a histogram (at + most 65 possible patterns for posit8, so the histogram itself is + small/printable) via a throwaway `.hoon` file (NOT a raw dojo one- + liner -- a long single-line multi-`=/`/`|-` dojo command got silently + truncated in transmission mid-session and left the dojo edit buffer + stuck with an unclosed expression; recovered with Ctrl-E then Ctrl-U, + not Ctrl-U alone, since the cursor was stuck at the line's start; + lesson: bulk/long computations belong in a deployed `.hoon` file run + via `-build-file`, not a giant single dojo line). Result: 61/65 + patterns hit (the 4 misses are the lowest-probability patterns near + zero, expected counts well under 1 at this N), binned into 20 + quantile bins (some of the highest-probability patterns near 1.0 + individually exceed 1/32 of the mass, so `--bins 32` yields fewer + than 32 actual bins -- expected, see the script's own comment), + chi-square = 26.86 on 19 dof, p = 0.108 -- consistent with the exact + expected distribution at any standard significance level. + - DONE (added post-milestone-9, for a paper on this construction -- + see `~/Documents/journal-posit-prng/`): exhaustive input-space + verification at posit8. `librand/tools/posit8_exhaustive.c` is a + from-scratch C re-derivation of the encode logic (not a + transliteration of `posit_unit_check.py`'s own `encode()` -- + deliberately independent, so a translation bug in one isn't + invisible to the other), enumerating all 2^32 possible numerators u + and tallying per-pattern counts exactly. Cross-validated against the + Python `encode()` on 2,000,000 random u plus the domain edges before + being trusted (zero mismatches), then run to completion (~6s at + `-O3`) and compared against `expected_probabilities('posit8')`'s + numerators EXACTLY -- 65/65 patterns match, all 2^32 draws + accounted for. This eliminates sampling error entirely at posit8; + the chi-square result above is now a corroborating footnote, not + the load-bearing evidence. Run via + `python3 librand/tools/posit_unit_check.py exhaustive8` (compiles + and runs the C harness, then does the comparison, all in one + command). + - DONE (also added post-milestone-9, for the paper): posit16 chi-square + run at N=1,000,000. posit16 has 16,385 reachable patterns -- far too + many to safely dump a raw per-pattern histogram out of a live dojo + session (confirmed the hard way: the terminal scrollback silently + truncates the printed map well before the full histogram appears, + even after raising tmux's `history-limit`; a first attempt at fixing + this by embedding a 16,385-entry pattern-to-bin lookup table as one + big Hoon literal crashed the ship outright, twice). Fix: have the + SHIP do the binning itself into the SAME 64 quantile bins the chi- + square already uses, via a small (64-entry) list of bin-boundary + PATTERNS rather than a full lookup table -- posit16's reachable + patterns are nonneg-valued, so they sort identically by raw integer + and by decoded value, meaning each quantile bin is already a + contiguous range of raw pattern integers and only needs its lower + boundary recorded. `librand/tools/posit_unit_check.py gen-binned + posit16 out.hoon --bins 64 --count 1.000.000` generates the + deployable harness; `chi2-binned` consumes the captured ship output. + Two Hoon lessons paid for getting the harness itself right (both + baked into `gen_binned_hoon()`'s own header comment so they don't + recur): `~[a b c ...]` list literals infer a FIXED-SHAPE tuple type, + not the general recursive `(list @)` mold, so using one as a `|-` + trap's loop variable across recursive `$(...)` calls fails with a + confusing `mint-vain` (the shape genuinely differs each iteration, + shrinking one element per recursion) -- always `^- (list @)` cast + list literals used this way; and an unrelated off-by-one in the + first `bin-of` draft (incrementing the running index BEFORE + checking the boundary it names, not after) produced a wildly wrong + first chi-square (15,700 on 63 dof) that would have been mistaken + for a real RNG bug if not caught by verifying `bin-of` against known + boundary values first. Corrected result: N=1,000,000, χ²=71.78 on 63 + dof, p=0.210 -- consistent with the exact expected distribution. + The bit-exact cross-check against `posit_unit_check.py`'s `encode()` + at posit16/32 (done for several seeds, in `+test-posit-unit`) + remains a stronger per-draw guarantee than either statistical test, + since it's the same deterministic code path at every width. + - Decided: the "sample at @rd, quantize/convert" pattern mentioned in + rand-spec.md 12.1/12.4 (fixed-point and posit *distributions*) is + NOT shipped as dedicated per-distribution wrapper arms (that would be + ~20 nearly-identical one-liners across fixedrand/unumrand doing + nothing beyond composing two already-existing arms). Instead: document + the one-line composition (`i754rand`'s dist arm -> the target type's + own `+from-rd`) in each adapter's own doccord, and add ONE test per + adapter proving the composition actually works (this also validates + `+from-rd` once it's added to `/lib/fixed`). Revisit if real callers + end up wanting the same composition in more than one or two places. + - Split into two passes: twocrand + fixedrand + complexrand first + (mechanical, same ship-verification rhythm as prior milestones -- + ALL THREE ARE DONE); then unumrand alone, since `++posit-unit`'s + value-uniformity claim needs a + NEW mpmath oracle script (`librand/tools/posit_unit_check.py`, + alongside `unum_cheb_check.py`) before it can be trusted, unlike + everything shipped so far which either had a real external KAT or was + ship-verified against a hand/Python-traced computation. +8. DONE: Saloon `+rand-ray` (counter-window design) + replay-equality tests. + Lands directly in `/lib/saloon`'s `+sa` core (a new `+| %rand` section), + not a separate librand file -- the spec's own framing ("In /lib/saloon, + a core taking a Lagoon meta and an rng") describes per-arm shape, not a + literal nested sub-core; Saloon's existing convention is one flat `+sa` + door with `+|` section markers, so `+rand-ray` follows that rather than + introducing complexrand/unumrand-style nesting. + - `+fill-uniform` (single non-rejecting draw/element, %phil gets + ctr0+i directly), `+fill-normal`/`+fill-expon` (rejection-based, + %phil gets the ctr0+i*2^32 window with an explicit crash on + exhaustion), `+fill-below` (%uint rays via Lemire, also windowed). + `%i754` only, bloq 5/6 (`@rs`/`@rd`) -- matches rand-spec.md's stated + v1 scope ("posit rays deferred"). + - Blocked on a pre-existing, unrelated bug: `+sa`'s scalar transcendental + dispatch (`+trans-scalar`, `+fadd`/etc) called `/lib/math`'s doors as + `[rnd rtol]` (2-tuple) where they need `[r rtol atol]` (3-tuple) -- + `saloon.hoon` didn't `-build-file` AT ALL before this, on ANY branch, + confirmed by testing the unmodified file directly. Fixed and shipped + as its OWN PR (#77, `sigilante/saloon-scalar-rtol-fix` off `main`), + since it's also needed upstream in `urbit/urbit` independent of + librand. That fix is ALSO carried on this branch (duplicated, not + rebased) so `+rand-ray` has something to build against; reconcile via + rebase once #77 merges to main. + - Two Hoon footguns hit writing `+rand-ray` itself, both the "narrowing + doesn't survive a recursive `$(...)` rebind" class already known in + this codebase: (a) mutating `.r` (`r(ctr.p ...)`) inside a `?=(%phil + -.r)`-narrowed recursive trap lost the narrowing on recursive re- + entry -- fixed by building FRESH `[%phil key0 ctrN]` literals instead + of mutating; (b) the non-%phil ("sequential engine") branch's `.r` + is narrowed to EXCLUDE %phil by the same `?:`, but `+draw`'s return + type is the full `rng` union, so the trap's first entry (narrow) and + recursive re-entries (widened by `=^ v r (draw r)`) disagreed -- + fixed by explicitly widening `` `rng:rand`r `` once before the trap + so every entry matches. + - Every arm bit-exact cross-checked against a DIRECT call to the + underlying `/lib/rand`/`/lib/i754rand` primitive at the expected + counter (not just "doesn't crash") -- see `tests/lib/saloon-rand- + ray.hoon` in the `saloon` desk (not `librand`, since the code lives + in Saloon). +9. DONE: full README rewrite + this file's final pass (this section). + +## Deferred to v2 + +Six items, all deliberately out of v1 scope per `rand-spec.md`. None of +these are bugs or gaps in what's shipped — they're follow-on work a future +session can pick up independently, in no particular order: + +- **Ziggurat.** `+normal`/`+expon` use Marsaglia polar / inversion, which + are simple and correctly rounded but do a `+sqt`/`+log` call per + accepted sample (Marsaglia's rejection rate is ~21.5%, so ~1.27 draws + and one transcendental call per deviate on average). Ziggurat avoids + the transcendental call entirely in the common case via precomputed + layer tables, at the cost of real implementation complexity (building + and validating the tables) for a speedup that mostly matters once these + arms are jetted. Not worth it before then. +- **BTPE for `+binomial`.** The current inversion-by-CDF-accumulation + method crashes above `n*min(p,1-p) >= 30` (the recurrence gets slow and + numerically risky past that point). Kachitvichyanukul & Schmeiser's + BTPE algorithm (Binomial, Triangle, Parallelogram, Exponential regions) + handles the large-n regime in O(1) expected time, but it's a + substantially more involved algorithm than anything else in `++dist` — + a real follow-on project, not a quick add. +- **Buffered Philox.** `+step`'s `%phil` branch keeps only the low 64 + bits of each 128-bit Philox4x32-10 block and discards the high 64 bits + (`/lib/rand`'s own `+step` doc comment: "wasting them is acceptable at + v1... a buffered variant that reuses both halves is a NEXT-STEPS + item"). A buffered variant would cache the unused half and serve it on + the *next* `+step` call instead of computing a fresh block, roughly + doubling sequential throughput with zero extra Philox rounds. Doesn't + change any output bit-for-bit; purely a performance follow-up, and one + that interacts with jetting (the cache would need to live in the `rng` + noun itself, changing `+$phil`'s shape) so it's more natural to do + alongside milestone 11's jetting work than before it. +- **Posit rays.** Saloon's `+rand-ray` (milestone 8) only fills `%i754` + rays. Extending `+fill-uniform`/`+fill-normal`/etc. to `%unum`-kind rays + (posit-valued Lagoon arrays), using `/lib/unumrand`'s `+posit-unit`/ + `+posit-lattice` as the per-element generator, is mechanically similar + to the existing `%i754` path but was out of scope for the milestone + that shipped it. +- **Quire Monte Carlo.** Noted when `/lib/unumrand` shipped: `/lib/unum`'s + quire (`+fdp`, Type III unums' exact fixed-point accumulator) makes + sample *sums* singly-rounded — a capability hardware IEEE floats simply + don't have (every float addition in an accumulation loop rounds; a + quire-based reduction rounds exactly once, at the very end). This is a + genuinely novel capability worth designing a Monte-Carlo-reduction API + around, but it's a design project in its own right, not a small + addition to an existing arm. +- **`@rh`/`@rq` distributions.** `/lib/i754rand`'s `++dist` only has + `++rd` (reference, double)/`++rs` (single) mirrors, matching + rand-spec.md's own scoping ("each arm exists at @rd (reference) and + @rs"). Half (`@rh`) and quad (`@rq`) precision distributions were never + in v1 scope; adding them is a mechanical re-instantiation of the same + algorithms (matching how `++rs` itself was built as a mirror of + `++rd`), gated only by whether a real caller needs sub-single or + above-double precision sampling. diff --git a/librand/README.md b/librand/README.md new file mode 100644 index 0000000..aa47858 --- /dev/null +++ b/librand/README.md @@ -0,0 +1,216 @@ +# `/lib/rand`: deterministic pseudo-random number generation for Urbit + +Reproducible random number generation for numerical work — Monte Carlo, +ML initialization, shuffles, simulation. **Not cryptographic**: entropy +acquisition is Arvo's job (`eny`), cryptographic randomness is Zuse's job. +Given the same seed, every arm here produces the same output on any +ship, any time, forever — that's the entire point. + +Full design is in `rand-spec.md` (repo root). All nine of its milestones +are done: three engines, integer and float uniform deviates, a full +distribution suite, sampling utilities, four non-float output adapters, +and a Lagoon array-filling layer in Saloon. Hoon reference implementation +first; jets are a follow-up (`rand-spec.md` section 11). + +## Layout: seven libraries, not one + +The whole surface is split by concern, mirroring how `/lib/twoc`/`fixed`/ +`complex`/`unum` are already kept separate from `/lib/math` in this +codebase: + +- **`/sur/rand`** — just the types (`+$rng`, `+$phil`, `+$sm64`, + `+$pcg64`), so anything that only needs the type gets it via a + lightweight `/-` import instead of pulling in the whole library. +- **`/lib/rand`** — the **plumbing**: the three engines, seeding, + engine-dispatched draw/fork, integer-only uniform deviates, generic + sampling. Every other library here imports this one. +- **`/lib/i754rand`** — the **porcelain**: IEEE-754 float generation, + Vose's alias method, and the full distribution suite. This is + `/lib/rand`'s own "math.hoon" — kept separate so the plumbing and the + non-float adapters below never need to pull in `/lib/math`. +- **`/lib/twocrand`**, **`/lib/fixedrand`**, **`/lib/complexrand`**, + **`/lib/unumrand`** — non-float output adapters (two's-complement + integers, fixed-point, complex numbers, posits). Siblings of + `i754rand`, not dependents of it: only `complexrand` needs floats for + its own uniform-generation arms (every arm draws floats directly); the + others only optionally reach for `i754rand` to compose a "sample at + `@rd`, then quantize/convert" distribution (documented as a pattern, + not shipped as dedicated wrapper arms — see `NEXT-STEPS.md`). +- **`/lib/saloon`** (a different desk) — the array-filling layer, + `+rand-ray`, is not here; it lives directly in Saloon's `+sa` core (see + below). + +## What's here + +### `/lib/rand` (plumbing) + +- **Engines**: `++philox` (Philox4x32-10, counter-based, the primary + engine — its whole design point is that element `i` of an array only + needs counter `ctr0+i`, so a jet can fill array elements in parallel + and land bit-identical results regardless of thread scheduling), + `++split-mix` (SplitMix64, a small sequential generator and the + seed-mixing primitive every engine's seeding/`+fork` path depends on), + `++pcg` (PCG64 XSL-RR, with `O(log n)` `+jump`). All three KAT-checked + against independent reference implementations (Random123's + `kat_vectors`, Vigna's public-domain SplitMix64 C, pcg-c's own + seed=42/seq=54 demo convention). +- **`++seed`** — `+from-atom`, `+from-eny`, `+fold-wide`, `+mix`. +- **`+step`/`+fork`** — the generic engine-dispatched draw and the + path-sensitive key derivation (JAX-style: fork by path instead of + threading one sequential state, so a draw inserted in one branch of a + computation never perturbs a sibling branch) across all three engine + shapes. `++gen` is a thin door facade over both. +- **`++uni`** (integer only) — `+bits`, `+below` (Lemire's unbiased + method — the primitive every bounded-integer draw in this codebase + uses; modulo bias is never an option here), `+between`. +- **`++sample`** — `+shuffle`/`+permutation`/`+choice`/`+choices`/ + `+sample-n`/`+reservoir` (Algorithm R). + +### `/lib/i754rand` (porcelain) + +- **`++uni`** — the four float auras (`+rs`/`+rd`/`+rh`/`+rq`) plus + open-open variants (`+rs-oo`/`+rd-oo`), exact bit constructions checked + against an independent Python IEEE-754 encoder. +- **`++alias`** — Vose's alias method (`+build`/`+draw`), moved here + from `++sample` because `+draw` needs an actual float draw to decide + accept-vs-redirect. +- **`++dist`**, nested `++rd` (reference, double precision) / `++rs` + (single precision, a mechanical re-instantiation of the same + algorithms) — `+normal` (Marsaglia polar), `+normal-mv`, `+expon` + (inversion), `+gamma` (Marsaglia-Tsang), `+beta`, `+chi2`, + `+student-t`, `+bernoulli`, `+geometric`, `+categorical` (`@rd` + only — the alias table's probabilities are fixed at that precision), + `+poisson` (Knuth for `lambda<10`, Hörmann's PTRS above — verified + against NumPy's `random_poisson_ptrs`), `+binomial` (inversion by CDF + accumulation; crashes above `n*min(p,1-p) >= 30`, where BTPE would be + needed — see `NEXT-STEPS.md`), `+dirichlet`. Moment-tested (mean/ + variance regression at 50k draws, fixed seed) for normal/expon/gamma. + `@rh`/`@rq` distributions are out of scope for v1 (see `NEXT-STEPS.md`). + +### The four non-float adapters (`rand-spec.md` section 12) + +- **`/lib/twocrand`** — `+twoc-full` (raw-bit passthrough) and + `+twoc-between` (Lemire-unbiased inclusive range in two's-complement + order, via `/lib/twoc`'s width-keyed `+twid` door). Depends only on + `/lib/rand` + `/lib/twoc`. +- **`/lib/fixedrand`** — `+fixed`, `+fixed-unit` (both raw-bit + passthroughs — a fixed-point lattice is uniform by construction), + `+fixed-between` (delegates to `twocrand`). Needed a new `+from-rd` + added to `/lib/fixed`, mirroring its existing `+from-rs`. +- **`/lib/complexrand`** — `+cuniform`, `+normal-parts`, `+cnormal` + (deliberately distinct from `+normal-parts`: "complex Gaussian" means + different things to signal-processing and statistics users), + `+on-circle`, `+in-disk`. One arm set per component-width door (`+cd` + reference precision, `+cs` mirror), matching `/lib/complex`'s own ship + order. Every arm draws floats directly, so this is the one adapter + that depends on `/lib/i754rand`, `/lib/math`, and `/lib/complex`. +- **`/lib/unumrand`** — two *inequivalent* uniform semantics over posits, + named so they can't be confused (posits are tapered: consecutive bit + patterns are not evenly spaced in value). `+posit-lattice` is uniform + over raw bit patterns excluding NaR (a fuzzing/property-testing + primitive — the induced value distribution is only approximately + log-uniform); it ships at all five width doors (`rpb`/`rph`/`rps`/ + `rpd`/`rpq`), with no bit-count subtlety at any width. `+posit-unit` is + uniform over *values* on `[0,1)`, exact; it's scoped to posit8/16/32 + only (see `NEXT-STEPS.md` for the `k=4n` bit-count derivation behind + its exactness claim, and why posit64/128 aren't included). Needs no new + `/lib/unum` plumbing — `+from-rh/rs/rd/rq` already exist at every width + door. Verified against a new exact-rational oracle + (`tools/posit_unit_check.py`) via chi-square at posit8 (100,000 + ship-drawn draws, p=0.108) and posit16 (1,000,000 draws quantile-binned + on-ship into 64 groups — posit16 has 16,385 reachable patterns, too many + to safely extract a raw histogram from dojo — χ²=71.78 on 63 dof, + p=0.210), an exhaustive input-space check at posit8 (all 2^32 possible + draws, via an independent C harness, `tools/posit8_exhaustive.c` — + every one of the 65 reachable patterns' counts matches the oracle's + numerators exactly, eliminating sampling error entirely), plus + bit-exact cross-checks at every in-scope width. + +None of the four ship dedicated per-distribution wrapper arms for the +"sample at `@rd`, quantize/convert" pattern (fixed-point and posit +distributions) — that would be a dozen-plus nearly-identical one-liners +composing two already-existing arms. Instead each adapter documents the +one-line composition and has one test proving it works end to end. + +### `+rand-ray` (in `/lib/saloon`, `rand-spec.md` section 8) + +Not part of this library's own file tree — it's a new `+| %rand` section +in Saloon's `+sa` core, filling a Lagoon `$ray` element-by-element in +row-major (C) order: `+fill-uniform`, `+fill-normal`, `+fill-expon` +(`%i754` rays, bloq 5/6 — `@rs`/`@rd` only; posit rays are deferred), +`+fill-below` (`%uint` rays via Lemire). Philox elements get per-element +counter treatment — the entire reason Philox is the primary engine: +`+fill-uniform`'s single non-rejecting draw per element assigns element +`i` counter `ctr0+i` directly, so the whole fill decomposes into `n` +independent draws a future jet can parallelize. The rejection-based fills +(`+fill-normal`/`+fill-expon`/`+fill-below`) get a wider per-element +counter *window* instead (`ctr0 + i*2^32`): each element's rejection loop +walks freely within its own window, and the returned rng's counter is +forced to `ctr0 + n*2^32` regardless of how many sub-draws any element +actually used, so the post-state is a pure function of `n` — with an +explicit crash if a single element's rejection loop ever walks past its +own window (astronomically improbable, but checked rather than silently +overflowing into the next element's window). Non-Philox engines have no +comparable jump primitive and aren't the ones a jet would parallelize +anyway, so they just thread sequentially. + +Every arm is bit-exact cross-checked against a direct call to the +underlying `/lib/rand`/`/lib/i754rand` primitive at the expected counter +— see `saloon/desk/tests/lib/saloon-rand-ray.hoon` (in the `saloon` desk, +not here). + +## Testing philosophy + +Every arm is verified against ship-computed output before being trusted — +this codebase does not hand-derive expected values for anything involving +floating-point transcendentals, bit-level float/posit construction, or +RNG state transitions, because a hand or Python derivation isn't +guaranteed to match this codebase's own kernels bit-for-bit. Known-answer +tests come from independent reference implementations wherever one +exists (Random123, Vigna's SplitMix64, pcg-c, NumPy's PTRS, an exact- +rational posit oracle); everything else is a ship-computed regression +value pinned into the test file. Two real algorithm bugs (a reversed +worklist in the alias method, a wrong-accumulator comparison in +`+binomial`) and two real pre-existing bugs in dependencies (a malformed +`@rs` float literal in `/lib/math`'s `+invsqt2`; a sample-shape mismatch +in Saloon's scalar transcendental dispatch) were caught exactly this way. + +``` +-test %/tests/lib/rand ~ +-test %/tests/lib/i754rand ~ +-test %/tests/lib/twocrand ~ +-test %/tests/lib/fixedrand ~ +-test %/tests/lib/complexrand ~ +-test %/tests/lib/unumrand ~ +-test %/tests/lib/saloon-rand-ray ~ :: in the saloon desk +``` + +## Layout + +``` +librand/ + README.md + NEXT-STEPS.md + rand-spec.md :: (repo root) the governing design spec + desk/sur/rand.hoon :: +$rng, +$phil, +$sm64, +$pcg64 + desk/lib/rand.hoon :: the plumbing + desk/lib/i754rand.hoon :: the IEEE-754 porcelain + desk/lib/twocrand.hoon :: two's-complement integer adapter + desk/lib/fixedrand.hoon :: fixed-point adapter + desk/lib/complexrand.hoon :: complex-number adapter + desk/lib/unumrand.hoon :: posit adapter + desk/tests/lib/rand.hoon :: -test %/tests/lib/rand ~ + desk/tests/lib/i754rand.hoon :: -test %/tests/lib/i754rand ~ + desk/tests/lib/twocrand.hoon :: -test %/tests/lib/twocrand ~ + desk/tests/lib/fixedrand.hoon :: -test %/tests/lib/fixedrand ~ + desk/tests/lib/complexrand.hoon :: -test %/tests/lib/complexrand ~ + desk/tests/lib/unumrand.hoon :: -test %/tests/lib/unumrand ~ + tools/posit_unit_check.py :: +posit-unit's exact-rational oracle + +saloon/desk/lib/saloon.hoon :: +rand-ray lives in +sa's +| %rand section +saloon/desk/tests/lib/saloon-rand-ray.hoon +``` + +See `NEXT-STEPS.md` for what's deliberately deferred past v1 (ziggurat, +BTPE, buffered Philox, posit rays, quire Monte Carlo, `@rh`/`@rq` +distributions) and the full development history/decision log. diff --git a/librand/desk/lib/complexrand.hoon b/librand/desk/lib/complexrand.hoon new file mode 100644 index 0000000..a662970 --- /dev/null +++ b/librand/desk/lib/complexrand.hoon @@ -0,0 +1,188 @@ +/+ rand, i754rand, complex, math +:::: /lib/complexrand -- complex-number output adapter (rand-spec.md +:::: section 12.3) +:: +:: One arm set per component-width door: `++cd` (double, @cd, the +:: reference precision) first, `++cs` (single, @cs, a mechanical +:: re-instantiation of the same arms) second -- matching /lib/i754rand's +:: own `++rd`-then-`++rs` convention, NOT /lib/complex's file order +:: (which actually defines `++cs` before `++cd`). +:: Every arm here draws floats directly (there is no raw-bit passthrough +:: the way twocrand/fixedrand have), so -- unlike those two -- this +:: adapter DOES depend on /lib/i754rand (for +rd:uni/+rd-oo:uni and +:: ++dist's +normal), plus /lib/math (for +cos/+sin/+tau/+invsqt2) and +:: /lib/complex (for +pak/+re/+im, to assemble/inspect the packed atom). +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:: +cd: complex-double (@cd) adapters. +:: +:: Uses /lib/math's rd door and /lib/i754rand's @rd generators +:: throughout. +:: Source +++ cd + |% + ++ m ~(. rd:math [%n .~1e-13 .~0]) + :: +cuniform: rng -> [@cd rng] + :: + :: Uniform over the unit square [0,1) x [0,1): two independent [0,1) + :: component draws, packed. + :: Examples + :: > `@ux`out:(cuniform:cd:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3fe8.9e6a.a1b9.65f4.3f95.072f.63b9.b5e0 + :: Source + ++ cuniform + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (rd:uni:i754rand r) + =^ v r (rd:uni:i754rand r) + [(~(pak cd:complex %n) u v) r] + :: +normal-parts: rng -> [@cd rng] + :: + :: Independent N(0,1) real and imaginary components (two draws from + :: /lib/i754rand's +normal, one per component). + :: Examples + :: > `@ux`out:(normal-parts:cd:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3ff2.e308.2bd4.88f8.bfee.55f7.4af6.d8b9 + :: Source + ++ normal-parts + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (normal:rd:dist:i754rand r) + =^ v r (normal:rd:dist:i754rand r) + [(~(pak cd:complex %n) u v) r] + :: +cnormal: rng -> [@cd rng] + :: + :: Standard circularly-symmetric complex normal CN(0,1): components + :: N(0, 1/2), i.e. +normal-parts scaled by 1/sqrt(2). Deliberately + :: distinct from +normal-parts (N(0,1) components) -- "complex + :: Gaussian" means different things to signal-processing users (who + :: usually mean CN(0,1), total variance 1) and statistics users (who + :: usually mean iid N(0,1) parts, total variance 2). Both are named + :: unambiguously so neither reading silently gets the wrong one. + :: Examples + :: > `@ux`out:(cnormal:cd:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3fea.b5c4.88e8.bf41.bfe5.735e.01a0.2b7b + :: Source + ++ cnormal + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ z r (normal-parts r) + =/ u (mul:m invsqt2:m (~(re cd:complex %n) z)) + =/ v (mul:m invsqt2:m (~(im cd:complex %n) z)) + [(~(pak cd:complex %n) u v) r] + :: +on-circle: rng -> [@cd rng] + :: + :: Uniform on the unit circle: theta = tau*u for u drawn uniform + :: [0,1), components (cos theta, sin theta) via /lib/math's + :: bit-specified @rd kernels (jet parity holds, unlike this adapter's + :: own naive Taylor loops it borrows nothing from). + :: Examples + :: > `@ux`out:(on-circle:cd:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3fc0.7838.f0db.1ef6.3fef.bbe7.a757.e0e1 :: re^2+im^2 = 1 + :: Source + ++ on-circle + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (rd:uni:i754rand r) + =/ theta (mul:m tau:m u) + [(~(pak cd:complex %n) (cos:m theta) (sin:m theta)) r] + :: +in-disk: rng -> [@cd rng] + :: + :: Uniform on the unit disk via rejection from the enclosing square: + :: draw x,y uniform on (-1,1) (the same +rd:uni-mapped-2x-1 trick + :: +normal:rd:dist:i754rand uses for its own rejection setup), reject + :: if x^2+y^2 >= 1, else accept [x y] as the point. Variable- + :: consumption (rejection loop): acceptance probability pi/4 (~78.5%), + :: astronomically bounded in practice -- same rand-spec.md section 5 + :: policy as every other rejection-based arm in this codebase. + :: Examples + :: > `@ux`out:(in-disk:cd:complexrand (from-atom:seed:rand %sm64 0)) + :: 0xbfd1.1d5e.36cd.f850.bfe7.45ce.ffed.7562 :: re^2+im^2 < 1 + :: Source + ++ in-disk + |= r=rng:rand + ^- [out=@ r=rng:rand] + |- ^- [out=@ r=rng:rand] + =^ u1 r (rd:uni:i754rand r) + =^ u2 r (rd:uni:i754rand r) + =/ x (sub:m (mul:m .~2 u1) .~1) + =/ y (sub:m (mul:m .~2 u2) .~1) + ?: (gte:m (add:m (mul:m x x) (mul:m y y)) .~1) + $ + [(~(pak cd:complex %n) x y) r] + -- +:: +cs: complex-single (@cs) adapters. +:: +:: The @rs/@cs mirror of +cd, same five arms, same algorithms -- a +:: mechanical re-instantiation over /lib/math's rs door and /lib/ +:: i754rand's @rs generators, per rand-spec.md's "each arm exists at +:: @rd (reference) and @rs" convention (already used throughout +:: /lib/i754rand's own ++dist). +:: Source +++ cs + |% + ++ m ~(. rs:math [%n .1e-5 .0]) + :: +cuniform: rng -> [@cs rng] (see +cuniform:cd for algorithm) + :: Examples + :: > `@ux`out:(cuniform:cs:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3f39.65f4.3dee.6d78 + :: Source + ++ cuniform + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (rs:uni:i754rand r) + =^ v r (rs:uni:i754rand r) + [(~(pak cs:complex %n) u v) r] + :: +normal-parts: rng -> [@cs rng] (see +normal-parts:cd for algorithm) + :: Examples + :: > `@ux`out:(normal-parts:cs:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x4009.430e.bf17.e80f + :: Source + ++ normal-parts + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (normal:rs:dist:i754rand r) + =^ v r (normal:rs:dist:i754rand r) + [(~(pak cs:complex %n) u v) r] + :: +cnormal: rng -> [@cs rng] (see +cnormal:cd for algorithm) + :: Examples + :: > `@ux`out:(cnormal:cs:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3fc2.1e20.bed6.d405 + :: Source + ++ cnormal + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ z r (normal-parts r) + =/ u (mul:m invsqt2:m (~(re cs:complex %n) z)) + =/ v (mul:m invsqt2:m (~(im cs:complex %n) z)) + [(~(pak cs:complex %n) u v) r] + :: +on-circle: rng -> [@cs rng] (see +on-circle:cd for algorithm) + :: Examples + :: > `@ux`out:(on-circle:cs:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3f2b.0087.3f3e.82b8 :: re^2+im^2 = 1 + :: Source + ++ on-circle + |= r=rng:rand + ^- [out=@ r=rng:rand] + =^ u r (rs:uni:i754rand r) + =/ theta (mul:m tau:m u) + [(~(pak cs:complex %n) (cos:m theta) (sin:m theta)) r] + :: +in-disk: rng -> [@cs rng] (see +in-disk:cd for algorithm) + :: Examples + :: > `@ux`out:(in-disk:cs:complexrand (from-atom:seed:rand %sm64 0)) + :: 0x3ee5.97d0.bf44.64a2 :: re^2+im^2 < 1 + :: Source + ++ in-disk + |= r=rng:rand + ^- [out=@ r=rng:rand] + |- ^- [out=@ r=rng:rand] + =^ u1 r (rs:uni:i754rand r) + =^ u2 r (rs:uni:i754rand r) + =/ x (sub:m (mul:m .2 u1) .1) + =/ y (sub:m (mul:m .2 u2) .1) + ?: (gte:m (add:m (mul:m x x) (mul:m y y)) .1) + $ + [(~(pak cs:complex %n) x y) r] + -- +-- diff --git a/librand/desk/lib/fixedrand.hoon b/librand/desk/lib/fixedrand.hoon new file mode 100644 index 0000000..180df8b --- /dev/null +++ b/librand/desk/lib/fixedrand.hoon @@ -0,0 +1,83 @@ +/+ rand, twocrand, fx=fixed +:::: /lib/fixedrand -- fixed-point output adapter (rand-spec.md section +:::: 12.1) +:: +:: Routes /lib/rand's engines through /lib/fixed's two's-complement +:: fixed-point representation. No new engine work, and no floats: a +:: fixed-point lattice is UNIFORM, so uniform generation is exact by +:: construction (stronger than the float case in rand-spec.md section +:: 5.2 -- no rounding argument needed here at all). +:: +:: Imports /lib/fixed aliased to .fx, not .fixed: this file's own +:: +fixed arm (named to match rand-spec.md's public API) would +:: otherwise shadow the library face of the same name within its own +:: body and every later sibling arm's -- the exact footgun /lib/fixed +:: itself hit and fixed by renaming +twoc to +neg (see its own +neg +:: doc comment). Here the arm name is the spec-mandated one, so the +:: import is what gets renamed instead. +:: +:: Distributions (rand-spec.md 12.1: "sample at @rd, quantize round-to- +:: nearest") are NOT shipped as dedicated per-distribution wrapper arms +:: here -- that would be a dozen-plus nearly-identical one-liners doing +:: nothing beyond composing two already-existing arms. The composition +:: is exactly: draw a value from /lib/i754rand's ++dist, then quantize +:: via /lib/fixed's +from-rd (added alongside this file, mirroring +:: +from-rs, since it didn't exist yet). See tests/lib/fixedrand.hoon +:: for that composition proven end to end, and librand/NEXT-STEPS.md +:: for the reasoning. +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:: +fixed: [r=rng p=prec:fx] -> [@ rng] +:: +:: Full-range uniform: draw N=(wid p) raw bits, return as the two's- +:: complement pattern at that width. Like +twoc-full, this is the +:: IDENTITY on the raw bits -- a fixed-point value's on-the-wire +:: representation at a given precision IS just a two's-complement bit +:: pattern, so a uniform draw over raw bits already is a uniform draw +:: over fixed-point values at that precision. +:: Examples +:: > (fixed:fixedrand (from-atom:seed:rand %sm64 0) [8 8]) +:: [out=118.191 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] +:: Source +++ fixed + |= [r=rng:rand p=prec:fx] + ^- [out=@ r=rng:rand] + (bits:uni:rand r (wid:fx p)) +:: +fixed-unit: [r=rng p=prec:fx] -> [@ rng] +:: +:: Uniform in [0,1): draw .b.p bits as the fraction field, with the +:: integer field (and sign) forced to zero. Since the upper bits of a +:: Hoon atom are implicitly zero, drawing exactly .b.p bits already +:: produces the correct full-width pattern -- no shifting or masking +:: needed on top. +:: Examples +:: > (fixed-unit:fixedrand (from-atom:seed:rand %sm64 0) [8 8]) +:: [out=175 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] +:: Source +++ fixed-unit + |= [r=rng:rand p=prec:fx] + ^- [out=@ r=rng:rand] + (bits:uni:rand r b.p) +:: +fixed-between: [r=rng p=prec:fx a=@ b=@] -> [@ rng] +:: +:: Inclusive range [a,b], where .a and .b are width-N fixed-point +:: PATTERNS (N = wid p) in the SAME sense +twoc-between's .a/.b are raw +:: two's-complement patterns, not decoded values -- convert numeric +:: bounds to patterns yourself first (e.g. via /lib/fixed's +from-s/ +:: +from-rd) if that's what you have. A thin delegation to +:: +twoc-between at width N: fixed-point's bit-level representation IS +:: two's-complement, so once the width matches there's nothing left to +:: "reinterpret" -- the same bits mean a fixed-point value instead of a +:: plain integer purely by the caller's own convention. Crashes unless +:: a<=b in twoc order (see +twoc-between). +:: Examples +:: > (fixed-between:fixedrand (from-atom:seed:rand %sm64 0) [8 8] 0x1.ff00 0x300) +:: :: range [-1.0, 3.0] in q8.8 +:: [out=649 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] +:: Source +++ fixed-between + |= [r=rng:rand p=prec:fx a=@ b=@] + ^- [out=@ r=rng:rand] + (twoc-between:twocrand r (wid:fx p) a b) +-- diff --git a/librand/desk/lib/i754rand.hoon b/librand/desk/lib/i754rand.hoon new file mode 100644 index 0000000..b3aeba7 --- /dev/null +++ b/librand/desk/lib/i754rand.hoon @@ -0,0 +1,985 @@ +/+ rand, math +:::: /lib/i754rand -- IEEE-754 float generation + distributions +:: +:: The porcelain layer on top of /lib/rand's plumbing: uniform floats +:: (++uni: +rs/+rd/+rh/+rq/+rs-oo/+rd-oo), Vose's alias method (++alias, +:: moved here from /lib/rand's ++sample since +draw needs an actual +:: uniform @rd float draw), and the continuous/discrete distributions +:: built on both (++dist). Named for the numeric domain it serves (IEEE +:: 754), mirroring how /lib/twoc/fixed/complex/unum are kept separate +:: from /lib/math -- this is /lib/rand's own "math.hoon". Non-float +:: output adapters (twocrand, fixedrand, complexrand, unumrand) are this +:: file's siblings, not its dependents: none of them need floats for +:: their own core uniform-generation arms, only (optionally) for a +:: "sample at @rd, then quantize/convert" distribution wrapper -- see +:: NEXT-STEPS.md. +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:::: ++uni :: (5.2) float generation +:: +:: NAMING FOOTGUN (same class as /lib/rand's ++seed +mix): this core +:: defines arms named +rs/+rd/+rh/+rq, which shadow the Hoon standard +:: library's own +rs/+rd/+rh/+rq IEEE-754 doors. /lib/math has the exact +:: same collision (its own +rs/+rd/+rh/+rq doors wrap the stdlib ones) and +:: handles it the same way: any reference to the real stdlib door from +:: inside ++uni uses ^rs/^rd/^rh/^rq, e.g. `(~(mul ^rs %n) a b)`. This +:: core also shadows /lib/rand's OWN ++uni (the integer-only one) by +:: name -- they're siblings in DIFFERENT files, not nested, so there is +:: no actual collision, but the float generators below reach the integer +:: primitives they're built on via `bits:uni:rand` (qualified through the +:: imported `rand` face), never bare `bits`. +:: +++ uni + |% + :: +rs: rng -> [@rs rng] + :: + :: Uniform in [0,1): draw 24 bits, multiply by the exact constant + :: 2^-24. +sun of a 24-bit unsigned integer and multiplying by an + :: exact power of two are BOTH exact IEEE-754 operations (no rounding, + :: regardless of mode), so every representable output is a multiple of + :: 2^-24 -- exactly uniform over that lattice. 0 is possible; 1 is not. + :: Examples + :: > (rs:uni:i754rand (from-atom:seed:rand %sm64 0)) + :: [out=.0.1164197 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] + :: Source + ++ rs + |= r=rng:rand + ^- [out=@rs r=rng:rand] + =^ b r (bits:uni:rand r 24) + [(~(mul ^rs %n) (~(sun ^rs %n) b) `@rs`0x3380.0000) r] + :: +rd: rng -> [@rd rng] + :: + :: Same as +rs at double precision: 53 bits, times the exact constant + :: 2^-53. + :: Source + ++ rd + |= r=rng:rand + ^- [out=@rd r=rng:rand] + =^ b r (bits:uni:rand r 53) + [(~(mul ^rd %n) (~(sun ^rd %n) b) `@rd`0x3ca0.0000.0000.0000) r] + :: +rh: rng -> [@rh rng] + :: + :: Same as +rs at half precision: 11 bits, times the exact constant + :: 2^-11. + :: Source + ++ rh + |= r=rng:rand + ^- [out=@rh r=rng:rand] + =^ b r (bits:uni:rand r 11) + [(~(mul ^rh %n) (~(sun ^rh %n) b) `@rh`0x1000) r] + :: +rq: rng -> [@rq rng] + :: + :: Same as +rs at quad precision: 113 bits, times the exact constant + :: 2^-113. + :: Source + ++ rq + |= r=rng:rand + ^- [out=@rq r=rng:rand] + =^ b r (bits:uni:rand r 113) + [(~(mul ^rq %n) (~(sun ^rq %n) b) `@rq`0x3f8e.0000.0000.0000.0000.0000.0000.0000) r] + :: +rs-oo: rng -> [@rs rng] + :: + :: Uniform in the OPEN interval (0,1): draw the same 24 bits as +rs, but + :: construct (2*bits+1) * 2^-25 instead of bits * 2^-24. With bits + :: ranging over [0,2^24), (2*bits+1) ranges over the odd integers in + :: [1,2^25), so the result spans (0,1) on a lattice of spacing 2^-24, + :: symmetric about 0.5, never touching either endpoint. Needed by + :: Box-Muller/log-based transforms (milestone 5) so log(0) never fires. + :: (2*bits+1) is exact @ arithmetic; the multiply by 2^-25 is exact for + :: the same reason it is in +rs. + :: Source + ++ rs-oo + |= r=rng:rand + ^- [out=@rs r=rng:rand] + =^ b r (bits:uni:rand r 24) + [(~(mul ^rs %n) (~(sun ^rs %n) +((mul 2 b))) `@rs`0x3300.0000) r] + :: +rd-oo: rng -> [@rd rng] + :: + :: Same construction as +rs-oo at double precision: 53 bits, (2*bits+1) + :: * 2^-54. + :: Source + ++ rd-oo + |= r=rng:rand + ^- [out=@rd r=rng:rand] + =^ b r (bits:uni:rand r 53) + [(~(mul ^rd %n) (~(sun ^rd %n) +((mul 2 b))) `@rd`0x3c90.0000.0000.0000) r] + -- +:: +:::: ++alias :: (6.1) alias method +:: +:: Vose's alias method, moved here from /lib/rand's ++sample: +draw needs +:: an actual uniform @rd float draw to decide accept-vs-redirect, so it +:: is not float-free the way the rest of ++sample is. +:: +:: $alias-table: [n=@ prob=(list @rd) alias=(list @ud)] +:: +:: Vose's alias method (1991) table: O(n) to build, O(1) to draw from. +:: .prob and .alias are indexed 0..n-1, parallel to the input weight +:: list's own order. ++$ alias-table [n=@ prob=(list @rd) alias=(list @ud)] +:: +++ alias + |% + :: +mth: the shared @rd math door for +build's arithmetic, forced %n + :: per section 6.1's "require @rd %n arithmetic only." + ++ mth ~(. rd:math [%n .~1e-13 .~0]) + :: +build: (list @rd) -> alias-table + :: + :: Vose's O(n) construction. Pure (no rng). Weights need not sum to + :: 1 -- they're normalized here (scaled so their mean is 1, which is + :: what Vose's small/large partition compares against). Weights must + :: be non-negative with at least one positive; crashes otherwise. + :: + :: The small/large worklist (the classic place for off-by-one and + :: float-compare bugs, per section 6.1): +build-loop pops one + :: under-weight index (l, scaled prob < 1) and one over-weight index + :: (g, scaled prob >= 1) each iteration. l's own probability becomes + :: its final +prob entry, and l's +alias entry is set to g -- so + :: drawing l falls through to g whenever the per-draw uniform check + :: fails (see +draw). g gives up exactly the mass l was short by + :: (1 - wt.l), and whatever's left of g's weight goes back on + :: whichever worklist it now belongs to. When either worklist runs + :: dry, everything left on the other is a floating-point leftover + :: (mathematically it should be exactly weight 1, off by rounding) -- + :: +build-cleanup gives each remaining index probability 1 outright, + :: with its own index as a harmless unused alias. + :: Examples + :: > (build:alias:i754rand ~[.~1 .~1 .~2]) + :: [n=3 prob=~[.~0.75 .~0.75 .~1] alias=~[2 2 2]] + :: Source + ++ build + |= weights=(list @rd) + ^- alias-table + =/ n (lent weights) + ~| %rand-empty-list + ?> (gth n 0) + =/ total (roll weights add:mth) + ~| %rand-bad-prob + ?> (gth:mth total .~0) + =/ scale (div:mth (sun:mth n) total) + =/ scaled (turn weights |=(w=@rd (mul:mth w scale))) + =/ idxd ^- (list [idx=@ud wt=@rd]) + (turn (gulf 0 (dec n)) |=(k=@ud [idx=k wt=(snag k scaled)])) + =/ parted (skid idxd |=([idx=@ud wt=@rd] (lth:mth wt .~1))) + =/ fin (build-loop -.parted +.parted *(map @ud @rd) *(map @ud @ud)) + :+ n + (turn (gulf 0 (dec n)) |=(k=@ud (~(got by probm.fin) k))) + (turn (gulf 0 (dec n)) |=(k=@ud (~(got by aliasm.fin) k))) + ++ build-loop + |= $: small=(list [idx=@ud wt=@rd]) + large=(list [idx=@ud wt=@rd]) + probm=(map @ud @rd) + aliasm=(map @ud @ud) + == + ^- [probm=(map @ud @rd) aliasm=(map @ud @ud)] + ?: |(?=(~ small) ?=(~ large)) + (build-cleanup (weld small large) probm aliasm) + =/ l i.small + =/ g i.large + =. probm (~(put by probm) idx.l wt.l) + =. aliasm (~(put by aliasm) idx.l idx.g) + =/ pg (sub:mth (add:mth wt.g wt.l) .~1) + :: pg<1: g is now under-weight -> joins small. pg>=1: g stays + :: over-weight -> joins large. (Reversing these was a real bug + :: caught during on-ship verification: prob/alias came out wrong + :: for a 3-element table, traced back to g landing in the opposite + :: worklist from where its own weight said it belonged.) + ?: (lth:mth pg .~1) + (build-loop [[idx=idx.g wt=pg] t.small] t.large probm aliasm) + (build-loop t.small [[idx=idx.g wt=pg] t.large] probm aliasm) + ++ build-cleanup + |= $: rest=(list [idx=@ud wt=@rd]) + probm=(map @ud @rd) + aliasm=(map @ud @ud) + == + ^- [probm=(map @ud @rd) aliasm=(map @ud @ud)] + ?~ rest + [probm aliasm] + %= $ + rest t.rest + probm (~(put by probm) idx.i.rest .~1) + aliasm (~(put by aliasm) idx.i.rest idx.i.rest) + == + :: +draw: [t=alias-table r=rng] -> [out=@ud rng] + :: + :: O(1): pick an index uniformly, then a coin flip (biased by + :: .prob.t at that index) decides whether to keep it or redirect to + :: its alias. + :: Source + ++ draw + |= [t=alias-table r=rng:rand] + ^- [out=@ud r=rng:rand] + =^ i r (below:uni:rand r n.t) + =^ u r (rd:uni r) + :- ?:((lth:mth u (snag i prob.t)) i (snag i alias.t)) + r + -- +:: +:::: ++dist :: (6) distributions +:: +:: Nested ++rd / ++rs sub-cores, one per precision (see each for detail). +:: Internally forced to round-to-nearest (%n) via +m in each, matching +:: /lib/math's own doors -- callers never see or choose a rounding mode +:: here. chi2/student-t aren't named in rand-spec.md's milestone punch +:: list but are one-line compositions of gamma/normal with no new +:: machinery, so they land here too rather than waiting on an unscheduled +:: slot. +:: +:: Crashes (`?>` with `~|` tags) on out-of-domain parameters, per section +:: 9's uniform error policy: domain violations are programmer error, not +:: data, so no NaN-returning "soft" errors. +:: +++ dist + |% + :: ++rd: the reference distributions, at @rd (double precision) + :: + :: Everything rand-spec.md section 6 specifies, at the reference + :: precision. See ++rs below for the routed-through-the-same- + :: algorithm single-precision mirror (rand-spec.md: "each arm exists + :: at @rd (reference) and @rs"). + ++ rd + |% + :: +m: the shared @rd math door, forced %n, rtol=1e-13 (the same "sane + :: default" convergence tolerance Saloon's own +feps uses for @rd). + :: Internal helper, mirroring /lib/fixed's +ng pattern. + ++ m ~(. rd:math [%n .~1e-13 .~0]) + :: +mz: same door, forced %z (truncate toward zero) -- used only by + :: +geometric's ceiling, since +toi respects the door's rounding mode + :: and %n would round instead of truncate. + ++ mz ~(. rd:math [%z .~1e-13 .~0]) + :: +normal: rng -> [@rd rng] + :: + :: Standard normal N(0,1) via the Marsaglia polar method (the Box-Muller + :: polar variant): draw u,v uniform on (-1,1) (via +rd:uni mapped + :: 2x-1), reject if s=u^2+v^2 is >=1 or =0, else return u*sqrt(-2 ln(s)/s). + :: Returns ONE deviate and discards the pair-mate v*sqrt(-2 ln(s)/s) -- + :: caching it would make the rng state opaque (the next +normal call + :: would need to "remember" a pending mate outside the plain +$rng + :: noun), so v's factor is thrown away at v1. A documented waste, not a + :: bug; Ziggurat is a NEXT-STEPS optimization. Variable-consumption + :: (rejection loop): expected iterations ~1.27 (rejection probability + :: 1 - pi/4), astronomically bounded in practice. + :: Examples + :: > (normal:rd:dist:i754rand (from-atom:seed:rand %sm64 0)) + :: [out=.~-0.9479938949723624 r=[%sm64 s=0x78dd.e6e5.fd29.f054]] + :: Source + ++ normal + |= r=rng:rand + ^- [out=@rd r=rng:rand] + |- ^- [out=@rd r=rng:rand] + =^ u1 r (rd:uni r) + =^ u2 r (rd:uni r) + =/ u (sub:m (mul:m .~2 u1) .~1) + =/ v (sub:m (mul:m .~2 u2) .~1) + =/ s (add:m (mul:m u u) (mul:m v v)) + ?: |((gte:m s .~1) (equ:m s .~0)) + $ + =/ factor (sqt:m (div:m (mul:m .~-2 (log:m s)) s)) + [(mul:m u factor) r] + :: +normal-mv: [r=rng mu=@rd sigma=@rd] -> [@rd rng] + :: + :: N(mu, sigma^2): mu + sigma*z where z ~ N(0,1). Crashes if sigma < 0. + :: Source + ++ normal-mv + |= [r=rng:rand mu=@rd sigma=@rd] + ^- [out=@rd r=rng:rand] + ~| %rand-bad-sigma + ?> !(lth:m sigma .~0) + =^ z r (normal r) + [(add:m mu (mul:m sigma z)) r] + :: +expon: [r=rng lambda=@rd] -> [@rd rng] + :: + :: Exponential(lambda) via inversion: -ln(u)/lambda, u drawn from the + :: OPEN (0,1) (+rd-oo, not +rd) specifically so log(0) never fires -- + :: this is exactly the case rand-spec.md section 5.2 built +rd-oo for. + :: Crashes if lambda <= 0. + :: Source + ++ expon + |= [r=rng:rand lambda=@rd] + ^- [out=@rd r=rng:rand] + ~| %rand-bad-rate + ?> (gth:m lambda .~0) + =^ u r (rd-oo:uni r) + [(div:m (neg:m (log:m u)) lambda) r] + :: +gamma: [r=rng alpha=@rd] -> [@rd rng] + :: + :: Gamma(alpha, scale=1) via Marsaglia-Tsang (2000). alpha>=1 direct + :: (+gamma-ge1); alpha<1 via the standard boost gamma(alpha) = + :: gamma(alpha+1) * u^(1/alpha), u drawn from the open (0,1) so the + :: u=0 lattice point (probability 2^-53, not truly 0 as it would be for + :: a continuous uniform) never manufactures a spurious exact-zero + :: sample. Crashes if alpha <= 0. + :: Source + ++ gamma + |= [r=rng:rand alpha=@rd] + ^- [out=@rd r=rng:rand] + ~| %rand-bad-shape + ?> (gth:m alpha .~0) + ?: (gte:m alpha .~1) + (gamma-ge1 r alpha) + =^ g r (gamma-ge1 r (add:m alpha .~1)) + =^ u r (rd-oo:uni r) + [(mul:m g (pow:m u (div:m .~1 alpha))) r] + :: +gamma-ge1: Marsaglia-Tsang squeeze for alpha>=1. d=alpha-1/3, + :: c=1/sqrt(9d); draw x~N(0,1), v=(1+cx)^3 (reject if v<=0), draw + :: u~(0,1) open (so log(u) never fires on 0), accept d*v if + :: ln(u) < x^2/2 + d - d*v + d*ln(v), else reject and redraw both x,u. + :: Variable-consumption (rejection loop), astronomically bounded. + ++ gamma-ge1 + |= [r=rng:rand alpha=@rd] + ^- [out=@rd r=rng:rand] + =/ d (sub:m alpha (div:m .~1 .~3)) + =/ c (div:m .~1 (sqt:m (mul:m .~9 d))) + |- ^- [out=@rd r=rng:rand] + =^ x r (normal r) + =/ t (add:m .~1 (mul:m c x)) + =/ v (mul:m t (mul:m t t)) + ?: !(gth:m v .~0) + $ + =^ u r (rd-oo:uni r) + =/ rhs + %+ add:m + (add:m (mul:m .~0.5 (mul:m x x)) d) + (sub:m (mul:m d (log:m v)) (mul:m d v)) + ?: (lth:m (log:m u) rhs) + [(mul:m d v) r] + $ + :: +beta: [r=rng a=@rd b=@rd] -> [@rd rng] + :: + :: Beta(a,b) via two independent gammas: x/(x+y), x~gamma(a), y~gamma(b). + :: Crashes if a <= 0 or b <= 0 (via +gamma's own precondition). + :: Source + ++ beta + |= [r=rng:rand a=@rd b=@rd] + ^- [out=@rd r=rng:rand] + =^ x r (gamma r a) + =^ y r (gamma r b) + [(div:m x (add:m x y)) r] + :: +chi2: [r=rng k=@] -> [@rd rng] + :: + :: Chi-squared with k degrees of freedom: gamma(k/2, scale=2) == + :: 2*gamma(k/2, scale=1) (+gamma is scale=1, so the factor of 2 is + :: applied directly -- a standard gamma scaling property). Crashes if + :: k = 0. + :: Source + ++ chi2 + |= [r=rng:rand k=@] + ^- [out=@rd r=rng:rand] + ~| %rand-bad-df + ?> !=(k 0) + =^ x r (gamma r (div:m (sun:m k) .~2)) + [(mul:m .~2 x) r] + :: +student-t: [r=rng k=@] -> [@rd rng] + :: + :: Student's t with k degrees of freedom: z / sqrt(chi2(k)/k), z~N(0,1). + :: Crashes if k = 0. + :: Source + ++ student-t + |= [r=rng:rand k=@] + ^- [out=@rd r=rng:rand] + ~| %rand-bad-df + ?> !=(k 0) + =^ z r (normal r) + =^ c r (chi2 r k) + [(div:m z (sqt:m (div:m c (sun:m k)))) r] + :: +bernoulli: [r=rng p=@rd] -> [? rng] + :: + :: %.y with probability p, else %.n: draw u~[0,1), return u &((gte:m p .~0) (lte:m p .~1)) + =^ u r (rd:uni r) + [(lth:m u p) r] + :: +geometric: [r=rng p=@rd] -> [@ud rng] + :: + :: Number of Bernoulli(p) trials up to and including the first success: + :: ceil(ln(u)/ln(1-p)), u drawn from the open (0,1) (avoids ln(0)). + :: p=1 is a special-cased edge that returns 1 directly (ln(1-p) would + :: divide by ln(0) = -inf otherwise). ceil is computed via +toi under + :: a SEPARATE %z-rounding (truncate-toward-zero) door instance, +mz -- + :: +toi respects the door's own rounding mode, and this core's shared + :: +m door is forced %n (round-to-nearest), which would round instead + :: of truncate. Crashes unless 0 < p <= 1. + :: Source + ++ geometric + |= [r=rng:rand p=@rd] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-prob + ?> &((gth:m p .~0) (lte:m p .~1)) + ?: (equ:m p .~1) + [1 r] + =^ u r (rd-oo:uni r) + =/ raw (div:m (log:m u) (log:m (sub:m .~1 p))) + =/ fl (abs:si (need (toi:mz raw))) + [?:(=(raw (sun:m fl)) fl +(fl)) r] + :: +categorical: [r=rng t=alias-table] -> [@ud rng] + :: + :: One draw from a pre-built alias table (rand-spec.md section 6.1's + :: ++alias, a sibling core in this file). Deliberately takes an + :: ALREADY-BUILT table, + :: not raw weights: alias's whole point is O(n) build, O(1) per draw, + :: amortized over many draws from the same distribution -- rebuilding + :: the table on every single draw would defeat that. Build once via + :: `(build:alias weights)`, draw many times via this arm (or + :: `draw:alias` directly, which this is a thin rename of). + :: Lives ONLY here, not mirrored in ++rs: alias-table's .prob is fixed + :: at @rd (rand-spec.md section 6.1), so this arm is precision- + :: invariant -- an "@rs categorical" would be byte-for-byte identical + :: code, not a real variant. + :: Source + ++ categorical + |= [r=rng:rand t=alias-table] + ^- [out=@ud r=rng:rand] + (draw:alias t r) + :: +poisson: [r=rng lambda=@rd] -> [@ud rng] + :: + :: Poisson(lambda): Knuth's product method for lambda<10 (simple, + :: O(lambda) expected multiplications -- fine at this scale, unusable + :: above it); Hörmann's PTRS (1993) transformed rejection for + :: lambda>=10 (O(1) expected, needed because Knuth's method's cost + :: scales linearly with lambda and would silently become the wrong + :: choice above the threshold rather than crashing -- so this arm + :: switches automatically rather than leaving the choice to the + :: caller). Crashes if lambda <= 0. + :: Source + ++ poisson + |= [r=rng:rand lambda=@rd] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-rate + ?> (gth:m lambda .~0) + ?: (lth:m lambda .~10) + (poisson-knuth r lambda) + (poisson-ptrs r lambda) + :: +poisson-knuth: k=0, p=1; loop: k+=1, p*=uniform(0,1); accept + :: (return k-1, i.e. the count BEFORE the last increment) once + :: p <= exp(-lambda). Variable-consumption (rejection-free here, but + :: a variable number of multiplications -- expected lambda of them). + ++ poisson-knuth + |= [r=rng:rand lambda=@rd] + ^- [out=@ud r=rng:rand] + =/ bigl (exp:m (neg:m lambda)) + =/ k 0 + =/ p .~1 + |- ^- [out=@ud r=rng:rand] + =^ u r (rd:uni r) + =/ p2 (mul:m p u) + ?: (lte:m p2 bigl) + [k r] + $(k +(k), p p2) + :: +poisson-ptrs: Hörmann 1993's transformed rejection with squeeze, + :: verified against NumPy's random_poisson_ptrs (src/distributions/ + :: distributions.c) -- an independent re-implementation of the exact + :: same reference this spec cites, not reconstructed from memory. + :: Draws a candidate k from a scaled/shifted uniform (the "transformed" + :: part), accepts immediately if a cheap squeeze test passes (avoiding + :: the exact log-probability computation on the common path), and + :: falls back to the exact accept/reject test (via +loggam, log(k!)) + :: otherwise. Rejection loop, astronomically bounded in practice. + ++ poisson-ptrs + |= [r=rng:rand lambda=@rd] + ^- [out=@ud r=rng:rand] + =/ slam (sqt:m lambda) + =/ loglam (log:m lambda) + =/ b (add:m .~0.931 (mul:m .~2.53 slam)) + =/ a (add:m .~-0.059 (mul:m .~0.02483 b)) + =/ invalpha (add:m .~1.1239 (div:m .~1.1328 (sub:m b .~3.4))) + =/ vr (sub:m .~0.9277 (div:m .~3.6224 (sub:m b .~2))) + |- ^- [out=@ud r=rng:rand] + =^ u1 r (rd:uni r) + =^ v r (rd:uni r) + =/ u (sub:m u1 .~0.5) + =/ us (sub:m .~0.5 (abs:m u)) + =/ kk (dfloor (add:m (add:m (mul:m (add:m (div:m (mul:m .~2 a) us) b) u) lambda) .~0.43)) + ?: &((gte:m us .~0.07) (lte:m v vr)) + [(abs:si kk) r] + ?: |(=(-1 (cmp:si kk --0)) &((lth:m us .~0.013) (gth:m v us))) + $ + =/ lhs (sub:m (add:m (log:m v) (log:m invalpha)) (log:m (add:m (div:m a (mul:m us us)) b))) + =/ rhs (sub:m (add:m (mul:m .~-1 lambda) (mul:m (san:m kk) loglam)) (loggam (add:m (san:m kk) .~1))) + ?: (lte:m lhs rhs) + [(abs:si kk) r] + $ + :: +dfloor: @rd -> @s, floor (round toward -infinity, NOT toward + :: zero) -- +poisson-ptrs's candidate k can be negative, where + :: truncate-toward-zero (what +toi under +mz gives) would round the + :: wrong way. Positive/exact values: truncation IS the floor. + :: Negative non-integers: truncation rounds up, so subtract 1. + ++ dfloor + |= x=@rd + ^- @s + =/ t (need (toi:mz x)) + ?: |((gte:m x .~0) (equ:m x (san:m t))) + t + (dif:si t --1) + :: +loggam: @rd -> @rd, natural log of the gamma function (so log(k!) + :: = loggam(k+1)). Stirling's series with a small-argument recurrence + :: for x<7 (shifts x up by whole integers until >=7, where the series + :: is accurate, then subtracts the shifted-away log-factors back out) + :: -- transcribed from NumPy's random_loggam (src/distributions/ + :: distributions.c), itself from Zhang & Jin's SPECFUN algorithm, not + :: derived independently: log-gamma accuracy is exactly the kind of + :: detail worth matching a vetted reference bit-for-bit rather than + :: reconstructing from the general shape of Stirling's series. + ++ loggam + |= x=@rd + ^- @rd + ?: |((equ:m x .~1) (equ:m x .~2)) + .~0 + =/ a ^- (list @rd) + :~ .~0.08333333333333333 .~-0.002777777777777778 .~0.0007936507936507937 + .~-0.0005952380952380952 .~0.0008417508417508417 .~-0.001917526917526918 + .~0.006410256410256411 .~-0.02955065359477124 .~0.1796443723688307 + .~-1.39243221690590 + == + =/ lg2pi .~1.8378770664093453 + =/ n ?:((lth:m x .~7) (abs:si (need (toi:mz (sub:m .~7 x)))) 0) + =/ x0 (add:m x (sun:m n)) + =/ x2 (div:m .~1 (mul:m x0 x0)) + =/ gl0 + =/ acc (snag 9 a) + =/ i 0 + |- ^- @rd + ?: =(i 9) + acc + $(acc (add:m (mul:m acc x2) (snag (sub 8 i) a)), i +(i)) + =/ gl + %+ sub:m + %+ add:m + (add:m (div:m gl0 x0) (mul:m .~0.5 lg2pi)) + (mul:m (sub:m x0 .~0.5) (log:m x0)) + x0 + ?. (lth:m x .~7) + gl + =/ fin + =/ gl2 gl + =/ x02 x0 + =/ k 1 + |- ^- @rd + ?: (gth k n) + gl2 + $(gl2 (sub:m gl2 (log:m (sub:m x02 .~1))), x02 (sub:m x02 .~1), k +(k)) + fin + :: +binomial: [r=rng n=@ p=@rd] -> [@ud rng] + :: + :: Binomial(n,p) via inversion by recurrence: P(X=0)=(1-p)^n, then + :: P(X=k+1) = P(X=k) * (n-k)/(k+1) * p/(1-p), accumulating the CDF + :: until it exceeds a single uniform draw. Only valid for + :: n*min(p,1-p) < 30 (Kachitvichyanukul & Schmeiser 1988) -- above + :: that, inversion needs too many terms and BTPE is the right + :: algorithm, deferred to NEXT-STEPS; this arm crashes rather than + :: running slow (many terms) or, worse, silently truncating. + :: Source + ++ binomial + |= [r=rng:rand n=@ p=@rd] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-prob + ?> &((gte:m p .~0) (lte:m p .~1)) + =/ q (sub:m .~1 p) + =/ qmin ?:((lth:m p .~0.5) p q) + ~| %rand-binomial-n-too-large + ?> (lth:m (mul:m (sun:m n) qmin) .~30) + =^ u r (rd:uni r) + =/ f (pow:m q (sun:m n)) + =/ cdf f + =/ k 0 + |- ^- [out=@ud r=rng:rand] + ?: (lth:m u cdf) + [k r] + =/ f2 %+ mul:m f + (mul:m (div:m (sun:m (sub n k)) (sun:m +(k))) (div:m p q)) + $(f f2, cdf (add:m cdf f2), k +(k)) + :: +dirichlet: [r=rng alphas=(list @rd)] -> [(list @rd) rng] + :: + :: Dirichlet(alphas): draw one gamma(alpha_i) per entry, normalize by + :: their sum. Crashes on an empty .alphas (and, via +gamma's own + :: precondition, on any non-positive alpha). + :: Source + ++ dirichlet + |= [r=rng:rand alphas=(list @rd)] + ^- [out=(list @rd) r=rng:rand] + ~| %rand-empty-list + ?> !=(~ alphas) + =/ gd + |- ^- [(list @rd) rng:rand] + ?~ alphas + [~ r] + =^ x r (gamma r i.alphas) + =^ rest r $(alphas t.alphas) + [[x rest] r] + =/ total (roll -.gd add:m) + [(turn -.gd |=(g=@rd (div:m g total))) +.gd] + -- + :: ++rs: the same distributions, routed through @rs (single + :: precision) + :: + :: Identical algorithms to ++rd, just built on ++uni's @rs generators + :: and /lib/math's @rs door instead of @rd's -- not a separate + :: derivation, a mechanical re-instantiation at the other precision + :: (mirroring how /lib/math itself has independent per-precision + :: doors rather than one generic implementation). +categorical has + :: no @rs mirror: alias-table's .prob is fixed at @rd, so it's + :: precision-invariant (see its doccord in ++rd). + ++ rs + |% + :: +m: the shared @rs math door, forced %n, rtol=1e-6 (the same "sane + :: default" convergence tolerance Saloon's own +feps uses for @rs). + :: Internal helper, mirroring /lib/fixed's +ng pattern. + ++ m ~(. rs:math [%n .1e-6 .0]) + :: +mz: same door, forced %z (truncate toward zero) -- used only by + :: +geometric's ceiling, since +toi respects the door's rounding mode + :: and %n would round instead of truncate. + ++ mz ~(. rs:math [%z .1e-6 .0]) + :: +normal: rng -> [@rs rng] + :: + :: Standard normal N(0,1) via the Marsaglia polar method (the Box-Muller + :: polar variant): draw u,v uniform on (-1,1) (via +rs:uni mapped + :: 2x-1), reject if s=u^2+v^2 is >=1 or =0, else return u*sqrt(-2 ln(s)/s). + :: Returns ONE deviate and discards the pair-mate v*sqrt(-2 ln(s)/s) -- + :: caching it would make the rng state opaque (the next +normal call + :: would need to "remember" a pending mate outside the plain +$rng + :: noun), so v's factor is thrown away at v1. A documented waste, not a + :: bug; Ziggurat is a NEXT-STEPS optimization. Variable-consumption + :: (rejection loop): expected iterations ~1.27 (rejection probability + :: 1 - pi/4), astronomically bounded in practice. + :: Examples + :: > (normal:rs:dist:i754rand (from-atom:seed:rand %sm64 0)) + :: [out=.-0.5933847 r=[%sm64 s=0x3c6e.f372.fe94.f82a]] + :: Source + ++ normal + |= r=rng:rand + ^- [out=@rs r=rng:rand] + |- ^- [out=@rs r=rng:rand] + =^ u1 r (rs:uni r) + =^ u2 r (rs:uni r) + =/ u (sub:m (mul:m .2 u1) .1) + =/ v (sub:m (mul:m .2 u2) .1) + =/ s (add:m (mul:m u u) (mul:m v v)) + ?: |((gte:m s .1) (equ:m s .0)) + $ + =/ factor (sqt:m (div:m (mul:m .-2 (log:m s)) s)) + [(mul:m u factor) r] + :: +normal-mv: [r=rng mu=@rs sigma=@rs] -> [@rs rng] + :: + :: N(mu, sigma^2): mu + sigma*z where z ~ N(0,1). Crashes if sigma < 0. + :: Source + ++ normal-mv + |= [r=rng:rand mu=@rs sigma=@rs] + ^- [out=@rs r=rng:rand] + ~| %rand-bad-sigma + ?> !(lth:m sigma .0) + =^ z r (normal r) + [(add:m mu (mul:m sigma z)) r] + :: +expon: [r=rng lambda=@rs] -> [@rs rng] + :: + :: Exponential(lambda) via inversion: -ln(u)/lambda, u drawn from the + :: OPEN (0,1) (+rd-oo, not +rd) specifically so log(0) never fires -- + :: this is exactly the case rand-spec.md section 5.2 built +rd-oo for. + :: Crashes if lambda <= 0. + :: Source + ++ expon + |= [r=rng:rand lambda=@rs] + ^- [out=@rs r=rng:rand] + ~| %rand-bad-rate + ?> (gth:m lambda .0) + =^ u r (rs-oo:uni r) + [(div:m (neg:m (log:m u)) lambda) r] + :: +gamma: [r=rng alpha=@rs] -> [@rs rng] + :: + :: Gamma(alpha, scale=1) via Marsaglia-Tsang (2000). alpha>=1 direct + :: (+gamma-ge1); alpha<1 via the standard boost gamma(alpha) = + :: gamma(alpha+1) * u^(1/alpha), u drawn from the open (0,1) so the + :: u=0 lattice point (probability 2^-53, not truly 0 as it would be for + :: a continuous uniform) never manufactures a spurious exact-zero + :: sample. Crashes if alpha <= 0. + :: Source + ++ gamma + |= [r=rng:rand alpha=@rs] + ^- [out=@rs r=rng:rand] + ~| %rand-bad-shape + ?> (gth:m alpha .0) + ?: (gte:m alpha .1) + (gamma-ge1 r alpha) + =^ g r (gamma-ge1 r (add:m alpha .1)) + =^ u r (rs-oo:uni r) + [(mul:m g (pow:m u (div:m .1 alpha))) r] + :: +gamma-ge1: Marsaglia-Tsang squeeze for alpha>=1. d=alpha-1/3, + :: c=1/sqrt(9d); draw x~N(0,1), v=(1+cx)^3 (reject if v<=0), draw + :: u~(0,1) open (so log(u) never fires on 0), accept d*v if + :: ln(u) < x^2/2 + d - d*v + d*ln(v), else reject and redraw both x,u. + :: Variable-consumption (rejection loop), astronomically bounded. + ++ gamma-ge1 + |= [r=rng:rand alpha=@rs] + ^- [out=@rs r=rng:rand] + =/ d (sub:m alpha (div:m .1 .3)) + =/ c (div:m .1 (sqt:m (mul:m .9 d))) + |- ^- [out=@rs r=rng:rand] + =^ x r (normal r) + =/ t (add:m .1 (mul:m c x)) + =/ v (mul:m t (mul:m t t)) + ?: !(gth:m v .0) + $ + =^ u r (rs-oo:uni r) + =/ rhs + %+ add:m + (add:m (mul:m .0.5 (mul:m x x)) d) + (sub:m (mul:m d (log:m v)) (mul:m d v)) + ?: (lth:m (log:m u) rhs) + [(mul:m d v) r] + $ + :: +beta: [r=rng a=@rs b=@rs] -> [@rs rng] + :: + :: Beta(a,b) via two independent gammas: x/(x+y), x~gamma(a), y~gamma(b). + :: Crashes if a <= 0 or b <= 0 (via +gamma's own precondition). + :: Source + ++ beta + |= [r=rng:rand a=@rs b=@rs] + ^- [out=@rs r=rng:rand] + =^ x r (gamma r a) + =^ y r (gamma r b) + [(div:m x (add:m x y)) r] + :: +chi2: [r=rng k=@] -> [@rs rng] + :: + :: Chi-squared with k degrees of freedom: gamma(k/2, scale=2) == + :: 2*gamma(k/2, scale=1) (+gamma is scale=1, so the factor of 2 is + :: applied directly -- a standard gamma scaling property). Crashes if + :: k = 0. + :: Source + ++ chi2 + |= [r=rng:rand k=@] + ^- [out=@rs r=rng:rand] + ~| %rand-bad-df + ?> !=(k 0) + =^ x r (gamma r (div:m (sun:m k) .2)) + [(mul:m .2 x) r] + :: +student-t: [r=rng k=@] -> [@rs rng] + :: + :: Student's t with k degrees of freedom: z / sqrt(chi2(k)/k), z~N(0,1). + :: Crashes if k = 0. + :: Source + ++ student-t + |= [r=rng:rand k=@] + ^- [out=@rs r=rng:rand] + ~| %rand-bad-df + ?> !=(k 0) + =^ z r (normal r) + =^ c r (chi2 r k) + [(div:m z (sqt:m (div:m c (sun:m k)))) r] + :: +bernoulli: [r=rng p=@rs] -> [? rng] + :: + :: %.y with probability p, else %.n: draw u~[0,1), return u &((gte:m p .0) (lte:m p .1)) + =^ u r (rs:uni r) + [(lth:m u p) r] + :: +geometric: [r=rng p=@rs] -> [@ud rng] + :: + :: Number of Bernoulli(p) trials up to and including the first success: + :: ceil(ln(u)/ln(1-p)), u drawn from the open (0,1) (avoids ln(0)). + :: p=1 is a special-cased edge that returns 1 directly (ln(1-p) would + :: divide by ln(0) = -inf otherwise). ceil is computed via +toi under + :: a SEPARATE %z-rounding (truncate-toward-zero) door instance, +mz -- + :: +toi respects the door's own rounding mode, and this core's shared + :: +m door is forced %n (round-to-nearest), which would round instead + :: of truncate. Crashes unless 0 < p <= 1. + :: Source + ++ geometric + |= [r=rng:rand p=@rs] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-prob + ?> &((gth:m p .0) (lte:m p .1)) + ?: (equ:m p .1) + [1 r] + =^ u r (rs-oo:uni r) + =/ raw (div:m (log:m u) (log:m (sub:m .1 p))) + =/ fl (abs:si (need (toi:mz raw))) + [?:(=(raw (sun:m fl)) fl +(fl)) r] + :: +poisson: [r=rng lambda=@rs] -> [@ud rng] + :: + :: Poisson(lambda): Knuth's product method for lambda<10 (simple, + :: O(lambda) expected multiplications -- fine at this scale, unusable + :: above it); Hörmann's PTRS (1993) transformed rejection for + :: lambda>=10 (O(1) expected, needed because Knuth's method's cost + :: scales linearly with lambda and would silently become the wrong + :: choice above the threshold rather than crashing -- so this arm + :: switches automatically rather than leaving the choice to the + :: caller). Crashes if lambda <= 0. + :: Source + ++ poisson + |= [r=rng:rand lambda=@rs] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-rate + ?> (gth:m lambda .0) + ?: (lth:m lambda .10) + (poisson-knuth r lambda) + (poisson-ptrs r lambda) + :: +poisson-knuth: k=0, p=1; loop: k+=1, p*=uniform(0,1); accept + :: (return k-1, i.e. the count BEFORE the last increment) once + :: p <= exp(-lambda). Variable-consumption (rejection-free here, but + :: a variable number of multiplications -- expected lambda of them). + ++ poisson-knuth + |= [r=rng:rand lambda=@rs] + ^- [out=@ud r=rng:rand] + =/ bigl (exp:m (neg:m lambda)) + =/ k 0 + =/ p .1 + |- ^- [out=@ud r=rng:rand] + =^ u r (rs:uni r) + =/ p2 (mul:m p u) + ?: (lte:m p2 bigl) + [k r] + $(k +(k), p p2) + :: +poisson-ptrs: Hörmann 1993's transformed rejection with squeeze, + :: verified against NumPy's random_poisson_ptrs (src/distributions/ + :: distributions.c) -- an independent re-implementation of the exact + :: same reference this spec cites, not reconstructed from memory. + :: Draws a candidate k from a scaled/shifted uniform (the "transformed" + :: part), accepts immediately if a cheap squeeze test passes (avoiding + :: the exact log-probability computation on the common path), and + :: falls back to the exact accept/reject test (via +loggam, log(k!)) + :: otherwise. Rejection loop, astronomically bounded in practice. + ++ poisson-ptrs + |= [r=rng:rand lambda=@rs] + ^- [out=@ud r=rng:rand] + =/ slam (sqt:m lambda) + =/ loglam (log:m lambda) + =/ b (add:m .0.931 (mul:m .2.53 slam)) + =/ a (add:m .-0.059 (mul:m .0.02483 b)) + =/ invalpha (add:m .1.1239 (div:m .1.1328 (sub:m b .3.4))) + =/ vr (sub:m .0.9277 (div:m .3.6224 (sub:m b .2))) + |- ^- [out=@ud r=rng:rand] + =^ u1 r (rs:uni r) + =^ v r (rs:uni r) + =/ u (sub:m u1 .0.5) + =/ us (sub:m .0.5 (abs:m u)) + =/ kk (dfloor (add:m (add:m (mul:m (add:m (div:m (mul:m .2 a) us) b) u) lambda) .0.43)) + ?: &((gte:m us .0.07) (lte:m v vr)) + [(abs:si kk) r] + ?: |(=(-1 (cmp:si kk --0)) &((lth:m us .0.013) (gth:m v us))) + $ + =/ lhs (sub:m (add:m (log:m v) (log:m invalpha)) (log:m (add:m (div:m a (mul:m us us)) b))) + =/ rhs (sub:m (add:m (mul:m .-1 lambda) (mul:m (san:m kk) loglam)) (loggam (add:m (san:m kk) .1))) + ?: (lte:m lhs rhs) + [(abs:si kk) r] + $ + :: +dfloor: @rs -> @s, floor (round toward -infinity, NOT toward + :: zero) -- +poisson-ptrs's candidate k can be negative, where + :: truncate-toward-zero (what +toi under +mz gives) would round the + :: wrong way. Positive/exact values: truncation IS the floor. + :: Negative non-integers: truncation rounds up, so subtract 1. + ++ dfloor + |= x=@rs + ^- @s + =/ t (need (toi:mz x)) + ?: |((gte:m x .0) (equ:m x (san:m t))) + t + (dif:si t --1) + :: +loggam: @rs -> @rs, natural log of the gamma function (so log(k!) + :: = loggam(k+1)). Stirling's series with a small-argument recurrence + :: for x<7 (shifts x up by whole integers until >=7, where the series + :: is accurate, then subtracts the shifted-away log-factors back out) + :: -- transcribed from NumPy's random_loggam (src/distributions/ + :: distributions.c), itself from Zhang & Jin's SPECFUN algorithm, not + :: derived independently: log-gamma accuracy is exactly the kind of + :: detail worth matching a vetted reference bit-for-bit rather than + :: reconstructing from the general shape of Stirling's series. + ++ loggam + |= x=@rs + ^- @rs + ?: |((equ:m x .1) (equ:m x .2)) + .0 + =/ a ^- (list @rs) + :~ .0.08333333333333333 .-0.002777777777777778 .0.0007936507936507937 + .-0.0005952380952380952 .0.0008417508417508417 .-0.001917526917526918 + .0.006410256410256411 .-0.02955065359477124 .0.1796443723688307 + .-1.39243221690590 + == + =/ lg2pi .1.8378770664093453 + =/ n ?:((lth:m x .7) (abs:si (need (toi:mz (sub:m .7 x)))) 0) + =/ x0 (add:m x (sun:m n)) + =/ x2 (div:m .1 (mul:m x0 x0)) + =/ gl0 + =/ acc (snag 9 a) + =/ i 0 + |- ^- @rs + ?: =(i 9) + acc + $(acc (add:m (mul:m acc x2) (snag (sub 8 i) a)), i +(i)) + =/ gl + %+ sub:m + %+ add:m + (add:m (div:m gl0 x0) (mul:m .0.5 lg2pi)) + (mul:m (sub:m x0 .0.5) (log:m x0)) + x0 + ?. (lth:m x .7) + gl + =/ fin + =/ gl2 gl + =/ x02 x0 + =/ k 1 + |- ^- @rs + ?: (gth k n) + gl2 + $(gl2 (sub:m gl2 (log:m (sub:m x02 .1))), x02 (sub:m x02 .1), k +(k)) + fin + :: +binomial: [r=rng n=@ p=@rs] -> [@ud rng] + :: + :: Binomial(n,p) via inversion by recurrence: P(X=0)=(1-p)^n, then + :: P(X=k+1) = P(X=k) * (n-k)/(k+1) * p/(1-p), accumulating the CDF + :: until it exceeds a single uniform draw. Only valid for + :: n*min(p,1-p) < 30 (Kachitvichyanukul & Schmeiser 1988) -- above + :: that, inversion needs too many terms and BTPE is the right + :: algorithm, deferred to NEXT-STEPS; this arm crashes rather than + :: running slow (many terms) or, worse, silently truncating. + :: Source + ++ binomial + |= [r=rng:rand n=@ p=@rs] + ^- [out=@ud r=rng:rand] + ~| %rand-bad-prob + ?> &((gte:m p .0) (lte:m p .1)) + =/ q (sub:m .1 p) + =/ qmin ?:((lth:m p .0.5) p q) + ~| %rand-binomial-n-too-large + ?> (lth:m (mul:m (sun:m n) qmin) .30) + =^ u r (rs:uni r) + =/ f (pow:m q (sun:m n)) + =/ cdf f + =/ k 0 + |- ^- [out=@ud r=rng:rand] + ?: (lth:m u cdf) + [k r] + =/ f2 %+ mul:m f + (mul:m (div:m (sun:m (sub n k)) (sun:m +(k))) (div:m p q)) + $(f f2, cdf (add:m cdf f2), k +(k)) + :: +dirichlet: [r=rng alphas=(list @rs)] -> [(list @rs) rng] + :: + :: Dirichlet(alphas): draw one gamma(alpha_i) per entry, normalize by + :: their sum. Crashes on an empty .alphas (and, via +gamma's own + :: precondition, on any non-positive alpha). + :: Source + ++ dirichlet + |= [r=rng:rand alphas=(list @rs)] + ^- [out=(list @rs) r=rng:rand] + ~| %rand-empty-list + ?> !=(~ alphas) + =/ gd + |- ^- [(list @rs) rng:rand] + ?~ alphas + [~ r] + =^ x r (gamma r i.alphas) + =^ rest r $(alphas t.alphas) + [[x rest] r] + =/ total (roll -.gd add:m) + [(turn -.gd |=(g=@rs (div:m g total))) +.gd] + -- +-- +-- diff --git a/librand/desk/lib/rand.hoon b/librand/desk/lib/rand.hoon new file mode 100644 index 0000000..92e0714 --- /dev/null +++ b/librand/desk/lib/rand.hoon @@ -0,0 +1,770 @@ +/- rand +/+ math +=+ rand +:::: /lib/rand -- deterministic pseudo-random number generation (plumbing) +:: +:: NOT CRYPTOGRAPHIC. This library produces deterministic pseudo-randomness +:: for numerical work (Monte Carlo, ML init, shuffles, simulation). It must +:: NOT be used for key material, nonces, or anything adversarial. Entropy +:: acquisition is Arvo's job (`eny`); cryptographic randomness is Zuse's job. +:: +:: This is the FOUNDATION every other rand library imports: the three +:: counter/sequential engines (++philox, ++split-mix, ++pcg), seeding +:: (++seed), the generic engine-dispatched draw/fork (+step, +fork), the +:: door facade (++gen), integer-only uniform deviates (++uni: +bits, +:: +below, +between -- no floats here), and generic sampling (++sample: +:: shuffle/permutation/choice/reservoir -- +alias moved out, see below). +:: +$rng and its component engine-state types live in /sur/rand, not +:: here, so any file that only needs the type (not this library's logic) +:: can get it via a lightweight `/-` import instead of pulling in this +:: whole file. +:: +:: IEEE-754 float generation and the distributions built on it (++dist) +:: live in a separate "porcelain" library, /lib/i754rand -- matching how +:: /lib/twoc, /lib/fixed, /lib/complex, and /lib/unum are already kept +:: separate from /lib/math in this codebase. ++sample's own +alias +:: (Vose's alias method) moved there too, alongside ++dist's +categorical +:: which uses it: +draw needs an actual uniform @rd float draw to decide +:: accept-vs-redirect, so it is NOT float-free the way the rest of +:: ++sample is, and belongs with the porcelain, not the plumbing. +:: Non-float output adapters (twocrand, fixedrand, complexrand, +:: unumrand) are i754rand's siblings, each importing this file for the +:: same engine/uni/sample primitives. +:: +:: Status: milestones 1-6 of rand-spec.md, split into this file plus +:: i754rand per the above. See rand-spec.md and each library's own +:: NEXT-STEPS.md for the full milestone order. +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:::: ++philox :: (3.1) Philox4x32-10 +:: +:: Reference: Salmon, Manuson, Jung, Shaw, "Parallel Random Numbers: As Easy +:: as 1, 2, 3" (SC'11); the Random123 library is the reference +:: implementation. The primary generator: output is a pure function of +:: (key, counter), so a Lagoon jet can fill a ray in parallel and still +:: match this sequential Hoon loop exactly (rand-spec.md section 8). +:: +++ philox + |% + :: +block: [key=@ ctr=@] -> @ + :: + :: The pure Philox4x32-10 function. .ctr packs four 32-bit + :: little-endian lanes c0/c1/c2/c3 (c0 = bits [0,32), ... c3 = bits + :: [96,128)); .key packs two 32-bit lanes k0/k1 (k0 = bits [0,32), k1 = + :: bits [32,64)). Ten rounds of the Philox S-box (mulhilo32 is exact: + :: a plain multiply, then hi/lo = high/low 32 bits of the product): + :: + :: hi0,lo0 = mulhilo32(M0, c0) ; hi1,lo1 = mulhilo32(M1, c2) + :: c0' = hi1 ^ c1 ^ k0 ; c1' = lo1 + :: c2' = hi0 ^ c3 ^ k1 ; c3' = lo0 + :: k0 += W0 ; k1 += W1 (mod 2^32) + :: + :: Output is the four final lanes repacked as one 128-bit atom in the + :: same little-endian lane order as .ctr (c0' lowest). KAT vectors + :: (Random123's kat_vectors file -- all-zero, all-0xffffffff, and the + :: pi-digits vector -- cross-checked against an independent Python + :: re-implementation of this exact algorithm) are in tests/lib/rand.hoon. + :: Examples + :: > (block:philox 0 0) + :: 0x9b00.dbd8.bc57.ac4c.e169.c58d.6627.e8d5 + :: Source + ++ block + ~/ %block + |= [key=@ ctr=@] + ^- @ + =/ c0 (cut 0 [0 32] ctr) + =/ c1 (cut 0 [32 32] ctr) + =/ c2 (cut 0 [64 32] ctr) + =/ c3 (cut 0 [96 32] ctr) + =/ k0 (cut 0 [0 32] key) + =/ k1 (cut 0 [32 32] key) + =/ n 0 + |- ^- @ + ?: =(n 10) + :(add c0 (lsh [0 32] c1) (lsh [0 64] c2) (lsh [0 96] c3)) + =/ m0 (mul 0xd251.1f53 c0) + =/ hi0 (rsh [0 32] m0) + =/ lo0 (end [0 32] m0) + =/ m1 (mul 0xcd9e.8d57 c2) + =/ hi1 (rsh [0 32] m1) + =/ lo1 (end [0 32] m1) + %= $ + c0 (mix hi1 (mix c1 k0)) + c1 lo1 + c2 (mix hi0 (mix c3 k1)) + c3 lo0 + k0 (end [0 32] (add k0 0x9e37.79b9)) + k1 (end [0 32] (add k1 0xbb67.ae85)) + n +(n) + == + :: +next: phil -> [out=@ p=phil] + :: + :: One Philox4x32-10 draw: out = (block key ctr), new state advances the + :: counter by 1 (mod 2^128), key unchanged. See +step (below) for how + :: callers reduce the 128-bit .out to a 64-bit draw. + :: Examples + :: > (next:philox [key=0 ctr=0]) + :: [out=0x9b00.dbd8.bc57.ac4c.e169.c58d.6627.e8d5 p=[key=0 ctr=1]] + :: Source + ++ next + ~/ %next + |= p=phil + ^- [out=@ p=phil] + [(block key.p ctr.p) [key.p (mod +(ctr.p) (bex 128))]] + -- +:: +:::: ++split-mix :: (3.2) SplitMix64 +:: +:: Reference: Steele, Lea, Flood (OOPSLA'14); Vigna's public-domain C. The +:: primary generator (++philox, milestone 2) is counter-based; SplitMix64 +:: here serves two roles: a small sequential generator in its own right +:: (%sm64 in +$rng), and the seed-mixing primitive every engine's +seed +:: path and (later) +fork depend on -- see +finalize below. +:: +++ split-mix + |% + :: +finalize: @ -> @ + :: + :: The z-mixing tail of SplitMix64's step function, exposed as its own + :: arm because +seed's +mix (below) reuses it on arbitrary atoms that are + :: NOT a sequential-generator state (no golden-ratio increment applies to + :: those inputs -- only +next's state-advance gets that step). Maps + :: 0 -> 0 (a fixed point of the finalizer alone, without the increment); + :: this is expected and is why +next always adds the increment before + :: finalizing, and why +mix's degenerate mix(0,0) = 0 is documented + :: rather than "fixed" (see +mix). + :: Examples + :: > (finalize:split-mix 0x9e37.79b9.7f4a.7c15) + :: 0xe220.a839.7b1d.cdaf + :: Source + ++ finalize + |= z=@ + ^- @ + =. z (end [0 64] (mul (mix z (rsh [0 30] z)) 0xbf58.476d.1ce4.e5b9)) + =. z (end [0 64] (mul (mix z (rsh [0 27] z)) 0x94d0.49bb.1331.11eb)) + (mix z (rsh [0 31] z)) + :: +next: sm64 -> [out=@ s=sm64] + :: + :: One SplitMix64 draw: advance state by the golden-ratio increment + :: (mod 2^64), then finalize. Matches Vigna's reference C bit-for-bit; + :: see the KAT vectors in tests/lib/rand.hoon. + :: Examples + :: > (next:split-mix 0) + :: [out=0xe220.a839.7b1d.cdaf s=0x9e37.79b9.7f4a.7c15] + :: Source + ++ next + ~/ %next + |= s=sm64 + ^- [out=@ s=sm64] + =/ s2 (end [0 64] (add s 0x9e37.79b9.7f4a.7c15)) + [(finalize s2) s2] + :: +split: sm64 -> [a=@ b=@ s=sm64] + :: + :: Two decorrelated child states from one parent state, via two +next + :: draws used directly as child seeds. Gamma-stepping per the original + :: paper is overkill for this library's needs and adds a second tunable + :: parameter for no benefit here; two draws is adequate and simpler. + :: Examples + :: > (split:split-mix 0) + :: [a=0xe220.a839.7b1d.cdaf b=0x6e78.9e6a.a1b9.65f4 s=0x3c6e.f372.fe94.f82a] + :: Source + ++ split + |= s=sm64 + ^- [a=@ b=@ s=sm64] + =^ a s (next s) + =^ b s (next s) + [a b s] + -- +:: +:::: ++seed :: (4) seeding +:: +:: No implicit global seeding. There is no "default rng"; callers own +:: their state. NOT CRYPTOGRAPHIC -- see the file-level warning above; +:: +from-eny carries the same warning explicitly since it is the arm most +:: likely to be reached for by someone wanting "real" randomness. +:: +++ seed + |% + :: +from-atom: [eng=?(%phil %sm64 %pcg) a=@] -> rng + :: + :: Seed any engine deterministically from a single atom, treating .a as + :: raw SplitMix64 state and drawing as many outputs as that engine's + :: state needs. %sm64 needs no pre-draw: SplitMix64's own state can be + :: seeded with any value (Vigna's reference explicitly allows this), and + :: +next's golden-ratio increment does the mixing on first use -- drawing + :: an extra output here would just be redundant hashing. %phil takes the + :: first output as .key with .ctr reset to 0. %pcg takes two outputs as + :: .state and .inc, forcing .inc odd (PCG64 requires an odd increment). + :: Examples + :: > (from-atom:seed %sm64 0) + :: [%sm64 s=0] + :: > (from-atom:seed %phil 0) + :: [%phil p=[key=0xe220.a839.7b1d.cdaf ctr=0]] + :: > (from-atom:seed %pcg 0) + :: [%pcg p=[state=0xe220.a839.7b1d.cdaf inc=0x6e78.9e6a.a1b9.65f5]] + :: Source + ++ from-atom + |= [eng=?(%phil %sm64 %pcg) a=@] + ^- rng + ?- eng + %sm64 [%sm64 s=(end [0 64] a)] + %phil + =^ key a (next:split-mix (end [0 64] a)) + [%phil p=[key=key ctr=0]] + %pcg + =^ st a (next:split-mix (end [0 64] a)) + =^ inc a (next:split-mix a) + [%pcg p=[state=st inc=(con inc 1)]] + == + :: +from-eny: [eng=?(%phil %sm64 %pcg) eny=@uvJ] -> rng + :: + :: Same pipeline as +from-atom, over Arvo entropy. NOT CRYPTOGRAPHIC -- + :: see the file-level warning. `eny` is 512 bits; +fold-wide below + :: compresses it to the 64 bits +from-atom expects. + :: Source + ++ from-eny + |= [eng=?(%phil %sm64 %pcg) eny=@uvJ] + ^- rng + (from-atom eng (fold-wide eny)) + :: +fold-wide: @ -> @ + :: + :: Compress an atom of ANY width to 64 bits: split into 64-bit words + :: (low word first, via +met bloq 6 for the word count) and fold them + :: through +mix in order (acc starts at 0; acc = (mix acc word) per + :: word). Not itself an rng seed -- feeds +from-atom (for @uvJ eny, + :: always 8 words) and +fork (for salts wider than 64 bits, section + :: 2.1c) alike, so there is exactly one fold algorithm in the library. + :: Examples + :: > (fold-wide:seed 0) + :: 0 + :: Source + ++ fold-wide + |= a=@ + ^- @ + =/ n (met 6 a) + =/ i 0 + =/ acc 0 + |- ^- @ + ?: =(i n) + acc + $(acc (mix acc (cut 6 [i 1] a)), i +(i)) + :: +mix: [a=@ b=@] -> @ + :: + :: Combine two 64-bit seed words (e.g. a ship + a per-agent salt) into + :: one. Fully specified as a two-word compression, NOT a single + :: finalizer call on the raw 128-bit concatenation of .a and .b -- the + :: finalizer is a function of one 64-bit word, so a wider input must be + :: folded, not fed in directly: + :: + :: w0 = a mod 2^64 ; w1 = b mod 2^64 + :: return finalize(finalize(w0) XOR w1) + :: + :: Order-sensitive by construction: (mix a b) differs from (mix b a) in + :: general, since .a is absorbed through one extra +finalize application + :: relative to .b. This is the primitive +fork (milestone 2) uses to + :: make (fork (fork r a) b) differ from (fork (fork r b) a): each nested + :: +fork call adds another +finalize layer at a different depth. + :: Deterministic, no rng threading (pure @ -> @). + :: + :: Degenerate case, documented rather than "fixed": mix(0,0) = 0, because + :: +finalize maps 0 -> 0 when no golden-ratio increment has been mixed in + :: first (see +finalize). Real seeds and salts are essentially never + :: both exactly zero, so this is not a practical concern. + :: + :: NAMING FOOTGUN: this arm's name collides with the Hoon standard + :: library's `++mix` (bitwise XOR). Inside `++seed`, an unqualified + :: `(mix a b)` call resolves to THIS arm, not stdlib XOR -- which is why + :: the body below reaches for `^mix` to get bitwise XOR. Any other code + :: added to `++seed` that wants stdlib XOR must do the same; callers + :: outside `++seed` are unaffected. + :: Examples + :: > (mix:seed 1 2) + :: 0xef30.b01c.2974.aeeb + :: > (mix:seed 2 1) + :: 0x3ec2.d42f.3a45.cc6e + :: Source + ++ mix + |= [a=@ b=@] + ^- @ + =/ w0 (end [0 64] a) + =/ w1 (end [0 64] b) + (finalize:split-mix (^mix (finalize:split-mix w0) w1)) + -- +:: +:::: +step :: (2) generic draw +:: +:: One 64-bit draw, engine-dispatched. Lives at the top level (not nested +:: under any one engine core) because it dispatches across the +$rng +:: tagged union -- this is what lets ++uni/++dist/++sample (later +:: milestones) be written once each, not once per engine. +:: +++ step + |= r=rng + ^- [out=@ r=rng] + ?- -.r + %sm64 + =^ out s.r (next:split-mix s.r) + [out r] + :: + %phil + :: +next:philox returns the full 128-bit block; +step keeps bits + :: [0,64) -- lanes c0'/c1', the low two lanes of the little-endian + :: repack in +block:philox. The top 64 bits (c2'/c3') are discarded; + :: wasting them is acceptable at v1 (rand-spec.md section 2), and a + :: buffered variant that reuses both halves is a NEXT-STEPS item. + =^ blk p.r (next:philox p.r) + [(end [0 64] blk) r] + :: + %pcg + =^ out p.r (next:pcg p.r) + [out r] + == +:: +:::: +fork :: (2.1c) key derivation +:: +:: Deterministic, independent child stream from a parent stream and a +:: salt -- the JAX-style key-derivation payoff of counter-based-first +:: (rand-spec.md section 2.1c). Callers fork by path instead of +:: threading one sequential state through a computation tree: draws are +:: independent of sibling evaluation order, and a draw inserted in one +:: branch doesn't perturb any other branch. +:: Examples +:: > =/ a (fork (from-atom:seed %phil 0) 1) +:: > =/ b (fork (from-atom:seed %phil 0) 2) +:: > =(a b) +:: %.n +:: Source +++ fork + |= [r=rng salt=@] + ^- rng + :: Salts up to 64 bits feed +mix directly (the common case: small + :: integer salts, e.g. per-event or per-layer indices). Wider salts + :: (e.g. a caller-combined `(cat 6 a b)`, or an @uvJ) are folded to 64 + :: bits first via +fold-wide, so no salt bits are silently dropped -- + :: +mix itself only ever consumes one 64-bit word per side. + =/ sw ?:((lte (met 6 salt) 1) salt (fold-wide:seed salt)) + ?- -.r + %phil r(key.p (mix:seed key.p.r sw), ctr.p 0) + :: + :: "child increment derived the same way [as %phil's key], forced + :: odd; state re-mixed" (section 2.1c) -- read here as: both fields + :: pushed through +mix with the same salt-derived word. + %pcg + r(state.p (mix:seed state.p.r sw), inc.p (con (mix:seed inc.p.r sw) 1)) + :: + :: "the SplitMix split construction" (section 2.1c) is read here as + :: reusing the same +mix primitive +split's decorrelation relies on, + :: salt-directed rather than sequential -- NOT a literal call to + :: +split (which takes no salt and can't be path-directed). This is + :: this implementation's resolution of that spec ambiguity. + %sm64 r(s (mix:seed s.r sw)) + == +:: +:::: ++gen :: (2.1b) door facade +:: +:: Ergonomic sugar over the functional core, never a storage format. +:: HARD RULE (section 2.1b): persisting this door in agent state pins a +:: battery across library upgrades (the classic +og misuse) -- persist +:: the +$rng noun, reconstruct the door locally. The facade adds no +:: capability; it is a strict wrapper over +step/+fork, so jets register +:: on those functional arms only. +:: +++ gen + |_ r=rng + :: +draw: gen -> [@ _..draw] + :: + :: Door wrapper over +step. `..draw` (not `+>`/`+>.$`) is the reliable + :: way to reference "this door, before any local rebinding" from an + :: arm's body regardless of whether the arm itself has a `|=` sample + :: (`+>`/`+>.$` axis arithmetic differs depending on that, which is a + :: footgun in itself -- `..` sidesteps it). + :: Examples + :: > =/ g ~(. gen (from-atom:seed:rand %sm64 0)) + :: > =^ x g draw:g + :: > x + :: 16.294.208.416.658.607.535 + :: Source + ++ draw + ^- [@ _..draw] + =^ out r (step r) + [out ..draw(r r)] + :: +fork: [gen @] -> _..fork + :: + :: Door wrapper over the functional +fork. NAMING FOOTGUN (same class + :: as ++seed's +mix, above): this arm's name shadows the top-level + :: +fork gate, so calling it unqualified from inside this arm's own + :: body would recurse into itself with the wrong arity -- the body + :: below reaches for `^fork` to reach the top-level gate. + :: Source + ++ fork + |= salt=@ + ^- _..fork + ..fork(r (^fork r salt)) + -- +:: +:: +:::: ++uni :: (5) uniform deviates (integer) +:: +:: +bits/+below/+between only -- the integer/bit-pattern primitives every +:: adapter (twocrand, fixedrand, unumrand) and i754rand's own float +:: generators are built on. IEEE-754 float generation (+rs/+rd/+rh/+rq/ +:: +rs-oo/+rd-oo) lives in i754rand's OWN ++uni instead of here, kept +:: separate so non-float adapters never need to import /lib/math. +:: +++ uni + |% + :: +bits: [r=rng n=@] -> [@ rng] + :: + :: .n uniform bits, assembled from ceil(n/64) +step draws, little-endian + :: (the first draw is bits [0,64) of the result, the second is bits + :: [64,128), and so on); the LAST draw is masked down to n's remainder + :: mod 64 when that remainder is nonzero (i.e. when n isn't itself a + :: multiple of 64). + :: Examples + :: > (bits:uni:rand (from-atom:seed:rand %sm64 0) 4) + :: [out=15 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] + :: Source + ++ bits + |= [r=rng n=@] + ^- [out=@ r=rng] + =/ full (div n 64) + =/ rem (mod n 64) + =/ words ?:(=(rem 0) full +(full)) + =/ i 0 + =/ acc 0 + |- ^- [out=@ r=rng] + ?: =(i words) + [acc r] + =^ out r (step r) + =/ bc ?:(&(=(i (dec words)) !=(rem 0)) rem 64) + %= $ + acc (add acc (lsh [0 (mul i 64)] (end [0 bc] out))) + i +(i) + == + :: +below: [r=rng n=@] -> [@ rng] + :: + :: Uniform in [0,n), UNBIASED, via Lemire's multiply-shift with + :: rejection (Lemire 2019, "Fast Random Integer Generation in an + :: Interval"). .t depends only on .n, so it's computed once outside the + :: rejection loop, not per draw. The primitive everything else + :: (shuffle, choice, +between) uses -- never ship modulo-bias + :: `(mod x n)` instead, here or anywhere else in this library. + :: Crashes (`?>`) on n=0. + :: Examples + :: > (below:uni:rand (from-atom:seed:rand %sm64 0) 10) + :: [out=8 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] + :: Source + ++ below + ~/ %below + |= [r=rng n=@] + ^- [out=@ r=rng] + ~| %rand-zero-modulus + ?> !=(n 0) + =/ t (mod (sub (bex 64) n) n) + |- ^- [out=@ r=rng] + =^ x r (step r) + =/ m (mul x n) + =/ l (end [0 64] m) + ?: &((lth l n) (lth l t)) + $ + [(rsh [0 64] m) r] + :: +between: [r=rng a=@s b=@s] -> [@s rng] + :: + :: Inclusive signed range [a,b], via +below on the span (b - a + 1), + :: offset back by .a. Crashes (`?>`) if a > b. + :: Examples + :: > (between:uni:rand (from-atom:seed:rand %sm64 0) -5 --5) + :: [out=--4 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] + :: Source + ++ between + |= [r=rng a=@s b=@s] + ^- [out=@s r=rng] + ~| %rand-bad-range + ?> !=(--1 (cmp:si a b)) + =/ span +((abs:si (dif:si b a))) + =^ off r (below r span) + [(sum:si a (sun:si off)) r] + -- +:::: ++pcg :: (3.3) PCG64 XSL-RR +:: +:: Reference: O'Neill 2014, pcg-random.org; the pcg-c library +:: (`include/pcg_variants.h`) is the reference implementation. +:: +:: **Spec correction, verified against pcg-c's `pcg_setseq_128_xsl_rr_64_ +:: random_r` and independently against NumPy's vendored `pcg64.orig.h` +:: (both consistent): the real reference ADVANCES the state first, THEN +:: computes the output permutation from the NEW state** -- not "output +:: the current state, then advance" as an earlier draft of rand-spec.md +:: claimed. That earlier claim was backwards; see rand-spec.md section +:: 3.3 for the correction. This is the implementation of the corrected +:: order: +next below steps first, then outputs. +:: +++ pcg + |% + :: +step-lcg: pcg64 -> pcg64 + :: + :: The bare LCG advance, state' = state * MULT + inc (mod 2^128), no + :: output. MULT is PCG's fixed 128-bit default multiplier; .inc is the + :: caller's odd stream increment (enforced odd at seed/fork time, not + :: here). + :: Source + ++ step-lcg + |= p=pcg64 + ^- pcg64 + %= p + state (end [0 128] (add (mul state.p 0x2360.ed05.1fc6.5da4.4385.df64.9fcc.f645) inc.p)) + == + :: +rotr64: [value=@ rot=@] -> @ + :: + :: 64-bit right-rotate: (value >> rot) | (value << (64-rot)), matching + :: pcg_rotr_64. .rot is always < 64 in practice (it's the top 6 bits of + :: a 128-bit state, so already in [0,63]); the formula is correct at + :: rot=0 too (the left-shift-by-64 term vanishes under the 64-bit mask, + :: leaving value unchanged) so no special case is needed. + :: Source + ++ rotr64 + |= [value=@ rot=@] + ^- @ + (mix (rsh [0 rot] value) (end [0 64] (lsh [0 (sub 64 rot)] value))) + :: +output: @ -> @ + :: + :: The xsl-rr output permutation of a 128-bit state to a 64-bit output: + :: rot = state >> 122; xored = (state >> 64) ^ (state & mask64); output + :: = rotr64(xored, rot). Matches pcg_output_xsl_rr_128_64. + :: Source + ++ output + |= state=@ + ^- @ + =/ hi (rsh [0 64] state) + =/ lo (end [0 64] state) + =/ rot (rsh [0 122] state) + (rotr64 (mix hi lo) rot) + :: +next: pcg64 -> [out=@ p=pcg64] + :: + :: One PCG64 draw. Per the spec correction above: advance first + :: (+step-lcg), THEN compute the output from the resulting NEW state -- + :: matching `pcg_setseq_128_xsl_rr_64_random_r`'s actual order, not the + :: order an earlier draft of this library's spec claimed. + :: Examples + :: > (next:pcg [state=0xde2b.ce05.be01.3be3.d3f6.c45a.41e5.4320 inc=0x6d]) + :: [out=0x86b1.da1d.7206.2b68 p=[state=0x10af.065f.4ea9.6e85.7bb2.a788.6ecb.d80d inc=0x6d]] + :: Source + ++ next + ~/ %next + |= p=pcg64 + ^- [out=@ p=pcg64] + =/ p2 (step-lcg p) + [(output state.p2) p2] + :: +advance: [p=pcg64 delta=@] -> pcg64 + :: + :: Advance .p by .delta LCG steps in O(log delta) via the standard PCG + :: skip-ahead trick: the affine map step(s) = s*mult + inc composes + :: under repeated squaring, so .delta steps can be applied as one + :: affine map (acc-mult, acc-plus) built by binary exponentiation, + :: rather than .delta sequential steps. +jump is the .delta = 2^64 + :: case this library actually needs (for stream partitioning); +advance + :: is exposed separately so the algorithm can be tested at a small, + :: tractable .delta instead of only at the astronomical one +jump uses: + :: `(advance p 3)` must equal three sequential `next:pcg` steps (see + :: tests/lib/rand.hoon). + :: Source + ++ advance + |= [p=pcg64 delta=@] + ^- pcg64 + =/ cur-mult 0x2360.ed05.1fc6.5da4.4385.df64.9fcc.f645 + =/ cur-plus inc.p + =/ acc-mult 1 + =/ acc-plus 0 + |- ^- pcg64 + ?: =(delta 0) + %= p + state (end [0 128] (add (mul acc-mult state.p) acc-plus)) + == + ?: =(1 (dis delta 1)) + %= $ + acc-mult (end [0 128] (mul acc-mult cur-mult)) + acc-plus (end [0 128] (add (mul acc-plus cur-mult) cur-plus)) + cur-plus (end [0 128] (mul (add cur-mult 1) cur-plus)) + cur-mult (end [0 128] (mul cur-mult cur-mult)) + delta (rsh [0 1] delta) + == + %= $ + cur-plus (end [0 128] (mul (add cur-mult 1) cur-plus)) + cur-mult (end [0 128] (mul cur-mult cur-mult)) + delta (rsh [0 1] delta) + == + :: +jump: pcg64 -> pcg64 + :: + :: Advance by exactly 2^64 steps, for stream partitioning (give + :: different logical streams non-overlapping windows of the same + :: underlying sequence). A thin +advance call at the one .delta this + :: library needs. + :: Source + ++ jump + |= p=pcg64 + ^- pcg64 + (advance p (bex 64)) + -- +:: +:::: ++sample :: (7) sampling +:: +:: Shuffles, permutations, choice, reservoir sampling. +below:uni is the +:: shared unbiased primitive throughout, per section 5.1's policy: never +:: modulo-bias. Vose's alias method (section 6.1) moved to i754rand: its +:: +draw needs an actual uniform @rd float draw (to decide accept vs. +:: redirect), so it isn't float-free the way the rest of this core is. +:: +++ sample + |% + :: +shuffle: [r=rng l=(list)] -> [(list) rng] + :: + :: Fisher-Yates, iterating from the end down: for i from n-1 to 1, swap + :: position i with a uniform position in [0,i] (+below:uni, inclusive + :: via i+1). Implemented over a (map @ud _elem) rather than repeated + :: +snag/+oust on the list itself, which would be O(n^2); the map + :: round-trip is O(n log n). Wet gate so the element type is + :: preserved for the caller (a dry `(list)` would erase it to `*`). + :: Examples + :: > (shuffle:sample (from-atom:seed:rand %sm64 0) ~[1 2 3 4 5]) + :: [~[3 4 1 2 5] [%sm64 s=0x78dd.e6e5.fd29.f054]] + :: Source + ++ shuffle + |* [r=rng l=(list)] + ^+ [l r] + =/ n (lent l) + ?: (lth n 2) + [l r] + =/ elem ?>(?=(^ l) i.l) + =/ m (~(gas by *(map @ud _elem)) (turn (gulf 0 (dec n)) |=(k=@ud [k (snag k l)]))) + =/ rr r + =/ ix (dec n) + |- ^+ [l r] + ?: =(ix 0) + [(turn (gulf 0 (dec n)) |=(k=@ud (~(got by m) k))) rr] + =^ j rr (below:uni rr +(ix)) + =/ vi (~(got by m) ix) + =/ vj (~(got by m) j) + %= $ + m (~(put by (~(put by m) ix vj)) j vi) + ix (dec ix) + == + :: +permutation: [r=rng n=@] -> [(list @ud) rng] + :: + :: A uniform random permutation of 0..n-1: +shuffle of (gulf 0 (dec n)). + :: Source + ++ permutation + |= [r=rng n=@] + ^- [(list @ud) rng] + ?: =(n 0) + [~ r] + (shuffle r (gulf 0 (dec n))) + :: +choice: [r=rng l=(list)] -> [* rng] + :: + :: One uniformly-chosen element. Wet gate (rand-spec.md section 7 + :: writes this arm's signature dry, `[* rng]` -- but a bare `*` return + :: loses the element's type at every call site for no benefit, and a + :: dry `(list)` argument runs into a real Hoon type-inference wall at + :: `+snag` -- mull-grow/nest-fail trying to prove a `(list *)` is + :: non-null after the `?~` guard. `|*` sidesteps both: each call site + :: gets its own precise element type, same reasoning as +shuffle. + :: Crashes on an empty list. + :: Source + ++ choice + |* [r=rng l=(list)] + ~| %rand-empty-list + ?> !=(~ l) + =^ i r (below:uni r (lent l)) + [(snag i l) r] + :: +choices: [r=rng n=@ l=(list)] -> [(list) rng] + :: + :: .n elements chosen uniformly WITH replacement (independent +choice + :: draws). Crashes on an empty .l if n > 0. + :: Source + ++ choices + |* [r=rng n=@ l=(list)] + ^+ [l r] + ?: =(n 0) + [~ r] + ~| %rand-empty-list + =/ elem ?>(?=(^ l) i.l) + =/ i 0 + =/ acc *(list _elem) + =/ rr r + |- ^+ [l r] + ?: =(i n) + [(flop acc) rr] + =^ x rr (choice rr l) + %= $ + i +(i) + acc [x acc] + == + :: +sample-n: [r=rng k=@ l=(list)] -> [(list) rng] + :: + :: .k elements WITHOUT replacement, via partial Fisher-Yates (the first + :: -- here, for implementation convenience, the LAST -- k positions of + :: a full shuffle): iterate i from n-1 down to n-k, swapping position i + :: with a uniform position in [0,i], same as +shuffle but stopping + :: early: this is the spec'd "first n of a permutation," just taken + :: from whichever end +shuffle itself iterates from, and either end is + :: equally uniform. NOT repeated-rejection sampling. Crashes if k > the + :: list's length. Wet gate, same reasoning as +shuffle. + :: Source + ++ sample-n + |* [r=rng k=@ l=(list)] + ^+ [l r] + =/ n (lent l) + ~| %rand-bad-count + ?> (lte k n) + ?: =(k 0) + [~ r] + =/ elem ?>(?=(^ l) i.l) + =/ m (~(gas by *(map @ud _elem)) (turn (gulf 0 (dec n)) |=(j=@ud [j (snag j l)]))) + =/ rr r + =/ ix (dec n) + =/ stop (sub n k) + |- ^+ [l r] + ?: =(ix stop) + [(turn (gulf stop (dec n)) |=(j=@ud (~(got by m) j))) rr] + =^ j rr (below:uni rr +(ix)) + =/ vi (~(got by m) ix) + =/ vj (~(got by m) j) + %= $ + m (~(put by (~(put by m) ix vj)) j vi) + ix (dec ix) + == + :: +reservoir: [r=rng n=@ l=(list)] -> [(list) rng] + :: + :: Algorithm R (Vitter 1985): a uniform sample of .n items from .l, + :: processed one at a time (the algorithm this library's own +below is + :: built on doesn't need true streaming, but the same algorithm serves + :: agents that DO stream events one at a time). The first .n items + :: seed the reservoir; each later item at index i replaces a uniformly + :: chosen reservoir slot with probability n/(i+1) (drawn as "is the + :: uniform draw in [0,i] less than n"). Crashes if n > the list's + :: length. Wet gate, same reasoning as +shuffle. + :: Source + ++ reservoir + |* [r=rng n=@ l=(list)] + =/ len (lent l) + ~| %rand-bad-count + ?> (lte n len) + ^+ [(scag n l) r] + ?: =(n 0) + [~ r] + =/ elem ?>(?=(^ l) i.l) + =/ m (~(gas by *(map @ud _elem)) (turn (gulf 0 (dec n)) |=(k=@ud [k (snag k l)]))) + =/ rr r + =/ ix n + |- ^+ [(scag n l) r] + ?: =(ix len) + [(turn (gulf 0 (dec n)) |=(k=@ud (~(got by m) k))) rr] + =^ j rr (below:uni rr +(ix)) + ?: (lth j n) + %= $ + m (~(put by m) j (snag ix l)) + ix +(ix) + == + $(ix +(ix)) + -- +-- diff --git a/librand/desk/lib/twocrand.hoon b/librand/desk/lib/twocrand.hoon new file mode 100644 index 0000000..b8504d9 --- /dev/null +++ b/librand/desk/lib/twocrand.hoon @@ -0,0 +1,62 @@ +/+ rand, twoc +:::: /lib/twocrand -- two's-complement integer output adapter (rand-spec.md +:::: section 12.2) +:: +:: Routes /lib/rand's engines through /lib/twoc's width-keyed two's- +:: complement representation. No new engine work: both arms are thin +:: compositions of /lib/rand's +bits/+below and /lib/twoc's +twid door. +:: This is /lib/rand's LEANEST adapter -- no floats anywhere, so no +:: dependency on /lib/i754rand or /lib/math, unlike /lib/complexrand. +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:: +twoc-full: [r=rng w=@] -> [@ rng] +:: +:: .w raw bits, interpreted directly as the width-.w two's-complement +:: pattern. This is the IDENTITY on the raw bits -- a two's-complement +:: value IS just a bit pattern, so a uniform draw over raw bits already +:: IS a uniform draw over two's-complement values at that width. No +:: encoding step needed (rand-spec.md section 12.2 marks this "Done"). +:: Examples +:: > (twoc-full:twocrand (from-atom:seed:rand %sm64 0) 8) +:: [out=175 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] +:: Source +++ twoc-full + |= [r=rng:rand w=@] + ^- [out=@ r=rng:rand] + (bits:uni:rand r w) +:: +twoc-between: [r=rng w=@ a=@ b=@] -> [@ rng] +:: +:: Inclusive range [a,b], where .a and .b are width-.w two's-complement +:: PATTERNS (raw bits, not @s) and the range is inclusive in TWOC order +:: (negatives below non-negatives, per /lib/twoc's own +gth/+lth). +:: Crashes unless a<=b in that order. +:: +:: Bias to unsigned: the span b-a+1 is computed via /lib/twoc's own +:: modular +sub, so it wraps correctly regardless of a/b's individual +:: signs; drawn uniformly via Lemire's +below (rand-spec.md section +:: 5.1); biased back via /lib/twoc's modular +add. /lib/twoc is the +:: primitive /lib/fixedrand's +fixed-between delegates to, per +:: rand-spec.md section 12.1. +:: +:: JET WARNING (rand-spec.md section 12.2, verbatim): the span of a +:: full-width-64 range (a=minint, b=maxint) needs 65 bits, not 64 -- +:: Hoon's arbitrary-precision @ doesn't care (this Hoon implementation +:: is correct at every width with no special case), but a future C jet +:: computing the span in a machine uint64_t WILL overflow for that one +:: case. The jet must special-case span=2^w, or compute the span in a +:: wider-than-w integer type. +:: Examples +:: > (twoc-between:twocrand (from-atom:seed:rand %sm64 0) 8 0xfb 0x5) +:: [out=4 r=[%sm64 s=0x9e37.79b9.7f4a.7c15]] +:: Source +++ twoc-between + |= [r=rng:rand w=@ a=@ b=@] + ^- [out=@ r=rng:rand] + =/ tw ~(. twid:twoc w) + ~| %rand-bad-range + ?> !(gth:tw a b) + =/ span +((sub:tw b a)) + =^ off r (below:uni:rand r span) + [(add:tw a off) r] +-- diff --git a/librand/desk/lib/unumrand.hoon b/librand/desk/lib/unumrand.hoon new file mode 100644 index 0000000..99ba524 --- /dev/null +++ b/librand/desk/lib/unumrand.hoon @@ -0,0 +1,135 @@ +/+ rand, unum +:::: /lib/unumrand -- posit output adapter (rand-spec.md section 12.4) +:: +:: Two inequivalent uniform semantics, named so they cannot be confused -- +:: posits are TAPERED, so consecutive bit patterns are NOT evenly spaced in +:: value: +:: +:: +posit-lattice -- uniform over raw BIT PATTERNS excluding NaR. The +:: induced VALUE distribution is only approximately log-uniform (dense +:: near +-1, sparse toward the extremes). A fuzzing/property-testing +:: primitive, not a statistics one. Ships at all five width doors +:: (rpb/rph/rps/rpd/rpq) -- no bit-count subtlety, exact at any width. +:: +:: +posit-unit -- uniform VALUE on [0,1), exact. Ships at posit8/16/32 +:: only (see librand/NEXT-STEPS.md's milestone-7-pass-2 entry for the +:: full k=4n derivation and the scope decision). Draw k raw bits u, +:: encode the dyadic u*2^-k -- the g-layer value [%p %.y -k u] (%z when +:: u=0) -- via /lib/unum's existing +bit (round-to-nearest-even, +:: saturating; confirmed no rounding-mode parameter needed). For this +:: to be EXACTLY the round-to-nearest image of continuous uniform (not +:: approximately), k must exceed the finest rounding-cell width anywhere +:: in [0,1) -- which sits near zero, at exponent -(4(n-2)+1), since +:: minpos = 2^-4(n-2). k=4n (32/64/128 for posit8/16/32) gives a +:: CONSTANT 7-bit safety margin at any width: 4n - (4(n-2)+1) = 7 always. +:: Verified against an exact mpmath-free rational oracle (every quantity +:: here is already dyadic, so plain Fraction arithmetic is exact) in +:: librand/tools/posit_unit_check.py -- see NEXT-STEPS.md for how to run +:: the chi-square check against ship-drawn draws. +:: +:: Distributions ("sample at @rd, convert") need NO new /lib/unum plumbing, +:: unlike fixedrand's +from-rd: /lib/unum already ships +from-rh/rs/rd/rq +:: at every width door. Not shipped as dedicated wrapper arms here, per +:: the same decision as fixedrand/complexrand -- documented composition + +:: one test proving it (tests/lib/unumrand.hoon). +:: +~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/math, /lib/twoc) +|% +:: +pl: shared +posit-lattice implementation, keyed on bloq -- mirrors +:: /lib/unum's own +pp idiom (one generic door, thin per-width forwarders +:: below), so the algorithm exists exactly once regardless of width. +++ pl + |= [r=rng:rand =bloq] + ^- [out=@ r=rng:rand] + =/ w ~(. pp:unum bloq) + |- ^- [out=@ r=rng:rand] + =^ b r (bits:uni:rand r n:w) + ?: =(b nar:w) + $ + [b r] +:: +pu: shared +posit-unit implementation, keyed on bloq and the drawn +:: bit-count k. Only called from rpb/rph/rps below (the widths in scope). +++ pu + |= [r=rng:rand =bloq k=@] + ^- [out=@ r=rng:rand] + =/ w ~(. pp:unum bloq) + =^ u r (bits:uni:rand r k) + :_ r + %- bit:w + ?: =(u 0) + [%z ~] + [%p %.y (dif:si --0 (sun:si k)) u] +:: +rpb: posit8 (n=8) adapters. +:: Source +++ rpb + |% + :: +posit-lattice: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-lattice:rpb:unumrand (from-atom:seed:rand %sm64 0)) + :: 0xaf + :: Source + ++ posit-lattice |=(r=rng:rand (pl r 3)) + :: +posit-unit: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-unit:rpb:unumrand (from-atom:seed:rand %sm64 0)) + :: 0x37 + :: Source + ++ posit-unit |=(r=rng:rand (pu r 3 32)) + -- +:: +rph: posit16 (n=16) adapters. +:: Source +++ rph + |% + :: +posit-lattice: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-lattice:rph:unumrand (from-atom:seed:rand %sm64 0)) + :: 0xcdaf + :: Source + ++ posit-lattice |=(r=rng:rand (pl r 4)) + :: +posit-unit: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-unit:rph:unumrand (from-atom:seed:rand %sm64 0)) + :: 0x3e22 + :: Source + ++ posit-unit |=(r=rng:rand (pu r 4 64)) + -- +:: +rps: posit32 (n=32) adapters. +:: Source +++ rps + |% + :: +posit-lattice: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-lattice:rps:unumrand (from-atom:seed:rand %sm64 0)) + :: 0x7b1d.cdaf + :: Source + ++ posit-lattice |=(r=rng:rand (pl r 5)) + :: +posit-unit: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-unit:rps:unumrand (from-atom:seed:rand %sm64 0)) + :: 0x35cf.13cd + :: Source + ++ posit-unit |=(r=rng:rand (pu r 5 128)) + -- +:: +rpd: posit64 (n=64) adapter -- +posit-lattice only (see header). +:: Source +++ rpd + |% + :: +posit-lattice: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-lattice:rpd:unumrand (from-atom:seed:rand %sm64 0)) + :: 0xe220.a839.7b1d.cdaf + :: Source + ++ posit-lattice |=(r=rng:rand (pl r 6)) + -- +:: +rpq: posit128 (n=128) adapter -- +posit-lattice only (see header). +:: Source +++ rpq + |% + :: +posit-lattice: rng -> [@ rng] + :: Examples + :: > `@ux`out:(posit-lattice:rpq:unumrand (from-atom:seed:rand %sm64 0)) + :: 0x6e78.9e6a.a1b9.65f4.e220.a839.7b1d.cdaf + :: Source + ++ posit-lattice |=(r=rng:rand (pl r 7)) + -- +-- diff --git a/librand/desk/sur/rand.hoon b/librand/desk/sur/rand.hoon new file mode 100644 index 0000000..e841984 --- /dev/null +++ b/librand/desk/sur/rand.hoon @@ -0,0 +1,23 @@ + :: /sur/rand +:::: Types for /lib/rand and its adapters (i754rand, twocrand, fixedrand, +:::: complexrand, unumrand): the engine states and the tagged +$rng union. +:: +|% +:: $phil: Philox key/counter state +:: +:: ctr is the 128-bit counter as a single @, key is the 64-bit key as a +:: single @. Both are plain atoms; width discipline is enforced by masking, +:: never by aura tricks. ++$ phil [key=@ ctr=@] +:: $sm64: SplitMix64 state (64 bits) ++$ sm64 @ +:: $pcg64: PCG64 state (128-bit state, 128-bit odd increment) ++$ pcg64 [state=@ inc=@] +:: $rng: a generic stream -- a tagged union so distribution code is +:: engine-agnostic. ++$ rng + $% [%phil p=phil] + [%sm64 s=sm64] + [%pcg p=pcg64] + == +-- diff --git a/librand/desk/tests/lib/complexrand.hoon b/librand/desk/tests/lib/complexrand.hoon new file mode 100644 index 0000000..cab091e --- /dev/null +++ b/librand/desk/tests/lib/complexrand.hoon @@ -0,0 +1,118 @@ + :: /tests/lib/complexrand +:::: +:: /lib/complexrand: value regressions for both width doors (+cd, the +:: @rd/@cd reference precision, and +cs, the @rs/@cs mirror) against +:: known-good ship output, plus sanity checks that don't depend on any +:: particular RNG draw: +on-circle lands exactly on the unit circle +:: (re^2+im^2 = 1) and +in-disk lands strictly inside the unit disk +:: (re^2+im^2 < 1). +:: +:: +cnormal:cs is also a regression for a real bug this file's own +:: construction caught: /lib/math's @rs +invsqt2 was missing a leading +:: `.0.` (parsed as the integer 70710677 instead of 0.70710677) -- see +:: libmath's tests/lib/math-constants.hoon for the fix itself. +:: +/+ *test, + rand, + i754rand, + complex, + math, + complexrand +|% +++ test-cuniform ^- tang + ;: weld + %+ expect-eq + !>(`@`0x3fe8.9e6a.a1b9.65f4.3f95.072f.63b9.b5e0) + !>(out:(cuniform:cd:complexrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3f39.65f4.3dee.6d78) + !>(out:(cuniform:cs:complexrand (from-atom:seed:rand %sm64 0))) + == +:: +++ test-normal-parts ^- tang + ;: weld + %+ expect-eq + !>(`@`0x3ff2.e308.2bd4.88f8.bfee.55f7.4af6.d8b9) + !>(out:(normal-parts:cd:complexrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x4009.430e.bf17.e80f) + !>(out:(normal-parts:cs:complexrand (from-atom:seed:rand %sm64 0))) + == +:: +++ test-cnormal ^- tang + ;: weld + %+ expect-eq + !>(`@`0x3fea.b5c4.88e8.bf41.bfe5.735e.01a0.2b7b) + !>(out:(cnormal:cd:complexrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3fc2.1e20.bed6.d405) + !>(out:(cnormal:cs:complexrand (from-atom:seed:rand %sm64 0))) + == +:: +++ test-on-circle ^- tang + ;: weld + %+ expect-eq + !>(`@`0x3fc0.7838.f0db.1ef6.3fef.bbe7.a757.e0e1) + !>(out:(on-circle:cd:complexrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3f2b.0087.3f3e.82b8) + !>(out:(on-circle:cs:complexrand (from-atom:seed:rand %sm64 0))) + == +:: +++ test-in-disk ^- tang + ;: weld + %+ expect-eq + !>(`@`0xbfd1.1d5e.36cd.f850.bfe7.45ce.ffed.7562) + !>(out:(in-disk:cd:complexrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3ee5.97d0.bf44.64a2) + !>(out:(in-disk:cs:complexrand (from-atom:seed:rand %sm64 0))) + == +:: +:: +on-circle lands exactly on the unit circle regardless of seed: no +:: rejection loop, no accumulated error beyond the underlying cos/sin's +:: own correctly-rounded precision. Checked across several distinct +:: seeds, not just the one regression value above. +++ on-circle-mag2 + |= s=@ + ^- @rd + =/ m ~(. rd:math [%n .~1e-13 .~0]) + =/ z out:(on-circle:cd:complexrand (from-atom:seed:rand %sm64 s)) + =/ re (~(re cd:complex %n) z) + =/ im (~(im cd:complex %n) z) + (add:m (mul:m re re) (mul:m im im)) +++ test-on-circle-unit-magnitude ^- tang + =/ m ~(. rd:math [%n .~1e-13 .~0]) + ;: weld + %+ expect-eq !>(%.y) !>((gth:m (on-circle-mag2 0) .~0.999999999)) + %+ expect-eq !>(%.y) !>((lth:m (on-circle-mag2 0) .~1.000000001)) + %+ expect-eq !>(%.y) !>((gth:m (on-circle-mag2 1) .~0.999999999)) + %+ expect-eq !>(%.y) !>((lth:m (on-circle-mag2 1) .~1.000000001)) + %+ expect-eq !>(%.y) !>((gth:m (on-circle-mag2 2) .~0.999999999)) + %+ expect-eq !>(%.y) !>((lth:m (on-circle-mag2 2) .~1.000000001)) + %+ expect-eq !>(%.y) !>((gth:m (on-circle-mag2 3) .~0.999999999)) + %+ expect-eq !>(%.y) !>((lth:m (on-circle-mag2 3) .~1.000000001)) + %+ expect-eq !>(%.y) !>((gth:m (on-circle-mag2 4) .~0.999999999)) + %+ expect-eq !>(%.y) !>((lth:m (on-circle-mag2 4) .~1.000000001)) + == +:: +:: +in-disk always lands strictly inside the unit disk by construction +:: (the rejection loop's own accept condition), across several seeds. +++ in-disk-mag2 + |= s=@ + ^- @rd + =/ m ~(. rd:math [%n .~1e-13 .~0]) + =/ z out:(in-disk:cd:complexrand (from-atom:seed:rand %sm64 s)) + =/ re (~(re cd:complex %n) z) + =/ im (~(im cd:complex %n) z) + (add:m (mul:m re re) (mul:m im im)) +++ test-in-disk-inside-unit-disk ^- tang + =/ m ~(. rd:math [%n .~1e-13 .~0]) + ;: weld + %+ expect-eq !>(%.y) !>((lth:m (in-disk-mag2 0) .~1)) + %+ expect-eq !>(%.y) !>((lth:m (in-disk-mag2 1) .~1)) + %+ expect-eq !>(%.y) !>((lth:m (in-disk-mag2 2) .~1)) + %+ expect-eq !>(%.y) !>((lth:m (in-disk-mag2 3) .~1)) + %+ expect-eq !>(%.y) !>((lth:m (in-disk-mag2 4) .~1)) + == +-- diff --git a/librand/desk/tests/lib/fixedrand.hoon b/librand/desk/tests/lib/fixedrand.hoon new file mode 100644 index 0000000..87ca169 --- /dev/null +++ b/librand/desk/tests/lib/fixedrand.hoon @@ -0,0 +1,55 @@ + :: /tests/lib/fixedrand +:::: +:: /lib/fixedrand: +fixed (a trivial re-export of +bits:uni:rand at the +:: precision's own width, checked for exact equality with it), +fixed-unit +:: and +fixed-between (value regressions against known-good ship output), +:: and one test proving the "sample at @rd via /lib/i754rand, quantize via +:: /lib/fixed's +from-rd" composition end to end (rand-spec.md section +:: 12.1's distribution story for fixed-point, per the "document + test, +:: don't ship dedicated wrapper arms" decision -- see librand/NEXT-STEPS.md). +:: +/+ *test, + rand, + twocrand, + fixedrand, + i754rand, + fixed +^| +=/ q88=prec:fixed [8 8] +|% +:: +fixed is exactly +bits:uni:rand at width (wid q88) -- same seed, same +:: result and resulting rng, by construction (not a coincidence to regress +:: against, but worth asserting so a future refactor that breaks the +:: passthrough is caught, same rationale as twocrand's own +twoc-full test). +++ test-fixed-is-bits ^- tang + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) (wid:fixed q88))) + !>((fixed:fixedrand (from-atom:seed:rand %sm64 0) q88)) +:: seed %sm64 0, q8.8: known-good ship output. +++ test-fixed-unit ^- tang + %+ expect-eq + !>(`@`175) + !>(out:(fixed-unit:fixedrand (from-atom:seed:rand %sm64 0) q88)) +:: range [-1.0, 3.0] in q8.8 (0x1.ff00, 0x300): known-good ship output, +:: delegates to +twoc-between at width 17 (rand-spec.md section 12.2). +++ test-fixed-between ^- tang + %+ expect-eq + !>(`@`649) + !>(out:(fixed-between:fixedrand (from-atom:seed:rand %sm64 0) q88 0x1.ff00 0x300)) +:: a > b in twoc order crashes, same as the +twoc-between it delegates to. +++ test-fixed-between-bad-range ^- tang + %- expect-fail + |.((fixed-between:fixedrand (from-atom:seed:rand %sm64 0) q88 0x300 0x1.ff00)) +:: The distribution composition rand-spec.md 12.1 describes but this library +:: doesn't wrap: draw a standard normal deviate at @rd via /lib/i754rand, +:: then quantize it to q8.8 via /lib/fixed's +from-rd. Both steps are +:: independently tested elsewhere (i754rand's own +test-normal, fixed's own +:: +test-rd-bridge); this proves the COMPOSITION against a known-good +:: ship-verified result, not either step in isolation. +++ test-fixed-dist-composition ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ z r (normal:rd:dist:i754rand r) + %+ expect-eq + !>(`@`130.829) + !>((from-rd:fixed z q88)) +-- diff --git a/librand/desk/tests/lib/i754rand.hoon b/librand/desk/tests/lib/i754rand.hoon new file mode 100644 index 0000000..bbf6a4f --- /dev/null +++ b/librand/desk/tests/lib/i754rand.hoon @@ -0,0 +1,300 @@ + :: /tests/lib/i754rand +:::: +:: /lib/i754rand (the porcelain): ++uni (float: +rs/+rd/+rh/+rq/+rs-oo/ +:: +rd-oo), ++alias (Vose's method), ++dist (++rd + ++rs). ++uni's +:: float arms are exact bit constructions (+sun then a power-of-two +:: multiply), so their expected values are checked against an +:: independent Python IEEE-754 encoder, not the Hoon float printer. +:: ++dist's value spot-checks (including the alias table) involve +:: /lib/math transcendentals (log/sqrt), so their expected values are +:: read off this same ship's actual output (sanity-checked by hand +:: against the closed-form algorithm first) rather than independently +:: rederived in Python, which isn't guaranteed to match this +:: codebase's correctly-rounded kernels to the last bit; ++dist's +:: moment tests are regression checks against the THEORETICAL targets +:: (mean/variance) at a fixed seed, with tolerances wide enough to +:: absorb real sampling noise at n=50000 but tight enough to catch an +:: actual algorithm bug. Two real bugs were caught by this on-ship +:: testing during development (alias-table's worklist assignment +:: reversed; +binomial comparing against the wrong accumulator) -- +:: see NEXT-STEPS.md. +:: +/+ *test, + rand, + i754rand, + math +|% +:: +rs/+rd/+rh/+rq: exact bit constructions, checked against an +:: independent Python IEEE-754 encoder (not the Hoon float printer). +++ test-float-rs ^- tang + %+ expect-eq + !>(`@rs`0x3dee.6d78) + !>(out:(rs:uni:i754rand (from-atom:seed:rand %sm64 0))) +++ test-float-rd ^- tang + %+ expect-eq + !>(`@rd`0x3f95.072f.63b9.b5e0) + !>(out:(rd:uni:i754rand (from-atom:seed:rand %sm64 0))) +++ test-float-rh ^- tang + %+ expect-eq + !>(`@rh`0x39af) + !>(out:(rh:uni:i754rand (from-atom:seed:rand %sm64 0))) +++ test-float-rq ^- tang + %+ expect-eq + !>(`@rq`0x3ffd.3cd5.4372.cbe9.c441.5072.f63b.9b5e) + !>(out:(rq:uni:i754rand (from-atom:seed:rand %sm64 0))) +++ test-float-rs-oo ^- tang + %+ expect-eq + !>(`@rs`0x3dee.6d7c) + !>(out:(rs-oo:uni:i754rand (from-atom:seed:rand %sm64 0))) +++ test-float-rd-oo ^- tang + %+ expect-eq + !>(`@rd`0x3f95.072f.63b9.b5f0) + !>(out:(rd-oo:uni:i754rand (from-atom:seed:rand %sm64 0))) +:: ++dist: value spot-checks against actual on-ship computation (these +:: involve /lib/math transcendentals -- log/sqrt -- so the expected +:: values below are read off this same ship's output, not independently +:: rederived in Python, which isn't guaranteed to match this codebase's +:: correctly-rounded kernels to the last bit). Each is sanity-checked +:: against the closed-form algorithm by hand (e.g. -ln(u)/lambda for +:: +expon) before being trusted as a regression constant. +++ test-normal ^- tang + %+ expect-eq + !>(`@rd`.~-0.9479938949723624) + !>(out:(normal:rd:dist:i754rand (from-atom:seed:rand %sm64 0))) +:: determinism: same seed -> same output (spot check; +normal's rejection +:: loop is the one arm in ++dist where a threading bug would most likely +:: show up as nondeterminism). +++ test-normal-determinism ^- tang + %+ expect-eq + !>((normal:rd:dist:i754rand (from-atom:seed:rand %phil 7))) + !>((normal:rd:dist:i754rand (from-atom:seed:rand %phil 7))) +++ test-normal-mv ^- tang + %+ expect-eq + !>(`@rd`.~8.104012210055275) + !>(out:(normal-mv:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~10 .~2)) +++ test-normal-mv-bad-sigma ^- tang + %- expect-fail + |.((normal-mv:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0 .~-1)) +:: -ln(u)/lambda, u = +rd-oo of the first SplitMix64 output (0.0205...): +:: -ln(0.0205352...)/2 ~ 1.9428 -- matches to displayed precision. +++ test-expon ^- tang + %+ expect-eq + !>(`@rd`.~1.942806871605541) + !>(out:(expon:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~2)) +++ test-expon-bad-rate ^- tang + %- expect-fail + |.((expon:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0)) +:: alpha=2 (>=1, direct Marsaglia-Tsang path). +++ test-gamma-ge1 ^- tang + %+ expect-eq + !>(`@rd`.~0.7179344171650754) + !>(out:(gamma:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~2)) +:: alpha=0.5 (<1, boost path: gamma(1.5) * u^(1/0.5)). +++ test-gamma-lt1 ^- tang + %+ expect-eq + !>(`@rd`.~0.05447897345403265) + !>(out:(gamma:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0.5)) +++ test-gamma-bad-shape ^- tang + %- expect-fail + |.((gamma:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0)) +++ test-beta ^- tang + %+ expect-eq + !>(`@rd`.~0.23622515961139348) + !>(out:(beta:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~2 .~3)) +++ test-chi2 ^- tang + %+ expect-eq + !>(`@rd`.~0.8261342997997718) + !>(out:(chi2:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 3)) +++ test-chi2-bad-df ^- tang + %- expect-fail + |.((chi2:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 0)) +++ test-student-t ^- tang + %+ expect-eq + !>(`@rd`.~-0.7137613421937223) + !>(out:(student-t:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 5)) +++ test-student-t-bad-df ^- tang + %- expect-fail + |.((student-t:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 0)) +:: u ~ 0.0205 < 0.5 -> %.y. +++ test-bernoulli ^- tang + %+ expect-eq + !>(`?`%.y) + !>(out:(bernoulli:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0.5)) +++ test-bernoulli-bad-prob ^- tang + %- expect-fail + |.((bernoulli:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~1.5)) +++ test-geometric ^- tang + %+ expect-eq + !>(`@ud`11) + !>(out:(geometric:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0.3)) +:: p=1 edge case: returns 1 without consuming a draw (rng unchanged). +++ test-geometric-p1 ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ out r (geometric:rd:dist:i754rand r .~1) + ;: weld + %+ expect-eq !>(`@ud`1) !>(out) + %+ expect-eq !>((from-atom:seed:rand %sm64 0)) !>(r) + == +++ test-geometric-bad-prob ^- tang + %- expect-fail + |.((geometric:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0)) +:: Moment tests (rand-spec.md section 10 item 4): sample mean/variance at +:: 50k draws, FIXED seed -- deterministic, so this is a regression test +:: against the theoretical constants, not a flaky statistical one. +:: Tolerances are generous relative to the actual sampling error at +:: n=50000 (stderr(mean) ~ sigma/sqrt(n), stderr(var) ~ sigma^2*sqrt(2/n)) +:: so a real algorithm bug (wrong constant, sign error, wrong formula) +:: fails loudly while ordinary run-to-run noise from a fixed seed does not +:: (there being only one fixed seed, "run-to-run" here really means +:: "robust to this exact computation shifting by a few ULP if the +:: underlying math kernels ever change"). +++ test-moments-normal ^- tang + =/ mv (moments (from-atom:seed:rand %phil 0) normal:rd:dist:i754rand 50.000) + ;: weld + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.02]) -.mv .~0)) + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.05]) +.mv .~1)) + == +++ test-moments-expon ^- tang + =/ mv (moments (from-atom:seed:rand %phil 0) |=(r=rng:rand (expon:rd:dist:i754rand r .~1)) 50.000) + ;: weld + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.05]) -.mv .~1)) + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.05]) +.mv .~1)) + == +++ test-moments-gamma ^- tang + =/ mv (moments (from-atom:seed:rand %phil 0) |=(r=rng:rand (gamma:rd:dist:i754rand r .~2)) 50.000) + ;: weld + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.05]) -.mv .~2)) + %+ expect-eq !>(%.y) !>((~(is-close rd:math [%n .~0 .~0.05]) +.mv .~2)) + == +:: +moments: draw .n samples via .f from .r0, return [mean=@rd var=@rd] +:: (population variance). Shared by the three moment tests above. +++ moments + |= [r0=rng:rand f=$-(rng:rand [@rd rng:rand]) n=@] + ^- [mean=@rd var=@rd] + =/ mth ~(. rd:math [%n .~1e-13 .~0]) + =/ r r0 + =/ i 0 + =/ sum .~0 + =/ ssq .~0 + |- ^- [@rd @rd] + ?: =(i n) + =/ ct (sun:mth n) + =/ mean (div:mth sum ct) + [mean (sub:mth (div:mth ssq ct) (mul:mth mean mean))] + =^ x r (f r) + %= $ + i +(i) + sum (add:mth sum x) + ssq (add:mth ssq (mul:mth x x)) + == +:: +alias +build: worked example matching rand-spec.md section 6.1's own +:: doc comment, hand-verified (weights [1,1,2]/4 -> true probabilities +:: [0.25,0.25,0.5]; this table's prob/alias reproduce that -- see the +:: arm's doccord for the by-hand trace). This test caught a real bug +:: during development: the small/large worklist assignment was reversed, +:: which +build-loop's own doc comment now explains. +++ test-alias-build ^- tang + %+ expect-eq + !>(`alias-table:i754rand`[n=3 prob=~[.~0.75 .~0.75 .~1] alias=~[2 2 2]]) + !>((build:alias:i754rand ~[.~1 .~1 .~2])) +:: +alias +build: the other worklist-draining order (large empties last). +++ test-alias-build-other-order ^- tang + %+ expect-eq + !>(`alias-table:i754rand`[n=4 prob=~[.~1 .~0.6666666666666666 .~0.6666666666666666 .~0.6666666666666666] alias=~[0 0 0 0]]) + !>((build:alias:i754rand ~[.~3 .~1 .~1 .~1])) +++ test-alias-build-empty ^- tang + %- expect-fail + |.((build:alias:i754rand ~)) +++ test-alias-build-bad-prob ^- tang + %- expect-fail + |.((build:alias:i754rand ~[.~0 .~0])) +:: +alias +draw: value regression, from the table above. +++ test-alias-draw ^- tang + =/ t (build:alias:i754rand ~[.~1 .~1 .~2]) + %+ expect-eq + !>(`[out=@ud rng:rand]`[2 [%sm64 s=0x3c6e.f372.fe94.f82a]]) + !>((draw:alias:i754rand t (from-atom:seed:rand %sm64 0))) +:: +categorical: thin wrapper over +draw:alias -- same table, same +:: result, by construction. +++ test-categorical ^- tang + =/ t (build:alias:i754rand ~[.~1 .~1 .~2]) + %+ expect-eq + !>((draw:alias:i754rand t (from-atom:seed:rand %sm64 0))) + !>((categorical:rd:dist:i754rand (from-atom:seed:rand %sm64 0) t)) +:: +poisson: Knuth branch (lambda<10) and PTRS branch (lambda>=10) value +:: regressions, plus determinism-implying mean sanity already covered by +:: hand-verified means during development (2000-draw sample means came +:: out within 0.5% of lambda for both 5 and 20 -- not re-asserted here as +:: an exact regression to keep this test fast). Crashes if lambda<=0. +++ test-poisson-knuth ^- tang + %+ expect-eq + !>(`@ud`2) + !>(-:(poisson:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~5)) +++ test-poisson-ptrs ^- tang + %+ expect-eq + !>(`@ud`14) + !>(-:(poisson:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~20)) +++ test-poisson-bad-rate ^- tang + %- expect-fail + |.((poisson:rd:dist:i754rand (from-atom:seed:rand %sm64 0) .~0)) +:: +binomial: value regression (also exercises the fix for a real bug +:: caught during development -- the accept test originally compared the +:: uniform draw against the individual pmf term instead of the +:: accumulated CDF, which either returned 0 too often or ran the +:: recurrence past n and crashed with subtract-underflow; a 2000-draw +:: sample mean check during development came out at 5.98 vs the +:: theoretical n*p=6). Crashes if n*min(p,1-p) >= 30. +++ test-binomial ^- tang + %+ expect-eq + !>(`@ud`2) + !>(-:(binomial:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 20 .~0.3)) +++ test-binomial-n-too-large ^- tang + %- expect-fail + |.((binomial:rd:dist:i754rand (from-atom:seed:rand %sm64 0) 1.000.000 .~0.5)) +:: +dirichlet: value regression (sums to 1 by construction; spot-check +:: the sum explicitly too, since that's the actual mathematical +:: invariant, not just a frozen regression number). Crashes on an +:: empty alphas list. +++ test-dirichlet ^- tang + =/ out `(list @rd)`-:(dirichlet:rd:dist:i754rand (from-atom:seed:rand %sm64 0) ~[.~1 .~2 .~3]) + =/ mth ~(. rd:math [%n .~1e-13 .~0]) + =/ total (add:mth (snag 0 out) (add:mth (snag 1 out) (snag 2 out))) + ;: weld + %+ expect-eq + !>(`(list @rd)`~[.~0.04688412483230968 .~0.4265186606662172 .~0.5265972145014731]) + !>(out) + %+ expect-eq !>(%.y) !>((is-close:mth total .~1)) + == +++ test-dirichlet-empty ^- tang + %- expect-fail + |.((dirichlet:rd:dist:i754rand (from-atom:seed:rand %sm64 0) ~)) +:: ++rs mirror: representative spot-checks (not exhaustive re-coverage +:: of every ++rd test) confirming the mechanical @rd->@rs +:: re-instantiation actually produces working, independently-computed +:: @rs output -- not just that it compiles. +++ test-rs-normal ^- tang + %+ expect-eq + !>(`@rs`.-0.5933847) + !>(out:(normal:rs:dist:i754rand (from-atom:seed:rand %sm64 0))) +++ test-rs-gamma ^- tang + %+ expect-eq + !>(`@rs`.1.0119848) + !>(out:(gamma:rs:dist:i754rand (from-atom:seed:rand %sm64 0) .2)) +++ test-rs-expon ^- tang + %+ expect-eq + !>(`@rs`.1.0752765) + !>(out:(expon:rs:dist:i754rand (from-atom:seed:rand %sm64 0) .2)) +++ test-rs-poisson-ptrs ^- tang + %+ expect-eq + !>(`@ud`14) + !>(out:(poisson:rs:dist:i754rand (from-atom:seed:rand %sm64 0) .20)) +++ test-rs-binomial ^- tang + %+ expect-eq + !>(`@ud`4) + !>(out:(binomial:rs:dist:i754rand (from-atom:seed:rand %sm64 0) 20 .0.3)) +++ test-rs-dirichlet ^- tang + =/ out `(list @rs)`out:(dirichlet:rs:dist:i754rand (from-atom:seed:rand %sm64 0) ~[.1 .2 .3]) + =/ mth ~(. rs:math [%n .1e-6 .0]) + =/ total (add:mth (snag 0 out) (add:mth (snag 1 out) (snag 2 out))) + %+ expect-eq !>(%.y) !>((is-close:mth total .1)) +-- diff --git a/librand/desk/tests/lib/rand.hoon b/librand/desk/tests/lib/rand.hoon new file mode 100644 index 0000000..55c50d8 --- /dev/null +++ b/librand/desk/tests/lib/rand.hoon @@ -0,0 +1,368 @@ + :: /tests/lib/rand +:::: +:: /lib/rand (the plumbing): ++split-mix, ++philox, ++seed, +fork, +:: ++gen, ++uni (integer: +bits/+below/+between), ++pcg, ++sample +:: (shuffle/permutation/choice/reservoir -- +alias moved to +:: /lib/i754rand, tested there instead). The SplitMix64 and +:: Philox4x32-10 KAT vectors are checked against their published +:: reference values (Vigna's reference C for SplitMix64; the Random123 +:: kat_vectors file -- all-zero, all-0xffffffff, and the pi-digits +:: vector -- for Philox), each cross-checked against an independent +:: Python re-implementation before being transcribed here. ++pcg's +:: KAT (seed=42, seq=54, pcg-c's own demo convention) is checked +:: against an independent Python re-implementation of pcg-c's +:: reference C (include/pcg_variants.h) -- this also exercises the +:: CORRECTED advance-then-output order (rand-spec.md section 3.3; an +:: earlier spec draft had it backwards). +from-atom, +fold-wide, +:: +mix, +fork, +bits, and +below have no external reference (they +:: are this library's own design, or Lemire's algorithm applied to +:: this library's own +step), so those are checked by direct +:: computation of the spec'd formula and by determinism/order- +:: sensitivity/path-sensitivity/unbiasedness properties instead. +:: +/+ *test, + rand +|% +:: SplitMix64 KAT: first 5 outputs from seed 0, against Vigna's reference C. +++ test-splitmix-kat-seed-0 ^- tang + =/ s0=@ 0 + =^ o0 s0 (next:split-mix:rand s0) + =^ o1 s0 (next:split-mix:rand s0) + =^ o2 s0 (next:split-mix:rand s0) + =^ o3 s0 (next:split-mix:rand s0) + =^ o4 s0 (next:split-mix:rand s0) + ;: weld + %+ expect-eq !>(`@`0xe220.a839.7b1d.cdaf) !>(o0) + %+ expect-eq !>(`@`0x6e78.9e6a.a1b9.65f4) !>(o1) + %+ expect-eq !>(`@`0x6c4.5d18.8009.454f) !>(o2) + %+ expect-eq !>(`@`0xf88b.b8a8.724c.81ec) !>(o3) + %+ expect-eq !>(`@`0x1b39.896a.51a8.749b) !>(o4) + == +:: SplitMix64 KAT: first 5 outputs from seed 0xdeadbeef. +++ test-splitmix-kat-seed-deadbeef ^- tang + =/ s0=@ 0xdead.beef + =^ o0 s0 (next:split-mix:rand s0) + =^ o1 s0 (next:split-mix:rand s0) + =^ o2 s0 (next:split-mix:rand s0) + =^ o3 s0 (next:split-mix:rand s0) + =^ o4 s0 (next:split-mix:rand s0) + ;: weld + %+ expect-eq !>(`@`0x4adf.b90f.68c9.eb9b) !>(o0) + %+ expect-eq !>(`@`0xde58.6a31.41a1.0922) !>(o1) + %+ expect-eq !>(`@`0x21f.bc2f.8e1c.fc1d) !>(o2) + %+ expect-eq !>(`@`0x7466.ce73.7be1.6790) !>(o3) + %+ expect-eq !>(`@`0x3bfa.8764.f685.bd1c) !>(o4) + == +:: +split is two +next draws used directly as child seeds. +++ test-split ^- tang + =/ r (split:split-mix:rand 0) + ;: weld + %+ expect-eq !>(`@`0xe220.a839.7b1d.cdaf) !>(a.r) + %+ expect-eq !>(`@`0x6e78.9e6a.a1b9.65f4) !>(b.r) + %+ expect-eq !>(`@`0x3c6e.f372.fe94.f82a) !>(s.r) + == +:: +from-atom %sm64: state = seed directly, no pre-draw. +++ test-from-atom-sm64 ^- tang + %+ expect-eq + !>(`rng:rand`[%sm64 s=0]) + !>((from-atom:seed:rand %sm64 0)) +:: +from-atom %phil: key = first SplitMix output, ctr reset to 0. +++ test-from-atom-phil ^- tang + %+ expect-eq + !>(`rng:rand`[%phil p=[key=0xe220.a839.7b1d.cdaf ctr=0]]) + !>((from-atom:seed:rand %phil 0)) +:: +from-atom %pcg: state/inc from two outputs, inc forced odd. +++ test-from-atom-pcg ^- tang + %+ expect-eq + !>(`rng:rand`[%pcg p=[state=0xe220.a839.7b1d.cdaf inc=0x6e78.9e6a.a1b9.65f5]]) + !>((from-atom:seed:rand %pcg 0)) +:: same seed -> identical rng noun, for all three engines (determinism). +++ test-from-atom-determinism ^- tang + ;: weld + %+ expect-eq !>((from-atom:seed:rand %sm64 42)) !>((from-atom:seed:rand %sm64 42)) + %+ expect-eq !>((from-atom:seed:rand %phil 42)) !>((from-atom:seed:rand %phil 42)) + %+ expect-eq !>((from-atom:seed:rand %pcg 42)) !>((from-atom:seed:rand %pcg 42)) + == +:: +mix: known values from the spec'd two-word compression, plus the +:: documented mix(0,0) = 0 degenerate point (see +mix's doc comment). +++ test-mix-values ^- tang + ;: weld + %+ expect-eq !>(`@`0xef30.b01c.2974.aeeb) !>((mix:seed:rand 1 2)) + %+ expect-eq !>(`@`0x3ec2.d42f.3a45.cc6e) !>((mix:seed:rand 2 1)) + %+ expect-eq !>(`@`0x0) !>((mix:seed:rand 0 0)) + == +:: order-sensitivity: (mix a b) != (mix b a) in general -- the whole point +:: of the two-word compression (rand-spec.md section 4), load-bearing for +:: +fork's path-sensitivity property once +fork lands in milestone 2. +++ test-mix-order-sensitive ^- tang + %+ expect-eq !>(%.n) !>(=((mix:seed:rand 42 7) (mix:seed:rand 7 42))) +:: +from-eny: determinism, and a +fold-wide KAT against the same reference +:: fold computed independently in Python over 8 known 64-bit words. +++ test-from-eny ^- tang + =/ e=@uvJ + `@uvJ`0x8888.8888.8888.8888.7777.7777.7777.7777.6666.6666.6666.6666.5555.5555.5555.5555.4444.4444.4444.4444.3333.3333.3333.3333.2222.2222.2222.2222.1111.1111.1111.1111 + ;: weld + %+ expect-eq !>((from-eny:seed:rand %sm64 e)) !>((from-eny:seed:rand %sm64 e)) + %+ expect-eq !>(`@`0xdeb.8a1b.ec35.d57d) !>((fold-wide:seed:rand e)) + == +:: Philox4x32-10 KAT: Random123 kat_vectors, all-zero case. +++ test-philox-kat-zero ^- tang + %+ expect-eq + !>(`@`0x9b00.dbd8.bc57.ac4c.e169.c58d.6627.e8d5) + !>((block:philox:rand 0 0)) +:: Philox4x32-10 KAT: Random123 kat_vectors, all-0xffffffff case. +++ test-philox-kat-ones ^- tang + %+ expect-eq + !>(`@`0x6d54.51fd.a20b.c7c6.41c8.3b0e.408f.276d) + !>((block:philox:rand 0xffff.ffff.ffff.ffff 0xffff.ffff.ffff.ffff.ffff.ffff.ffff.ffff)) +:: Philox4x32-10 KAT: Random123 kat_vectors, pi-digits case. +++ test-philox-kat-pi ^- tang + %+ expect-eq + !>(`@`0x2412.6ea1.5001.e420.94fd.cceb.d16c.fe09) + !>((block:philox:rand 0x299f.31d0.a409.3822 0x370.7344.1319.8a2e.85a3.08d3.243f.6a88)) +:: +next:philox: out = (block key ctr), ctr advances by 1, key unchanged. +++ test-philox-next ^- tang + %+ expect-eq + !>(`[out=@ p=phil:rand]`[0x9b00.dbd8.bc57.ac4c.e169.c58d.6627.e8d5 [key=0 ctr=1]]) + !>((next:philox:rand [key=0 ctr=0])) +:: +step dispatches to the right engine and keeps only bits [0,64) of a +:: Philox block. +++ test-step-phil ^- tang + =/ r (from-atom:seed:rand %phil 0) + =^ out r (step:rand r) + ;: weld + %+ expect-eq !>(`@`0x24cb.d2fb.a9e3.9636) !>(out) + %+ expect-eq !>(`rng:rand`[%phil p=[key=0xe220.a839.7b1d.cdaf ctr=1]]) !>(r) + == +++ test-step-sm64 ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ out r (step:rand r) + ;: weld + %+ expect-eq !>(`@`0xe220.a839.7b1d.cdaf) !>(out) + %+ expect-eq !>(`rng:rand`[%sm64 s=0x9e37.79b9.7f4a.7c15]) !>(r) + == +:: +fork: deterministic (same parent + salt -> same child). +++ test-fork-determinism ^- tang + =/ r (from-atom:seed:rand %phil 0) + %+ expect-eq + !>((fork:rand r 1)) + !>((fork:rand r 1)) +:: +fork: different salts -> different children, for all three engines. +++ test-fork-distinct-salts ^- tang + ;: weld + %+ expect-eq !>(%.n) + !> =((fork:rand (from-atom:seed:rand %phil 0) 1) (fork:rand (from-atom:seed:rand %phil 0) 2)) + %+ expect-eq !>(%.n) + !> =((fork:rand (from-atom:seed:rand %pcg 0) 1) (fork:rand (from-atom:seed:rand %pcg 0) 2)) + %+ expect-eq !>(%.n) + !> =((fork:rand (from-atom:seed:rand %sm64 0) 1) (fork:rand (from-atom:seed:rand %sm64 0) 2)) + == +:: +fork: %phil's ctr always resets to 0 in the child, regardless of the +:: parent's ctr. +++ test-fork-phil-ctr-reset ^- tang + =/ r [%phil p=[key=0x2a ctr=99]] + =/ child (fork:rand r 7) + ?> ?=(%phil -.child) + %+ expect-eq !>(0) !>(ctr.p.child) +:: +fork: %pcg's child increment is always forced odd. +++ test-fork-pcg-inc-odd ^- tang + =/ r (from-atom:seed:rand %pcg 0) + =/ child (fork:rand r 5) + ?> ?=(%pcg -.child) + %+ expect-eq !>(1) !>((dis 1 inc.p.child)) +:: Nesting rule (rand-spec.md section 2.1c): the mix is genuinely +:: path-sensitive. (fork (fork r a) b) must differ from +:: (fork (fork r b) a) and from (fork r (cat 6 a b)). +++ test-fork-path-sensitive ^- tang + =/ r (from-atom:seed:rand %phil 0) + =/ a 11 + =/ b 22 + =/ fab (fork:rand (fork:rand r a) b) + =/ fba (fork:rand (fork:rand r b) a) + =/ fcat (fork:rand r (cat 6 a b)) + ;: weld + %+ expect-eq !>(%.n) !>(=(fab fba)) + %+ expect-eq !>(%.n) !>(=(fab fcat)) + %+ expect-eq !>(%.n) !>(=(fba fcat)) + == +:: ++gen door facade: +draw mirrors the functional +step, +fork mirrors +:: the functional +fork. +++ test-gen-draw ^- tang + =/ g ~(. gen:rand (from-atom:seed:rand %sm64 0)) + =^ x g draw:g + %+ expect-eq !>(`@`0xe220.a839.7b1d.cdaf) !>(x) +++ test-gen-fork ^- tang + =/ g ~(. gen:rand (from-atom:seed:rand %phil 0)) + %+ expect-eq + !>((fork:rand (from-atom:seed:rand %phil 0) 3)) + !>(r:(fork:g 3)) +:: +bits: n=4 -> low 4 bits of the first +step draw (0xe220a8397b1dcdaf). +++ test-bits ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ out r (bits:uni:rand r 4) + ;: weld + %+ expect-eq !>(`@`15) !>(out) + %+ expect-eq !>(`rng:rand`[%sm64 s=0x9e37.79b9.7f4a.7c15]) !>(r) + == +:: +bits: n=130 (>64, spans 3 draws: two full 64-bit words, one 2-bit +:: remainder) round-trips against a manual reconstruction from the same +:: three +step draws taken independently. +++ test-bits-wide ^- tang + =/ r0 (from-atom:seed:rand %sm64 0) + =/ result (bits:uni:rand r0 130) + =^ w0 r0 (step:rand r0) + =^ w1 r0 (step:rand r0) + =^ w2 r0 (step:rand r0) + =/ expected :(add w0 (lsh [0 64] w1) (lsh [0 128] (end [0 2] w2))) + ;: weld + %+ expect-eq !>(expected) !>(out.result) + %+ expect-eq !>(r0) !>(r.result) + == +:: +below: Lemire, n=10, against an independently-computed reference +:: (same algorithm, computed in Python before transcription here). +++ test-below ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ out r (below:uni:rand r 10) + ;: weld + %+ expect-eq !>(`@`8) !>(out) + %+ expect-eq !>(`rng:rand`[%sm64 s=0x9e37.79b9.7f4a.7c15]) !>(r) + == +:: +below: unbiasedness smoke test, n=6, 1200 draws, chi-square below a +:: fixed threshold (section 10 item 3: a smoke test, not TestU01). Uses +:: ++fork to get 1200 independent draws without threading one rng +:: sequentially (each fork is its own single-draw +below call). Chi- +:: square is computed as an exact integer comparison (sum-sq < k*expected) +:: rather than a division, so this needs no float/math import. +++ test-below-unbiased ^- tang + =/ base (from-atom:seed:rand %phil 0) + =/ n 6 + =/ reps 1.200 + =/ expected (div reps n) + =/ counts + %+ roll (gulf 0 (dec reps)) + |= [i=@ counts=(map @ @)] + =/ child (fork:rand base i) + =^ draw child (below:uni:rand child n) + (~(put by counts) draw +((fall (~(get by counts) draw) 0))) + =/ sum-sq + %+ roll (gulf 0 (dec n)) + |= [k=@ acc=@] + =/ obs (fall (~(get by counts) k) 0) + =/ diff (abs:si (dif:si (sun:si obs) (sun:si expected))) + (add acc (mul diff diff)) + :: chi-square = sum-sq / expected; threshold chosen generously (df=5, + :: p~0.001 critical value is ~20.5) since this is a smoke test, not a + :: rigorous statistical test -- compared as sum-sq < 30*expected to + :: avoid a division. + %+ expect-eq !>(%.y) !>((lth sum-sq (mul 30 expected))) +:: +between: span=11 (a=-5,b=--5), same underlying draw as +test-below. +++ test-between ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ out r (between:uni:rand r -5 --5) + %+ expect-eq !>(`@s`--4) !>(out) +:: +between: crashes if a > b. +++ test-between-bad-range ^- tang + %- expect-fail + |.((between:uni:rand (from-atom:seed:rand %sm64 0) --5 -5)) +:: +below: crashes on n=0. +++ test-below-zero ^- tang + %- expect-fail + |.((below:uni:rand (from-atom:seed:rand %sm64 0) 0)) +:: PCG64 KAT: pcg-c's own seeding (pcg_setseq_128_srandom_r, seed=42, +:: seq=54 -- the library's classic demo parameters), first 5 outputs, +:: against an independent Python re-implementation of pcg-c's reference +:: C (include/pcg_variants.h). This exercises the CORRECTED operation +:: order (advance then output -- see rand-spec.md section 3.3), not the +:: order an earlier draft of the spec claimed. Seeded directly with the +:: literal state/inc here rather than via +from-atom, since this KAT +:: checks +next:pcg's own algorithm against pcg-c's native seeding, not +:: this library's (unrelated, SplitMix64-derived) +from-atom scheme. +++ test-pcg-kat ^- tang + =/ p [state=0xde2b.ce05.be01.3be3.d3f6.c45a.41e5.4320 inc=0x6d] + =^ o0 p (next:pcg:rand p) + =^ o1 p (next:pcg:rand p) + =^ o2 p (next:pcg:rand p) + =^ o3 p (next:pcg:rand p) + =^ o4 p (next:pcg:rand p) + ;: weld + %+ expect-eq !>(`@`0x86b1.da1d.7206.2b68) !>(o0) + %+ expect-eq !>(`@`0x1304.aa46.c985.3d39) !>(o1) + %+ expect-eq !>(`@`0xa367.0e9e.0dd5.0358) !>(o2) + %+ expect-eq !>(`@`0xf909.0e52.9a7d.ae00) !>(o3) + %+ expect-eq !>(`@`0xc85b.9fd8.3799.6f2c) !>(o4) + == +:: +advance: at a small, tractable delta (3), must equal 3 sequential +:: +next:pcg steps -- the property that justifies the O(log delta) +:: skip-ahead algorithm +jump relies on at the (untestable-by-brute- +:: force) delta = 2^64 scale. +++ test-pcg-advance ^- tang + =/ p0 [state=5 inc=0x6d] + =/ p1 (advance:pcg:rand p0 3) + =^ s1 p0 (next:pcg:rand p0) + =^ s2 p0 (next:pcg:rand p0) + =^ s3 p0 (next:pcg:rand p0) + %+ expect-eq !>(state.p0) !>(state.p1) +:: +jump: deterministic, and actually changes the state (catches a +:: no-op mistake). +++ test-pcg-jump ^- tang + =/ p [state=5 inc=0x6d] + ;: weld + %+ expect-eq !>((jump:pcg:rand p)) !>((jump:pcg:rand p)) + %+ expect-eq !>(%.n) !>(=(p (jump:pcg:rand p))) + == +:: +step dispatches %pcg correctly (next:pcg wired in, not the earlier +:: crash stub). +++ test-step-pcg ^- tang + =/ r (from-atom:seed:rand %pcg 0) + =^ out r (step:rand r) + ;: weld + %+ expect-eq !>(`@`0x517a.36a9.6d93.79b8) !>(out) + %+ expect-eq + !>(`rng:rand`[%pcg p=[state=0x15ab.4f3e.beb3.372a.3aed.9a13.0cdc.0020 inc=0x6e78.9e6a.a1b9.65f5]]) + !>(r) + == +:: +shuffle: value regression, and the permutation property (same +:: multiset) over a fresh draw at a different seed. +++ test-shuffle ^- tang + %+ expect-eq + !>(`(list @ud)`~[3 4 1 2 5]) + !>(-:(shuffle:sample:rand (from-atom:seed:rand %sm64 0) `(list @ud)`~[1 2 3 4 5])) +++ test-shuffle-same-multiset ^- tang + =/ in `(list @ud)`~[10 20 30 40 50 60] + =/ out -:(shuffle:sample:rand (from-atom:seed:rand %phil 3) in) + %+ expect-eq !>((sort in lth)) !>((sort out lth)) +:: +permutation: value regression (delegates to +shuffle over gulf 0 4). +++ test-permutation ^- tang + %+ expect-eq + !>(`(list @ud)`~[2 3 0 1 4]) + !>(-:(permutation:sample:rand (from-atom:seed:rand %sm64 0) 5)) +:: +choice: value regression, and crash on empty. +++ test-choice ^- tang + %+ expect-eq + !>(30) + !>(-:(choice:sample:rand (from-atom:seed:rand %sm64 0) `(list @ud)`~[10 20 30])) +++ test-choice-empty ^- tang + %- expect-fail + |.((choice:sample:rand (from-atom:seed:rand %sm64 0) `(list @ud)`~)) +:: +choices: value regression (with replacement -- repeats are fine). +++ test-choices ^- tang + %+ expect-eq + !>(`(list @ud)`~[30 20 10 30]) + !>(-:(choices:sample:rand (from-atom:seed:rand %sm64 0) 4 `(list @ud)`~[10 20 30])) +:: +sample-n: value regression, and crash if k exceeds the list length. +++ test-sample-n ^- tang + %+ expect-eq + !>(`(list @ud)`~[3 2 5]) + !>(-:(sample-n:sample:rand (from-atom:seed:rand %sm64 0) 3 `(list @ud)`~[1 2 3 4 5])) +++ test-sample-n-bad-count ^- tang + %- expect-fail + |.((sample-n:sample:rand (from-atom:seed:rand %sm64 0) 6 `(list @ud)`~[1 2 3 4 5])) +:: +reservoir: value regression, and crash if n exceeds the list length. +++ test-reservoir ^- tang + %+ expect-eq + !>(`(list @ud)`~[8 2 5]) + !>(-:(reservoir:sample:rand (from-atom:seed:rand %sm64 0) 3 `(list @ud)`~[1 2 3 4 5 6 7 8])) +++ test-reservoir-bad-count ^- tang + %- expect-fail + |.((reservoir:sample:rand (from-atom:seed:rand %sm64 0) 9 `(list @ud)`~[1 2 3 4 5 6 7 8])) +-- diff --git a/librand/desk/tests/lib/twocrand.hoon b/librand/desk/tests/lib/twocrand.hoon new file mode 100644 index 0000000..e7bb434 --- /dev/null +++ b/librand/desk/tests/lib/twocrand.hoon @@ -0,0 +1,31 @@ + :: /tests/lib/twocrand +:::: +:: /lib/twocrand: +twoc-full (a trivial re-export of +bits:uni:rand, +:: checked for exact equality with it -- no independent value needed +:: since it IS that arm) and +twoc-between (value regression against a +:: by-hand Lemire trace, plus a crash test for a>b in twoc order). +:: +/+ *test, + rand, + twocrand +|% +:: +twoc-full is exactly +bits:uni:rand -- same seed, same width, same +:: result and same resulting rng, by construction (not a coincidence to +:: regress against, but worth asserting so a future refactor that +:: breaks the passthrough is caught). +++ test-twoc-full-is-bits ^- tang + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 8)) + !>((twoc-full:twocrand (from-atom:seed:rand %sm64 0) 8)) +:: w=8, a=0xfb (-5), b=0x5 (5): span=11, below(11)=9 (the same known +:: Lemire trace as ++uni's own +test-between), biased back by +9 from +:: -5 lands on 4. +++ test-twoc-between ^- tang + %+ expect-eq + !>(`@`4) + !>(out:(twoc-between:twocrand (from-atom:seed:rand %sm64 0) 8 0xfb 0x5)) +:: a > b in twoc order (0x5 = 5, 0xfb = -5) crashes. +++ test-twoc-between-bad-range ^- tang + %- expect-fail + |.((twoc-between:twocrand (from-atom:seed:rand %sm64 0) 8 0x5 0xfb)) +-- diff --git a/librand/desk/tests/lib/unumrand.hoon b/librand/desk/tests/lib/unumrand.hoon new file mode 100644 index 0000000..40c6d34 --- /dev/null +++ b/librand/desk/tests/lib/unumrand.hoon @@ -0,0 +1,91 @@ + :: /tests/lib/unumrand +:::: +:: /lib/unumrand: +posit-lattice regressions at all five width doors (one +:: proving the identity-on-non-NaR-bits passthrough, one proving the +:: NaR-rejection redraw actually redraws and threads rng state correctly +:: through both draws), +posit-unit regressions at posit8/16/32 (two +:: seeds each), and one test proving the "sample at @rd via i754rand, +:: convert via /lib/unum's existing +from-rd" distribution composition +:: end to end. +:: +:: The RIGOROUS chi-square verification against librand/tools/ +:: posit_unit_check.py's exact oracle is a separate, offline check (see +:: NEXT-STEPS.md) -- not reproduced here, since embedding a 256- or +:: 65536-row expected-probability table as Hoon literals isn't practical +:: for a self-contained regression test. +:: +/+ *test, + rand, + i754rand, + unum, + unumrand +|% +:: +posit-lattice is the identity on non-NaR raw bits -- same seed, same +:: width, same result as a plain +bits:uni:rand draw, whenever that draw +:: doesn't happen to be NaR (as it isn't at seed 0, any width here). +++ test-posit-lattice-is-bits ^- tang + ;: weld + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 8)) + !>((posit-lattice:rpb:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 16)) + !>((posit-lattice:rph:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 32)) + !>((posit-lattice:rps:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 64)) + !>((posit-lattice:rpd:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>((bits:uni:rand (from-atom:seed:rand %sm64 0) 128)) + !>((posit-lattice:rpq:unumrand (from-atom:seed:rand %sm64 0))) + == +:: Seed %sm64 149's first 8-bit draw is EXACTLY NaR (0x80) -- ship-found by +:: brute-force search, not hand-picked. +posit-lattice:rpb must reject it, +:: redraw (0x9b, the second 8-bit draw from the post-first-draw state), and +:: return an rng state matching two draws consumed, not one. +++ test-posit-lattice-rejects-nar ^- tang + %+ expect-eq + !>([out=0x9b r=[%sm64 s=0x3c6e.f372.fe94.f8bf]]) + !>((posit-lattice:rpb:unumrand (from-atom:seed:rand %sm64 149))) +:: +posit-unit value regressions, two seeds per width (posit8/16/32 only -- +:: posit64/128 are out of scope, see NEXT-STEPS.md). Each value is exactly +:: encode(False, -k, u, n) for u = the corresponding k-bit +bits:uni:rand +:: draw, cross-checked against librand/tools/posit_unit_check.py's encode() +:: before being taken as ground truth here (not hand-derived). +++ test-posit-unit ^- tang + ;: weld + %+ expect-eq + !>(`@`0x37) + !>(out:(posit-unit:rpb:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x39) + !>(out:(posit-unit:rpb:unumrand (from-atom:seed:rand %sm64 1))) + %+ expect-eq + !>(`@`0x3e22) + !>(out:(posit-unit:rph:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3911) + !>(out:(posit-unit:rph:unumrand (from-atom:seed:rand %sm64 1))) + %+ expect-eq + !>(`@`0x35cf.13cd) + !>(out:(posit-unit:rps:unumrand (from-atom:seed:rand %sm64 0))) + %+ expect-eq + !>(`@`0x3bee.b8da) + !>(out:(posit-unit:rps:unumrand (from-atom:seed:rand %sm64 1))) + == +:: The distribution composition rand-spec.md 12.4 describes but this +:: library doesn't wrap: draw a standard normal deviate at @rd via +:: /lib/i754rand, then convert to posit32 via /lib/unum's EXISTING +:: +from-rd (no new /lib/unum plumbing needed, unlike fixedrand). Both +:: steps are independently tested elsewhere (i754rand's own +test-normal, +:: /lib/unum's own from-rd round-trip tests); this proves the COMPOSITION +:: against a known-good ship-verified result. +++ test-unumrand-dist-composition ^- tang + =/ r (from-atom:seed:rand %sm64 0) + =^ z r (normal:rd:dist:i754rand r) + %+ expect-eq + !>(`@`0xc0d5.045b) + !>((from-rd:rps:unum z)) +-- diff --git a/librand/tools/posit8_exhaustive.c b/librand/tools/posit8_exhaustive.c new file mode 100644 index 0000000..33f1d58 --- /dev/null +++ b/librand/tools/posit8_exhaustive.c @@ -0,0 +1,95 @@ +/* Exhaustive verification of +posit-unit at posit8: encode every one of + * the 2^32 possible k-bit numerators u (k=4n=32, n=8) via a C mirror of + * /lib/unum's +bit (the same encode() logic already validated in + * posit_unit_check.py), tally per-pattern counts, and print them. + * + * This is a from-scratch re-derivation, not a transliteration copy-paste + * of posit_unit_check.py's Python -- deliberately, so a translation bug + * in one isn't invisible to the other. Cross-validated against the + * Python encode() for a large random sample (see cross_check.py in this + * directory) before trusting the full 2^32 run. + * + * Since the k-bit dyadic construction only ever encodes nonnegative + * values (u/2^k is always in [0, 1 - 2^-k]), neg=0 always here -- see + * posit_unit_check.py's reachable_patterns() for why negative-signed + * patterns are unreachable by this construction. + * + * Build: cc -O3 -o posit8_exhaustive posit8_exhaustive.c + * Run: ./posit8_exhaustive > posit8_exhaustive_counts.txt + */ +#include +#include + +#define N 8 +#define MAXPOS ((1 << (N - 1)) - 1) /* 63 */ +#define MSK ((1 << N) - 1) /* 0xff */ + +/* Floor division/modulus for signed operands -- explicit, not relying on + * C's truncating "/" or implementation-defined behavior of ">>" on + * negative values, even though in practice most compilers do arithmetic + * shift. Matches Python's ">>"/"%" semantics on negative operands exactly. */ +static int64_t fdiv(int64_t a, int64_t b) { + int64_t q = a / b, r = a % b; + if (r != 0 && ((r < 0) != (b < 0))) q--; + return q; +} + +static int bitlen64(uint64_t a) { + int n = 0; + while (a) { n++; a >>= 1; } + return n; +} + +/* Mirrors posit_unit_check.py's encode(neg=0, e, a, n=8). */ +static unsigned encode8(int64_t e, uint64_t a) { + if (a == 0) return 0; + int lead = bitlen64(a) - 1; + int64_t x = e + (int64_t)lead; + uint64_t frac = a & ((lead < 64) ? ((1ULL << lead) - 1) : ~0ULL); + int64_t r = fdiv(x, 4); + int64_t elo = x - 4 * r; /* always in [0,4), matches Python's x - 4*r */ + + if (r >= N - 2) return (unsigned)MAXPOS; + if (r <= -(N - 1)) return 1u; + + int64_t regval, regwid; + if (r >= 0) { regval = ((1LL << (r + 1)) - 1) << 1; regwid = r + 2; } + else { regval = 1; regwid = -r + 1; } + + int64_t totw = regwid + 2 + lead; + uint64_t pay = ((uint64_t)regval << (2 + lead)) | ((uint64_t)elo << lead) | frac; + int64_t pw = N - 1; + uint64_t mag; + + if (totw <= pw) { + mag = pay << (pw - totw); + } else { + int64_t sh = totw - pw; + uint64_t keep = pay >> sh; + uint64_t guard = (pay >> (sh - 1)) & 1; + uint64_t sticky = (sh >= 2) ? ((pay & ((1ULL << (sh - 1)) - 1)) != 0) : 0; + uint64_t lsbit = keep & 1; + if (guard && (sticky || lsbit)) keep++; + if (keep > (uint64_t)MAXPOS) keep = (uint64_t)MAXPOS; + mag = keep; + } + return (unsigned)mag; +} + +int main(void) { + uint64_t hist[1 << N] = {0}; + const int64_t e = -32; /* k = 4n = 32 for posit8 */ + /* u ranges over exactly [0, 2^32) -- a 32-bit counter wrapping to 0 + * after UINT32_MAX covers this with a do/while, avoiding a 64-bit + * comparison against 1ULL<<32 on every iteration of a 4-billion- + * iteration hot loop. */ + uint32_t u = 0; + do { + hist[encode8(e, u)]++; + u++; + } while (u != 0); + + for (unsigned p = 0; p < (1u << N); p++) + if (hist[p]) printf("%u %llu\n", p, (unsigned long long)hist[p]); + return 0; +} diff --git a/librand/tools/posit_unit_check.py b/librand/tools/posit_unit_check.py new file mode 100644 index 0000000..37c37a0 --- /dev/null +++ b/librand/tools/posit_unit_check.py @@ -0,0 +1,468 @@ +#!/usr/bin/env python3 +"""Oracle for librand's /lib/unumrand +posit-unit (rand-spec.md section 12.4). + ++posit-unit draws k raw bits u and encodes the dyadic u/2**k (the g-layer +value [%p %.y -k u]) via /lib/unum's +bit (round-to-nearest-even, +saturating) to get a uniform value on [0,1) at posit precision. This script +is the independent oracle: it reuses libmath/tools/posit_check.py's own +decode/encode reference model (already cross-checked against SoftPosit) to +derive, EXACTLY -- via Fraction arithmetic, no floating point anywhere -- +the true probability every reachable posit pattern must receive under this +construction. That table is then chi-squared against empirical counts drawn +on-ship at a fixed seed. + +decode()/encode() below are a verbatim copy of libmath/tools/posit_check.py's +reference model (not imported, to keep this script standalone and runnable +from any directory) -- see that file's own header for its SoftPosit +cross-check history. Everything past that point (the expected-probability +derivation, binning, chi-square) is new, specific to +posit-unit. + +Usage: + # Inspect the exact expected-probability table: + python3 posit_unit_check.py table posit8 + python3 posit_unit_check.py table posit16 + + # Chi-square ship-drawn raw patterns (one integer per line, hex 0x.. or + # decimal) against the exact expected table, binned into equal-mass groups: + python3 posit_unit_check.py chi2 posit8 counts.txt --bins 32 + python3 posit_unit_check.py chi2 posit16 counts.txt --bins 64 + +Optional: scipy, for a p-value alongside the raw statistic. Falls back to +reporting just the statistic and degrees of freedom if scipy isn't +installed. +""" +import sys +import argparse +from fractions import Fraction + +try: + from scipy.stats import chi2 as _scipy_chi2 +except ImportError: + _scipy_chi2 = None + +WIDTHS = {'posit8': 8, 'posit16': 16, 'posit32': 32} +# k = 4n, per rand-spec.md 12.4: gives a constant 7-bit safety margin over +# the finest rounding-cell width in [0,1) at any n (that width sits at +# exponent -(4(n-2)+1), so k=4n exceeds it by exactly 4n - (4(n-2)+1) = 7). +KBITS = {'posit8': 32, 'posit16': 64, 'posit32': 128} + +# ---- verbatim reference model (see module docstring) ---- + +def decode(p, n): + msk = (1 << n) - 1; p &= msk; nar = 1 << (n - 1) + if p == 0: return ('z',) + if p == nar: return ('n',) + neg = (p >> (n - 1)) & 1 + mag = (1 << n) - p if neg else p + pw = n - 1; r0 = (mag >> (pw - 1)) & 1; k = 1 + while True: + if k == pw: + r = (k - 1) if r0 == 1 else -k + return ('p', neg, 4 * r, 1) + if (mag >> (pw - 1 - k)) & 1 == r0: k += 1; continue + break + r = (k - 1) if r0 == 1 else -k + remwid = pw - (k + 1); rem = mag & ((1 << remwid) - 1) + if remwid >= 2: elo = rem >> (remwid - 2); fw = remwid - 2 + elif remwid == 1: elo = rem << 1; fw = 0 + else: elo = 0; fw = 0 + frac = rem & ((1 << fw) - 1) + x = 4 * r + elo; a = (1 << fw) + frac + return ('p', neg, x - fw, a) + +def encode(neg, e, a, n): + if a == 0: return 0 + msk = (1 << n) - 1; maxpos = (1 << (n - 1)) - 1 + lead = a.bit_length() - 1; x = e + lead; frac = a & ((1 << lead) - 1) + r = x >> 2; elo = x - 4 * r + if r >= n - 2: return ((1 << n) - maxpos) & msk if neg else maxpos + if r <= -(n - 1): return ((1 << n) - 1) & msk if neg else 1 + if r >= 0: regval = ((1 << (r + 1)) - 1) << 1; regwid = r + 2 + else: regval = 1; regwid = -r + 1 + totw = regwid + 2 + lead + pay = (regval << (2 + lead)) | (elo << lead) | frac + pw = n - 1 + if totw <= pw: + mag = pay << (pw - totw) + else: + sh = totw - pw; keep = pay >> sh + guard = (pay >> (sh - 1)) & 1; low = pay & ((1 << (sh - 1)) - 1) + if guard and ((1 if low else 0) or (keep & 1)): keep += 1 + if keep > maxpos: keep = maxpos + mag = keep + return ((1 << n) - mag) & msk if neg else mag + +def value(p, n): + """Exact Fraction value of pattern p at width n, or None for NaR.""" + d = decode(p, n) + if d[0] == 'z': return Fraction(0) + if d[0] == 'n': return None + _, neg, e, a = d + v = Fraction(a) * (Fraction(2) ** e) + return -v if neg else v + +# ---- exact expected-probability table for +posit-unit ---- + +def reachable_patterns(n): + """Nonneg-valued patterns (0 plus positive posits) with value <= 1, + sorted by value ascending. +posit-unit can only ever produce one of + these: the input is always in [0, 1 - 2**-k], and RNE rounding of a + nonnegative value never flips sign or overshoots past 1.0's own + pattern.""" + entries = [] + for p in range(1 << n): + d = decode(p, n) + if d[0] == 'n': + continue + if d[0] == 'z': + entries.append((Fraction(0), p)) + continue + _, neg, e, a = d + if neg: + continue + v = Fraction(a) * (Fraction(2) ** e) + if v <= 1: + entries.append((v, p)) + entries.sort() + return entries + +def expected_probabilities(width_name): + """Exact Fraction probability for every reachable pattern. + + For each pattern, finds the exact range of numerators u in [0, 2**k) + that round to it, via binary search against the TRUSTED `encode` + (the same round-to-nearest-even logic /lib/unum's +bit implements) -- + not by re-deriving midpoint/tie-break math by hand. Nonneg posit + patterns are monotonic in both raw-integer and decoded-value order, so + comparing raw pattern integers directly (no re-decode) is enough. + """ + n = WIDTHS[width_name] + k = KBITS[width_name] + two_k = 1 << k + entries = reachable_patterns(n) + patterns = [p for _, p in entries] + m = len(patterns) + + def first_u_reaching(target_pattern): + lo, hi = 0, two_k + while lo < hi: + mid = (lo + hi) // 2 + out = encode(False, -k, mid, n) + if out >= target_pattern: + hi = mid + else: + lo = mid + 1 + return lo + + thresholds = [0] * (m + 1) + for i in range(1, m): + thresholds[i] = first_u_reaching(patterns[i]) + thresholds[m] = two_k + + probs = {} + for i in range(m): + count = thresholds[i + 1] - thresholds[i] + assert count >= 0, (i, patterns[i], thresholds[i], thresholds[i + 1]) + probs[patterns[i]] = Fraction(count, two_k) + assert sum(probs.values()) == 1 + return probs, entries + +# ---- binning + chi-square ---- + +def quantile_bins(entries, probs, num_bins): + """Partition the value-sorted pattern list into ~num_bins contiguous + groups of roughly equal expected mass (never splitting a pattern across + a bin boundary). Returns a list of (patterns_in_bin, total_prob).""" + bins = [] + cur_patterns = [] + cur_prob = Fraction(0) + target = Fraction(1, num_bins) + for _, pattern in entries: + cur_patterns.append(pattern) + cur_prob += probs[pattern] + if cur_prob >= target and len(bins) < num_bins - 1: + bins.append((cur_patterns, cur_prob)) + cur_patterns = [] + cur_prob = Fraction(0) + if cur_patterns: + bins.append((cur_patterns, cur_prob)) + return bins + +def chi_square_stat(counts_by_pattern, bins, n_draws): + stat = Fraction(0) + rows = [] + for pats, prob in bins: + observed = sum(counts_by_pattern.get(p, 0) for p in pats) + expected = prob * n_draws + term = (Fraction(observed) - expected) ** 2 / expected + stat += term + rows.append((observed, expected, term)) + return stat, rows + +def parse_counts(path): + """One raw pattern per line (hex 0x.. or decimal); returns pattern->count.""" + counts = {} + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + p = int(line, 0) + counts[p] = counts.get(p, 0) + 1 + return counts + +# ---- large-width (posit16+) binned harness ---- +# +# posit16 has 16,385 reachable patterns -- far too many to safely dump a +# raw per-pattern histogram out of a live dojo session (tried it: the +# terminal scrollback truncates silently, and even RAISING tmux's +# history-limit didn't help reliably; ship-crashed twice trying to embed +# a 16,385-entry pattern->bin lookup table as a single big Hoon literal). +# +# Fix: have the SHIP do the binning itself, using the SAME quantile bins +# this script already computes, so it only ever needs to print `--bins` +# (e.g. 64) numbers, not thousands. Two Hoon lessons paid for finding +# this the hard way, baked into gen_binned_hoon() below so they don't +# recur: +# 1. `~[a b c ...]` list literals infer a FIXED-SHAPE tuple type, not +# the general recursive (list @) mold -- using one as a `|-` trap's +# loop variable across recursive $(...) calls fails with a +# confusing "mint-vain" (the shape genuinely differs each +# iteration as the tuple gets shorter). Always `^- (list @)` cast +# list literals before using them this way. +# 2. Decimal literals over 3 digits NEED Hoon's dot-thousands- +# separator (`50.000`, not `50000`) -- a hard parse error otherwise. +# The bin-of gate below was verified against known boundary values +# (0,10,20,30 -> bin-of should give 0,0,0,1,1,1,2,2,2,3,3 for +# 0,5,9,10,15,19,20,25,29,30,100) before trusting a real run. + +def hoon_dec(n): + """Format a decimal integer with Hoon's required dot-thousands-separators.""" + s = str(n) + if len(s) <= 3: + return s + parts = [] + while len(s) > 3: + parts.insert(0, s[-3:]) + s = s[:-3] + parts.insert(0, s) + return '.'.join(parts) + +def bin_boundaries(entries, bins): + """The lowest pattern in each bin. Relies on +posit-unit's reachable + patterns (all nonneg-valued) sorting identically by raw integer and + by decoded value, so quantile_bins' value-sorted contiguous groups + are ALSO contiguous ranges of raw pattern integers -- confirmed for + this construction elsewhere in this project's own verification work.""" + boundaries = [min(pats) for pats, _ in bins] + assert boundaries == sorted(boundaries) + return boundaries + +def gen_binned_hoon(width_name, num_bins, count, seed): + n = WIDTHS[width_name] + door = {'posit8': 'rpb', 'posit16': 'rph', 'posit32': 'rps'}[width_name] + probs, entries = expected_probabilities(width_name) + bins = quantile_bins(entries, probs, num_bins) + boundaries = bin_boundaries(entries, bins) + lit = ' '.join(hoon_dec(b) for b in boundaries) + src = f"""\ +/+ rand, unumrand +=/ boundaries ^- (list @) ~[{lit}] +=/ bin-of + |= v=@ + ^- @ + =/ bs boundaries + =/ i 0 + |- ^- @ + ?~ bs (dec i) + ?: (gth i.bs v) + (dec i) + $(bs t.bs, i +(i)) +=/ count {hoon_dec(count)} +=/ seed 0x{seed:x} +=/ r (from-atom:seed:rand %sm64 seed) +=/ m *(map @ @) +=/ i 0 +|- ^- (map @ @) +?: =(i count) + m +=^ v r (posit-unit:{door}:unumrand r) +=/ bin (bin-of v) +=/ c (fall (~(get by m) bin) 0) +$(i +(i), m (~(put by m) bin +(c))) +""" + return src, bins + +def cmd_gen_binned(args): + src, bins = gen_binned_hoon(args.width, args.bins, args.count, args.seed) + with open(args.output, 'w') as f: + f.write(src) + print(f"wrote {args.output} ({len(bins)} bins, {args.count} draws, seed=0x{args.seed:x})") + print(f"Deploy to a ship's %base desk, |commit, then run bare (no `=face` prefix --") + print(f"that specific combination with -build-file on a raw-expression file has been") + print(f"unreliable in dojo):") + print(f" -build-file %/lib/{args.output.rsplit('/', 1)[-1].removesuffix('.hoon')}/hoon") + print(f"Capture the printed (map @ @) (bin -> count), then run:") + print(f" python3 {sys.argv[0]} chi2-binned {args.width} --bins {args.bins}") + +def parse_bin_counts(path): + """Parses a captured dojo pane dump of a printed (map @ @): lines like + '[p=N q=M]' (N=bin index, M=count), Hoon's own decimal-dot-grouped + display. Strips the dots before converting to int.""" + import re + counts = {} + with open(path) as f: + text = f.read() + for m in re.finditer(r'p=([\d.]+)\s+q=([\d.]+)', text): + bi = int(m.group(1).replace('.', '')) + c = int(m.group(2).replace('.', '')) + counts[bi] = c + return counts + +def cmd_chi2_binned(args): + probs, entries = expected_probabilities(args.width) + bins = quantile_bins(entries, probs, args.bins) + counts = parse_bin_counts(args.counts_file) + n_draws = sum(counts.values()) + stat = Fraction(0) + rows = [] + for bi, (pats, prob) in enumerate(bins): + observed = counts.get(bi, 0) + expected = prob * n_draws + term = (Fraction(observed) - expected) ** 2 / expected + stat += term + rows.append((observed, expected, term)) + dof = len(bins) - 1 + stat_f = float(stat) + print(f"{args.width}: N={n_draws} draws, {len(bins)} bins, dof={dof}") + print(f"chi-square statistic = {stat_f:.4f}") + if _scipy_chi2 is not None: + p_value = 1.0 - _scipy_chi2.cdf(stat_f, dof) + print(f"p-value = {p_value:.4f} (fail to reject H0 [uniform-on-[0,1)] at alpha=0.01 if p > 0.01)") + else: + print("(scipy not installed -- compare the statistic above to a " + f"chi-square critical value table at dof={dof} yourself)") + if args.verbose: + for i, ((pats, prob), (observed, expected, term)) in enumerate(zip(bins, rows)): + print(f" bin {i:3d}: patterns={len(pats):5d} observed={observed:6d} " + f"expected={float(expected):9.2f} term={float(term):.4f}") + +# ---- CLI ---- + +def cmd_table(args): + probs, entries = expected_probabilities(args.width) + n = WIDTHS[args.width] + print(f"{args.width}: {len(entries)} reachable patterns, k={KBITS[args.width]}") + for v, p in entries: + print(f" 0x{p:0{(n+3)//4}x} value={float(v):.10g} prob={probs[p]} ({float(probs[p]):.3e})") + +def cmd_chi2(args): + probs, entries = expected_probabilities(args.width) + bins = quantile_bins(entries, probs, args.bins) + counts = parse_counts(args.counts_file) + n_draws = sum(counts.values()) + stat, rows = chi_square_stat(counts, bins, n_draws) + dof = len(bins) - 1 + stat_f = float(stat) + print(f"{args.width}: N={n_draws} draws, {len(bins)} bins, dof={dof}") + print(f"chi-square statistic = {stat_f:.4f}") + if _scipy_chi2 is not None: + p_value = 1.0 - _scipy_chi2.cdf(stat_f, dof) + print(f"p-value = {p_value:.4f} (fail to reject H0 [uniform-on-[0,1)] at alpha=0.01 if p > 0.01)") + else: + print("(scipy not installed -- compare the statistic above to a " + f"chi-square critical value table at dof={dof} yourself)") + if args.verbose: + for i, ((pats, prob), (observed, expected, term)) in enumerate(zip(bins, rows)): + print(f" bin {i:3d}: patterns={len(pats):5d} observed={observed:6d} " + f"expected={float(expected):9.2f} term={float(term):.4f}") + +def cmd_exhaustive8(args): + """Compile and run posit8_exhaustive.c (an independent, from-scratch C + re-derivation of the encode logic, not a transliteration of this + file's own encode()), then compare its full 2**32-numerator tally + against expected_probabilities('posit8') EXACTLY -- eliminates + sampling error entirely, since it's a complete enumeration of the + construction's domain, not a statistical sample.""" + import subprocess, os, tempfile + here = os.path.dirname(os.path.abspath(__file__)) + src = os.path.join(here, 'posit8_exhaustive.c') + with tempfile.TemporaryDirectory() as td: + binf = os.path.join(td, 'posit8_exhaustive') + countsf = os.path.join(td, 'counts.txt') + print(f"compiling {src}...") + subprocess.run(['cc', '-O3', '-o', binf, src], check=True) + print("running (full 2^32 enumeration)...") + with open(countsf, 'w') as f: + subprocess.run([binf], stdout=f, check=True) + exhaustive = {} + total = 0 + with open(countsf) as f: + for line in f: + p, c = line.split() + exhaustive[int(p)] = int(c) + total += int(c) + assert total == 1 << 32, f"expected 2**32 total, got {total}" + + probs, entries = expected_probabilities('posit8') + mismatches = 0 + for pattern, frac in probs.items(): + assert (1 << 32) % frac.denominator == 0, (pattern, frac) + expected_count = frac.numerator * ((1 << 32) // frac.denominator) + got = exhaustive.get(pattern, 0) + if got != expected_count: + print(f" MISMATCH pattern=0x{pattern:02x} expected={expected_count} got={got}") + mismatches += 1 + for pattern, count in exhaustive.items(): + if pattern not in probs: + print(f" UNEXPECTED pattern=0x{pattern:02x} count={count}") + mismatches += 1 + print(f"patterns checked: {len(probs)}, total enumerated: {total} (== 2**32: {total == 1<<32})") + print("EXACT MATCH -- zero sampling error at posit8" if mismatches == 0 + else f"FAILED -- {mismatches} mismatches") + return mismatches + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + sub = ap.add_subparsers(dest='cmd', required=True) + + t = sub.add_parser('table', help='print the exact expected-probability table') + t.add_argument('width', choices=WIDTHS.keys()) + t.set_defaults(func=cmd_table) + + c = sub.add_parser('chi2', help='chi-square ship-drawn counts against the exact table') + c.add_argument('width', choices=['posit8', 'posit16']) + c.add_argument('counts_file') + c.add_argument('--bins', type=int, default=32) + c.add_argument('--verbose', action='store_true') + c.set_defaults(func=cmd_chi2) + + x = sub.add_parser('exhaustive8', + help='compile+run posit8_exhaustive.c, compare against the oracle exactly') + x.set_defaults(func=cmd_exhaustive8) + + g = sub.add_parser('gen-binned', + help='generate a .hoon harness that bins draws on-ship (for widths too ' + 'large to safely dump a raw per-pattern histogram out of dojo)') + g.add_argument('width', choices=['posit16', 'posit32']) + g.add_argument('output', help='output .hoon file path') + g.add_argument('--bins', type=int, default=64) + g.add_argument('--count', type=int, default=1_000_000) + g.add_argument('--seed', type=lambda s: int(s, 0), default=0x2001) + g.set_defaults(func=cmd_gen_binned) + + b = sub.add_parser('chi2-binned', + help='chi-square a captured dojo pane dump of gen-binned\'s printed ' + '(map @ @) bin counts against the exact table') + b.add_argument('width', choices=['posit16', 'posit32']) + b.add_argument('counts_file', help='captured pane text containing [p=bin q=count] entries') + b.add_argument('--bins', type=int, default=64) + b.add_argument('--verbose', action='store_true') + b.set_defaults(func=cmd_chi2_binned) + + args = ap.parse_args() + args.func(args) + +if __name__ == '__main__': + main() diff --git a/rand-spec.md b/rand-spec.md new file mode 100644 index 0000000..f46d6ae --- /dev/null +++ b/rand-spec.md @@ -0,0 +1,375 @@ +# SPEC: `/lib/rand` — Deterministic Random Number Generation for Urbit Numerics + +**Target repo:** `urbit/numerics`, new directory `librand/` (structured like `libmath/`), with a Saloon extension for ray-filling. +**Status:** Implemented. All nine milestones (section 13) are done — see `librand/README.md` and `librand/NEXT-STEPS.md` (including its "Deferred to v2" section for what's intentionally out of scope) for the current state. Jetting (section 11) remains a follow-up milestone. +**Dependencies:** `/lib/math` (float transcendentals), `/lib/twoc` (width-keyed modular integers), optionally `/lib/unum` and `/lib/fixed` for output adapters. Saloon layer depends on `/lib/lagoon`. + +--- + +## 0. Design principles (read before writing code) + +1. **The Hoon is the spec.** Every generator and every distribution transform must be defined such that a C jet can reproduce it **bit-for-bit**. This is the same discipline as `/lib/math`'s SoftFloat parity. Consequences: + - All float arithmetic in distribution transforms goes through the stdlib `@r*` doors or `/lib/math` kernels (which have jet-parity guarantees). No algorithm may depend on extended-precision intermediates, FMA availability, or evaluation order beyond what the Hoon states. + - All integer arithmetic is exact `@` arithmetic masked to width via `/lib/twoc` (`+twid` at width 32/64/128) or explicit `(mod . (bex w))`. Never rely on implicit truncation. + +2. **Counter-based first.** The primary generator is **Philox4x32-10** (Salmon et al., "Parallel Random Numbers: As Easy as 1, 2, 3", SC'11; the Random123 library is the reference). Rationale: + - Output is a pure function `(key, counter) -> 4 x u32`. No sequential state dependency, so a Lagoon jet can fill a ray in parallel and still match the sequential Hoon loop exactly. + - Reproducibility across ships, runtimes, and replays (Nockchain verification contexts) falls out for free. + - Stream splitting is trivial and collision-free (distinct keys = distinct streams). + +3. **Sequential generators are provided but secondary.** `+split-mix` (SplitMix64) for seed expansion and cheap draws; `+pcg` (PCG64 XSL-RR 128/64) for users who want a conventional stateful generator. **xoshiro is explicitly excluded** — nothing it offers matters here and its low-bit weaknesses are a support liability. + +4. **Not cryptographic. Say so loudly.** This library produces *deterministic pseudo-randomness for numerical work* (Monte Carlo, ML init, shuffles, simulation). It must NOT be used for key material, nonces, or anything adversarial. Entropy acquisition is Arvo's job (`eny`); cryptographic randomness is Zuse's job. The library docstring, the README, and each seeding arm that accepts `eny` must carry this warning. + +5. **State threading is explicit and total.** Every drawing arm returns `[value new-state]` over plain-data state (section 2.1 defines the full discipline: functional core, door facade, key derivation). No arm hides state mutation. Rejection loops are fine (they thread state a variable number of times) but must be documented as variable-consumption arms. Relationship to stdlib `+og`: we adopt its door *ergonomics* as a facade, not its engine (iterated SHA-256, slow, unjetted for this purpose) and not its historical misuse (persisting the core). + +6. **Match house style.** Doors and arm docs follow the existing repo conventions (see `/lib/math`): `:: +arm: type -> type` header comment, `Examples` block with dojo output, `Source` marker, `~/` jet hints with `~% %non ..part ~` registration at the core. + +--- + +## 1. Module layout + +**Revised during implementation (after milestone 6), superseding the +single-file layout originally specified below.** A single `/lib/rand` +covering sections 2-7 grew large enough, and IEEE-754-float-specific +enough in its `++uni`/`++dist` halves, that it was split to match how +`/lib/twoc`/`fixed`/`complex`/`unum` are already kept separate from +`/lib/math` in this codebase — see `librand/NEXT-STEPS.md` for the +rationale and the exact arm-by-arm split. Current layout: + +``` +librand/ + README.md + NEXT-STEPS.md + desk/sur/rand.hoon :: +$rng, +$phil, +$sm64, +$pcg64 + desk/lib/rand.hoon :: plumbing: ++philox/++split-mix/++pcg/ + :: ++seed/+step/+fork/++gen/++uni(integer + :: only)/++sample(minus ++alias) + desk/lib/i754rand.hoon :: porcelain: ++uni(float)/++alias/++dist + desk/lib/twocrand.hoon :: non-float adapter (section 12.2) + desk/lib/fixedrand.hoon :: non-float adapter (section 12.1) + desk/lib/complexrand.hoon :: non-float adapter (section 12.3) + desk/lib/unumrand.hoon :: non-float adapter (section 12.4) +saloon/desk/lib/saloon.hoon :: gains a +rand-ray core (section 7) +``` + +Original spec (single library file, one top-level core, sub-cores per +concern) — kept for historical context, no longer the actual layout: + +``` +/lib/rand +|% +++ philox :: counter-based engine (primary) +++ split-mix :: SplitMix64: seed expansion + cheap sequential generator +++ pcg :: PCG64 (XSL-RR 128/64) sequential generator +++ seed :: seeding utilities (from @, from eny, hash-mixing) +++ uni :: uniform deviates: bits, bounded ints, floats per aura +++ dist :: nonuniform distributions (normal, exp, gamma, ...) +++ sample :: shuffles, permutations, choice, reservoir, alias tables +-- +``` + +--- + +## 2. Core types + +```hoon +:: Philox key/counter state. ctr is the 128-bit counter as a single @, +:: key is the 64-bit key as a single @. Both are plain atoms; width +:: discipline is enforced by masking, never by aura tricks. ++$ phil [key=@ ctr=@] + +:: Sequential generator state. ++$ sm64 @ :: SplitMix64: 64-bit state ++$ pcg64 [state=@ inc=@] :: 128-bit state, 128-bit odd increment + +:: A generic stream: tagged union so distribution code is engine-agnostic. ++$ rng + $% [%phil p=phil] + [%sm64 s=sm64] + [%pcg p=pcg64] + == +``` + +Distribution and sampling arms take and return `rng`, calling a single internal `+step` gate, which lives at the top level of `/lib/rand` (not nested under any one engine core, since it dispatches across the `+$rng` tagged union): + +```hoon +:: +step: rng -> [u64=@ rng] (one 64-bit draw, engine-dispatched) +``` + +This keeps `+dist` and `+sample` written once, not per-engine. Resolved design (no buffering, no caching — simplicity and an unambiguous spec beat throughput at v1): for `%phil`, `+step` calls `next:philox`, which returns the full 128-bit block `out = (block key ctr)`; `+step` takes **bits `[0,64)` of that atom**, i.e. lanes `c0'` and `c1'` (the low two lanes of the little-endian repack — `out = c3'*2^96 + c2'*2^64 + c1'*2^32 + c0'`, so bits `[0,64)` = `c1'*2^32 + c0'`), and increments the counter by 1. The top 64 bits (`c2'`, `c3'`) are discarded; wasting them is acceptable at v1, and a buffered variant that reuses both halves is a documented follow-up in NEXT-STEPS. For `%sm64`, `+step` calls `next:split-mix` directly (already 64 bits). For `%pcg`, `+step` calls `next:pcg`, whose `xsl-rr` output is already a 64-bit rotr64 result. + +### 2.1 State discipline and API surface + +Three sanctioned usage patterns; the API serves all three. + +**(a) Functional core — canonical.** Every drawing arm is `[value new-rng]` over the plain-data `+$rng` noun. Call-site idiom is `=^`: + +```hoon +=^ x rng (rd:uni:rand rng) +=^ y rng (normal:dist:rand rng) +``` + +Every arm's `Examples` block demonstrates the `=^` form. Agent state stores the `+$rng` **noun**. + +**(b) Door facade — ergonomic sugar, never a storage format.** Provide `++gen`, an `+og`-style door whose sample is the `rng` noun and whose arms mirror the functional core, returning `[value _+>.$]`: + +```hoon +=/ g ~(. gen:rand (from-atom:seed:rand %phil 42)) +=^ x g rd:g +``` + +**HARD RULE, stated in the door's docstring:** the door is a call-site convenience. Persisting it in agent state pins a battery across library upgrades — the classic `og` misuse. Persist the `+$rng` data noun; reconstruct the door locally. The facade adds no capability; it is a strict wrapper over (a), so jets register on the functional arms only. + +**(c) Key derivation — JAX-style, first-class, the payoff of counter-based-first.** + +```hoon +:: +fork: [rng @] -> rng (deterministic, independent child stream) +++ fork |=([r=rng salt=@] ^-(rng)) +``` + +Per engine: `%phil` — child key = `(mix key salt)`, the two-word compression pinned down in section 4, counter reset to 0. `%pcg` — child increment derived the same way (`(mix inc salt)`), forced odd; state re-mixed. `%sm64` — the SplitMix split construction. Distinct salts give computationally independent streams. + +Instead of threading one sequential state through a computation tree, callers fork by path: per-layer keys for Maroon weight initialization, per-event keys in agents (`(fork base eny-of-event-id)`), per-cell keys in simulations. Properties sequential threading cannot provide: draws are independent of sibling evaluation order, and inserting a draw in one branch does not perturb any other branch. Document this as the *recommended* pattern for tree-structured and ML workloads; sequential threading (a) remains correct and simpler for linear code. + +Nesting rule to spec: `(fork (fork r a) b)` must differ from `(fork (fork r b) a)` and from `(fork r (cat 6 a b))` — i.e., the mix is genuinely path-sensitive. Test this. + + +--- + +## 3. Engines + +### 3.1 Philox4x32-10 (`++philox`) + +Reference: Random123. Constants (u32): + +``` +M0 = 0xD2511F53 M1 = 0xCD9E8D57 +W0 = 0x9E3779B9 W1 = 0xBB67AE85 (bumpkeys per round) +Rounds = 10 +``` + +Arms: + +- `++block |=([key=@ ctr=@] ^-(@))` — the pure Philox function. Split `ctr` into `c0 c1 c2 c3` (u32 little-endian lanes), `key` into `k0 k1`. Ten rounds of the Philox S-box: + ``` + hi0.lo0 = mulhilo32(M0, c0); hi1.lo1 = mulhilo32(M1, c2) + c0' = hi1 ^ c1 ^ k0 ; c1' = lo1 + c2' = hi0 ^ c3 ^ k1 ; c3' = lo0 + k0 += W0 ; k1 += W1 (mod 2^32) + ``` + Return the four output lanes repacked as one 128-bit atom, same lane order. + `mulhilo32` is exact: `p = (mul a b)`, `hi = (rsh [0 32] p)`, `lo = (end [0 32] p)`. +- `++next |=(p=phil ^-([out=@ p=phil]))` — `out = (block key ctr)`, new state `[key (mod +(ctr) (bex 128))]`. +- **Test vectors:** the spec MUST embed the Random123 known-answer test: `philox4x32-10(ctr = {0,0,0,0}, key = {0,0})`, `ctr = key = all 0xffffffff...`, and the pi-digits vector from the Random123 `kat_vectors` file. These go in `/tests/lib/rand.hoon` and are non-negotiable — an engine that fails KAT is wrong, full stop. + +### 3.2 SplitMix64 (`++split-mix`) + +Reference: Steele, Lea, Flood (OOPSLA'14); Vigna's public-domain C. + +``` +next(s): s += 0x9E3779B97F4A7C15 (mod 2^64) + z = s + z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9 (mod 2^64) + z = (z ^ (z >> 27)) * 0x94D049BB133111EB (mod 2^64) + return z ^ (z >> 31) +``` + +Arms: `++next |=(s=sm64 [out=@ s=sm64])`, plus `++split` returning two decorrelated child states (gamma-stepping per the paper is overkill; two `+next` draws as child seeds is adequate and simpler — document the choice). + +KAT: first 5 outputs from seed `0` and seed `0xDEADBEEF` against the reference C. + +### 3.3 PCG64 XSL-RR 128/64 (`++pcg`) + +Reference: O'Neill 2014, pcg-random.org. LCG multiplier `0x2360ED051FC65DA44385DF649FCCF645` (128-bit), user-chosen odd increment (stream id). Output: `xsl-rr`: `rot = state >> 122`; `xored = (state >> 64) ^ (state & mask64)`; output = rotr64(xored, rot). + +**Correction (implementation phase, verified against `pcg-c`'s `include/pcg_variants.h` — `pcg_setseq_128_xsl_rr_64_random_r` — and independently against NumPy's vendored `pcg64.orig.h`, both consistent): an earlier draft of this spec had the operation order backwards.** The real reference does `pcg_setseq_128_step_r(rng); return pcg_output_xsl_rr_128_64(rng->state);` — i.e. **advance the state FIRST, then compute the output permutation from the NEW (already-advanced) state.** (The draft text below, retained struck through for the record, claimed the opposite: "State update AFTER output extraction (i.e., output the *current* state's permutation, then advance)" — that is wrong; do not implement it that way.) Per this spec's own closing rule (section 13): the reference implementation wins over the spec text. This is exactly the classic off-by-one-draw divergence the original text worried about, just resolved in the other direction. + +Include `++jump` (advance by `2^64` steps via LCG skip-ahead: standard `O(log n)` modular matrix trick on `(a, c)`) for stream partitioning. + +KAT: `pcg64_random_r` demo output from the reference distribution, seed `42, 54`. + +--- + +## 4. Seeding (`++seed`) + +- `++from-atom |=(a=@ ^-(rng))` — hash `a` through SplitMix64 to fill whichever engine (parameterize: `|=([eng=?(%phil %sm64 %pcg) a=@] ...)`). Philox: key = first SplitMix output, ctr = 0. PCG: state/inc from two outputs, inc forced odd via `(con inc 1)`. +- `++from-eny |=(eny=@uvJ ^-(rng))` — same pipeline, over Arvo entropy. `eny` is 512 bits; fold it to 64 bits first by chunking into eight 64-bit words (low to high) and folding sequentially through `+mix` (`acc = 0; acc = (mix acc chunk)` for each chunk in order) before handing the result to `+from-atom`'s pipeline. Carries the not-for-crypto warning. +- `++mix |=([a=@ b=@] @)` — combine two 64-bit seed words (e.g., ship + per-agent salt) into one. **Fully specified as a two-word compression, not a single finalizer call on the raw 128-bit concatenation** (the finalizer below is a function of one 64-bit word; a 128-bit input must be folded, not fed in directly): + ``` + w0 = a mod 2^64 ; w1 = b mod 2^64 + z0 = finalize(w0) :: the z-mixing tail of +next (3.2): z=z0; z=(z^(z>>30))*C1; z=(z^(z>>27))*C2; return z^(z>>31) — no `s += golden` step, that increment is +next's state-advance, not part of the finalizer + return finalize(z0 XOR w1) + ``` + Order-sensitive by construction: `(mix a b)` differs from `(mix b a)` in general, since `a` is absorbed through one extra `finalize` application relative to `b`. This is what makes `(fork (fork r a) b)` differ from `(fork (fork r b) a)` (section 2.1c) — each nested `fork` call adds another `finalize` layer at a different depth. Deterministic, no rng threading (pure `@ -> @`). For inputs wider than 64 bits (e.g. `eny`), chunk into 64-bit words and fold sequentially as shown above for `+from-eny` — do not truncate. + **Naming note (implementation footgun):** this arm's name collides with the Hoon standard library's `++mix` (bitwise XOR). Inside `++seed`, an unqualified `(mix a b)` call resolves to this shadowing local arm; any internal use of stdlib XOR inside `++seed` must reference it explicitly (e.g. `^mix`). Callers outside `++seed` are unaffected. Note this in the arm's doc comment so nobody is bitten by it during implementation. + +No implicit global seeding. There is no "default rng"; callers own their state. + +--- + +## 5. Uniform deviates (`++uni`) + +All arms `|=([r=rng ...] [out new-rng])` unless noted. + +### 5.1 Integers + +- `++bits |=([r=rng n=@] [@ rng])` — `n` uniform bits, `n <= 64` per draw, looping for larger `n` (assemble little-endian, low draw first — specify this). +- `++below |=([r=rng n=@] [@ rng])` — uniform in `[0, n)`, **unbiased**, via Lemire's multiply-shift with rejection (Lemire 2019, "Fast Random Integer Generation in an Interval"): + ``` + draw x: u64; m = x * n (exact @ arithmetic, 128-bit product) + l = m mod 2^64 + if l < n: t = (2^64 - n) mod n; while l < t: redraw x, recompute + return m >> 64 + ``` + Crash (`?>`) on `n = 0`. This is the primitive everything else (shuffle, choice) uses. Do NOT ship modulo-bias `(mod x n)` anywhere, including tests. +- `++between |=([r=rng a=@s b=@s] [@s rng])` — inclusive signed range via `+below` on the width, offset. Uses `si` arithmetic; crash on `a > b`. + +### 5.2 Floats + +One arm per aura, matching the repo's four-precision pattern: + +- `++rs |=(r=rng [@rs rng])` — uniform in `[0,1)`: draw 24 bits, `out = (mul:rs (sun:rs bits) 2^-24)` — implemented as the exact bit construction `(bits << 0) * 0x1p-24`, i.e. build the float by `~(sun rs %n)` then multiply by the constant `0x3380.0000` (`2^-24`). Every representable output is a multiple of `2^-24`; distribution is exactly uniform over that lattice. Document that `0` is possible and `1` is not. +- `++rd` — 53 bits, times `2^-53` (`0x3CA0.0000.0000.0000`). +- `++rh` — 11 bits, times `2^-11`. +- `++rq` — 113 bits, times `2^-113`. +- `++rs-oo` / `++rd-oo` (open-open `(0,1)`): same but `(bits + 0.5) * 2^-w` via drawing `w` bits (the SAME bit count as the closed-open variant, not `w-1` — drawing `w-1` bits confines the output to `(0, 0.5)`, not `(0,1)`; this was an error in an earlier draft) and setting the low half-ulp: implement as `((2*bits + 1) * 2^-(w+1))` with exact integer construction. With `bits` ranging over `[0, 2^w)`, `2*bits+1` ranges over the odd integers in `[1, 2^(w+1))`, so the result spans `(0,1)` on a lattice of spacing `2^-w`, symmetric about 0.5, never touching either endpoint. Needed by Box-Muller/log-based transforms so `log(0)` never fires. +- Posit, fixed, twoc, and complex outputs: section 12. Posit spacing is non-uniform, so posit uniformity is a distinct design with two semantics — resolved there, not deferred. + +--- + +## 6. Distributions (`++dist`) + +All doors keyed the way `/lib/math` doors are where a rounding mode matters; internally force `%n` like the math kernels do and say so. Each arm exists at `@rd` (reference) and `@rs` (routed through the same algorithm at single precision). `@rh`/`@rq` deferred to NEXT-STEPS. + +Algorithm selections — chosen for jet-parity determinism (comparison-and-arithmetic only, or transcendentals that route through `/lib/math`'s bit-specified kernels): + +| Arm | Algorithm | Reference | Notes | +|---|---|---|---| +| `++normal` | **Polar Marsaglia** (Box-Muller polar variant) | Marsaglia & Bray 1964 | Rejection on `s >= 1` or `s == 0`; uses `+log`/`+sqt` from `/lib/math`. Returns ONE deviate, discards the pair-mate at v1 (caching the mate makes state opaque; document the waste). Ziggurat is a NEXT-STEPS optimization — it needs embedded tables that must be spec'd exactly, don't let Sonnet improvise them at v1. | +| `++normal-mv` | mean/sigma wrapper | — | `mu + sigma * z`. Crash on `sigma < 0`. | +| `++expon` | inversion: `-log(u)` over `(0,1)` draw | — | rate parameter `lambda`: divide. Crash `lambda <= 0`. | +| `++gamma` | **Marsaglia–Tsang** squeeze | Marsaglia & Tsang 2000 | `alpha >= 1` direct; `alpha < 1` via boost `gamma(alpha+1) * u^(1/alpha)` (uses `+pow` from math). Crash `alpha <= 0`. | +| `++beta` | two gammas: `x/(x+y)` | — | | +| `++chi2` | `gamma(k/2, 2)` | — | | +| `++student-t` | normal / sqrt(chi2/k) | — | | +| `++poisson` | Knuth product for `lambda < 10`; **PTRS** (Hörmann 1993) transformed rejection for `lambda >= 10` | Hörmann, "The transformed rejection method for generating Poisson random variables" | Knuth-only is O(lambda) and will bite someone; spec both branches and the exact switch point. | +| `++binomial` | inversion for `n*min(p,1-p) < 30`; **BTPE** deferred | Kachitvichyanukul & Schmeiser 1988 | v1 ships inversion + a documented crash (`?>`) above the threshold rather than a slow or wrong large-n path. BTPE in NEXT-STEPS. | +| `++geometric` | `ceil(log(u)/log(1-p))` | — | edge: `p = 1` returns 1; crash `p <= 0` or `p > 1`. | +| `++bernoulli` | `u < p` | — | returns `?`. | +| `++dirichlet` | k gammas, normalized | — | list in, list out. | +| `++categorical` | via `++alias` (section 6.1) | Vose 1991 | | + +### 6.1 Alias method (`++alias` in `++sample`) + +Vose's O(n) construction, O(1) draw. Table type: + +```hoon ++$ alias-table [n=@ prob=(list @rd) alias=(list @ud)] +``` + +Construction is pure (no rng); drawing takes `[t=alias-table r=rng]`. Weights in as `(list @rd)`, non-negative, at least one positive (crash otherwise). Normalization inside construction. Spec the small/large worklist algorithm precisely — it's the classic place for off-by-one and float-compare bugs; require the construction to use `@rd` `%n` arithmetic only. + +--- + +## 7. Sampling utilities (`++sample`) + +- `++shuffle |=([r=rng l=(list)] [(list) rng])` — **Fisher–Yates**, iterating from the END down (`i` from `n-1` to `1`, swap with `(below r +(i))`). Implement over a `(map @ud *)` or flopped-list index structure to avoid O(n²) `snag`/`oust`; O(n log n) via map is fine, note the tradeoff. Wet gate (`|*`) so element type is preserved. +- `++permutation |=([r=rng n=@] [(list @ud) rng])` — shuffle of `(gulf 0 (dec n))`. +- `++choice |=([r=rng l=(list)] [* rng])` — one element, uniform. Crash on empty. +- `++choices` — n with replacement. +- `++sample-n` — n WITHOUT replacement: partial Fisher–Yates (first n of a permutation), not repeated rejection. +- `++reservoir |=([r=rng n=@ l=(list)] ...)` — Algorithm R (Vitter). Included because agents streaming events want it. + +--- + +## 8. Saloon layer (`+rand-ray`) + +In `/lib/saloon`, a core taking a Lagoon `meta` and an `rng`, returning `[ray rng]`: + +- `++fill-uniform |=([=meta r=rng] [ray rng])` — uniform `[0,1)` at the meta's float aura/bloq. Iterates the counter per element in **row-major (C) order** — the order must be specified because the jet will parallelize and must land identical bits per element. For `%phil`, element `i` uses counter `ctr0 + i` and the post-state is `ctr0 + n`; this is the jet-parallelism payoff and the reason Philox is primary. For sequential engines, elements are drawn in row-major sequence. +- `++fill-normal`, `++fill-expon` — same pattern. **Constraint:** rejection-based transforms (polar normal) consume variable draws, which breaks per-element counter assignment. Resolve: per-element sub-counter space — element `i` owns counters `ctr0 + i*2^32 ..`, rejection walks within its window (window exhaustion is astronomically improbable; crash if it happens). Post-state after filling `n` elements is `ctr0 + n*2^32` — i.e. the base counter advances past the *entire* window block for every element, not just the sub-counters actually consumed within each element's rejection loop. This keeps the post-state a pure function of `n`, independent of how many draws each element's rejection happened to need, so replay and post-state composition (e.g. chaining another `+fill-*` call) stay simple. Spec this exactly; it is the one genuinely novel design element in the library. +- Integer rays: `++fill-below` for `%uint` rays via Lemire. + +Aura/type support at v1: `%r32`/`%r64` rays. Posit rays deferred. + +--- + +## 9. Errors and edge cases (uniform policy) + +- Domain violations crash via `?>` with `~|` tags (`%rand-empty-list`, `%rand-bad-prob`, `%rand-zero-modulus`), matching `/lib/fixed`'s style. No NaN-returning "soft" errors for parameter mistakes — parameters are programmer input, not data. +- Float-valued *data* edge cases (e.g., a weight list containing NaN) crash rather than propagate. +- All arms total over their asserted domain: no infinite loops. Rejection loops must have a documented, astronomically-bounded expected iteration count; the Saloon counter-window variant gets an explicit crash on window exhaustion. + +## 10. Testing + +`/tests/lib/rand.hoon`, runnable with `-test %/tests/lib/rand ~`: + +1. **Known-answer tests** (blocking): Philox KAT vectors (Random123), SplitMix64 first-N, PCG64 demo output, Lemire `+below` spot values against a reference trace. +2. **Determinism/replay**: same seed → identical sequence; `+fill-uniform` sequential Hoon equals per-element counter formula. +3. **Unbiasedness**: `+below` over small `n` (e.g., 3, 6, 7) for 10k draws — chi-square statistic below a fixed threshold (hardcode the threshold and the expected counts; this is a smoke test, not TestU01). +4. **Distribution moments**: sample mean/variance of normal, exponential, gamma within `+is-close:rd` tolerances at 50k draws with a FIXED seed (so it's a regression test, not a flaky statistical test — record expected values as constants). +5. **Shuffle**: permutation property (same multiset), and uniformity smoke test over all 6 permutations of a 3-list. +6. **Offline** (README instruction, not Hoon): pipe jet output into PractRand/dieharder. Hoon-side statistical testing beyond smoke tests is out of scope. + +## 11. Jetting plan (follow-up milestone, spec'd now) + +- Register under the existing `%non` jet chapter pattern. +- Jets for: `+block` (Philox), `+next` (SplitMix, PCG), `+below`, and the Saloon `+fill-*` arms (the ones that matter for throughput; the fill jets use per-element counters and OpenMP/threads freely because the spec fixed the bit-exact per-element mapping). +- Distribution transforms need no dedicated jets at first: they decompose into engine draws + `/lib/math` calls, both already jetted. +- Parity harness: like `libmath/tools/rq_check.c`, a `librand/tools/rand_check.c` that runs the C reference (Random123, PCG reference, Vigna's SplitMix) against dojo output over a seed sweep. + +## 12. Non-float output types (`++uni` and `++dist` adapters) + +Adapters over the same engines; no new engine work. Ordered by difficulty. + +### 12.1 Fixed-point (`/lib/fixed`) + +The fixed lattice is uniform, so uniform generation is **exact by construction** — stronger than the float case, no rounding argument required. + +- `++fixed |=([r=rng p=prec] [@ rng])` — full-range uniform: draw `N = (wid p)` raw bits, return as the two's-complement pattern. +- `++fixed-unit |=([r=rng p=prec] [@ rng])` — uniform `[0,1)`: draw `b.p` bits as the fraction field, integer field zero. +- `++fixed-between` — signed range via the twoc adapter (12.2) at width N, then reinterpret. +- Distributions: sample at `@rd`, quantize round-to-nearest (extend `/lib/fixed` with a `+from-rd` mirroring its `+from-rs` if absent). Single quantization step; document it. + +### 12.2 Twoc (`/lib/twoc`) + +- `++twoc-full |=([r=rng w=@] [@ rng])` — w raw bits, interpreted as the width-w two's-complement pattern. Done. +- `++twoc-between |=([r=rng w=@ a=@ b=@] [@ rng])` — a,b as width-w patterns, `a <= b` in twoc order (crash otherwise). Bias to unsigned (subtract a via twoc `+sub`, reinterpret as unsigned offset), Lemire `+below` on span `b - a + 1`, bias back via twoc `+add`. +- **JET WARNING (put this in a source comment):** the span of a width-64 range needs 65 bits. Hoon `@` doesn't care; a C jet using `uint64_t` will overflow. The jet must special-case span = 2^w or compute the span in 128-bit. +- Since `+twid` is keyed on arbitrary widths, this adapter is the primitive `++fixed-between` delegates to. + +### 12.3 Complex (`/lib/complex`) + +Component-wise over the packed representation; one arm set per width door (`cs`/`cd` first, matching that library's ship order). + +- `++cuniform` — uniform over the unit square: two `[0,1)` component draws, `+pak`. +- `++normal-parts` — iid N(0,1) real and imaginary components. +- `++cnormal` — standard circularly-symmetric complex normal CN(0,1): components N(0, 1/2), i.e. `normal-parts` scaled by `invsqt2`. **Both exist and both are named unambiguously** — "complex Gaussian" means different things to signal-processing and statistics users, and conflating them is a guaranteed bug report. +- `++on-circle` — uniform on the unit circle: `theta = tau * u`, components via `/lib/math` `+cos`/`+sin` (bit-specified kernels, jet parity holds). +- `++in-disk` — uniform on the unit disk: rejection from the square (variable-consumption arm, document per section 5 policy). + +### 12.4 Unum / posits (`/lib/unum`) + +Two inequivalent uniform semantics. Ship both, named so they cannot be confused: + +- `++posit-lattice |=(r=rng [@ rng])` per width door — uniform over **bit patterns** excluding NaR: draw n bits, redraw on the NaR pattern (single-pattern rejection, expected 1 + 2^-n draws). The induced *value* distribution is approximately log-uniform (posits taper: dense near ±1, sparse at extremes). This is a fuzzing/property-testing primitive, not a statistics primitive; the docstring says so. +- `++posit-unit |=(r=rng [@ rng])` — uniform **value** on `[0,1)`, exactly: draw k raw bits `u`, form the g-layer dyadic `[%p %.y (dif:si --0 (sun:si k)) u]` (`%z` when u = 0), encode via `+bit` — one RNE round. Correctness condition: every posit rounding boundary in `[0,1)` at width n is dyadic with exponent >= -(4(n-2)+1) (minpos = 2^-4(n-2)), so any k above that bound makes the output *exactly* the round-to-nearest image of uniform measure — not approximately. Fix **k = 32 (posit8), 64 (posit16), 128 (posit32)**. No float intermediary, no double rounding; this construction exists only because `@` is arbitrary-precision and `+bit` is a single rounding step. + - Documented standard-conformant wrinkle: posits never round nonzero to zero (`+bit` saturates underflow to minpos, per the 2022 standard), so P(zero) = 2^-k exactly and the mass below ~minpos lands on minpos. Not a bug; write the comment so nobody "fixes" it. +- Distributions: sample at `@rd`, convert via the existing `+from-rd`. Double rounding is harmless at p8/16/32 (float64's 52 fraction bits strictly dominate posit32's max 27-bit fraction) — conveniently the same widths where `/lib/unum`'s transcendentals are verified. **p64/p128 inherit that file's existing caveat verbatim: treat as unverified until the oracle sweep is extended.** +- Tests: verify `++posit-unit` at posit8 against an mpmath oracle that computes each posit8 value's exact rounding-interval measure; chi-square the empirical counts at a fixed seed. Lives alongside `unum_cheb_check.py` in `librand/tools/`. +- NEXT-STEPS entry: quire-accumulated Monte Carlo — `+fdp` makes sample *sums* singly-rounded, a capability hardware floats don't have. Design note only; out of scope for v1. + +## 13. Milestones for Sonnet + +1. `++split-mix` + `++seed` + KATs. (Small; validates the test harness.) +2. `++philox` + KATs. `++fork` across all engines + path-sensitivity tests. `++gen` facade (thin; lands here so later milestones can use it in Examples blocks). +3. `++uni` complete (bits/below/between/floats) + bias tests. +4. `++pcg` + jump. +5. `++dist` at `@rd`: normal, expon, gamma, beta, bernoulli, geometric. Moment tests. +6. `++sample`: shuffle/permutation/choice/alias. `++dist` categorical/poisson/binomial/dirichlet, `@rs` variants. +7. Non-float adapters (section 12): twoc + fixed (one sitting — twoc is the primitive), then complex, then posit lattice/unit + mpmath oracle test. +8. Saloon `+rand-ray` with the counter-window design + replay-equality tests. +9. README + NEXT-STEPS (ziggurat, BTPE, buffered Philox, posit rays, quire Monte Carlo note, `@rh`/`@rq` distributions). + +Each milestone lands with tests green before the next starts. If any reference constant in this spec conflicts with the cited reference implementation, **the reference implementation wins** — flag the discrepancy rather than silently choosing. diff --git a/saloon/README.md b/saloon/README.md index 69a98d0..ace6ef4 100644 --- a/saloon/README.md +++ b/saloon/README.md @@ -47,6 +47,22 @@ Linear algebra (over Lagoon arrays): Set rounding mode and tolerance with `++sake` before calling `++eig` (the bare `++sa` default `rtol` is unusable); `rtol`'s width must match the component. +Random array filling (`rand-spec.md` section 8, in `/Users/neal/urbit/numerics/librand/`): + +- `++fill-uniform`, `++fill-normal`, `++fill-expon` — fill an `%i754` ray + (bloq 5/6, `@rs`/`@rd` only) element-by-element in row-major order from + an `/lib/rand` engine. Philox elements get per-element counter + treatment so a future jet can parallelize across elements and land + identical bits regardless of thread scheduling: `++fill-uniform`'s + single non-rejecting draw assigns element `i` counter `ctr0+i` + directly; the rejection-based `++fill-normal`/`++fill-expon` instead + give each element a `ctr0 + i*2^32` counter *window* to walk freely + within, with an explicit crash if one element's rejection loop ever + exhausts its own window. +- `++fill-below` — fill a `%uint` ray via Lemire's unbiased method, same + per-element windowing. +- Posit (`%unum`) rays are deferred (see librand's `NEXT-STEPS.md`). + ## References - Milton Abramowitz & Irene Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_. 1964–2010. diff --git a/saloon/desk/lib/saloon.hoon b/saloon/desk/lib/saloon.hoon index 69d3c31..da65c85 100644 --- a/saloon/desk/lib/saloon.hoon +++ b/saloon/desk/lib/saloon.hoon @@ -9,7 +9,9 @@ /- ls=lagoon /+ *lagoon, math, - complex + complex, + rand, + i754rand :: :: :::: ++sa :: (2v) vector/matrix ops ~% %saloon ..part ~ @@ -399,59 +401,59 @@ ?+ bloq !! %7 ?- fun - %neg ~(neg rq:math [rnd rtol]) - %factorial ~(factorial rq:math [rnd rtol]) - %exp ~(exp rq:math [rnd rtol]) - %sin ~(sin rq:math [rnd rtol]) - %cos ~(cos rq:math [rnd rtol]) - %tan ~(tan rq:math [rnd rtol]) - %log ~(log rq:math [rnd rtol]) - %log-10 ~(log-10 rq:math [rnd rtol]) - %log-2 ~(log-2 rq:math [rnd rtol]) - %sqrt ~(sqrt rq:math [rnd rtol]) - %cbrt ~(cbrt rq:math [rnd rtol]) + %neg ~(neg rq:math [rnd rtol `@rq`0]) + %factorial ~(factorial rq:math [rnd rtol `@rq`0]) + %exp ~(exp rq:math [rnd rtol `@rq`0]) + %sin ~(sin rq:math [rnd rtol `@rq`0]) + %cos ~(cos rq:math [rnd rtol `@rq`0]) + %tan ~(tan rq:math [rnd rtol `@rq`0]) + %log ~(log rq:math [rnd rtol `@rq`0]) + %log-10 ~(log-10 rq:math [rnd rtol `@rq`0]) + %log-2 ~(log-2 rq:math [rnd rtol `@rq`0]) + %sqrt ~(sqrt rq:math [rnd rtol `@rq`0]) + %cbrt ~(cbrt rq:math [rnd rtol `@rq`0]) == :: fun %6 ?- fun - %neg ~(neg rd:math [rnd rtol]) - %factorial ~(factorial rd:math [rnd rtol]) - %exp ~(exp rd:math [rnd rtol]) - %sin ~(sin rd:math [rnd rtol]) - %cos ~(cos rd:math [rnd rtol]) - %tan ~(tan rd:math [rnd rtol]) - %log ~(log rd:math [rnd rtol]) - %log-10 ~(log-10 rd:math [rnd rtol]) - %log-2 ~(log-2 rd:math [rnd rtol]) - %sqrt ~(sqrt rd:math [rnd rtol]) - %cbrt ~(cbrt rd:math [rnd rtol]) + %neg ~(neg rd:math [rnd rtol `@rd`0]) + %factorial ~(factorial rd:math [rnd rtol `@rd`0]) + %exp ~(exp rd:math [rnd rtol `@rd`0]) + %sin ~(sin rd:math [rnd rtol `@rd`0]) + %cos ~(cos rd:math [rnd rtol `@rd`0]) + %tan ~(tan rd:math [rnd rtol `@rd`0]) + %log ~(log rd:math [rnd rtol `@rd`0]) + %log-10 ~(log-10 rd:math [rnd rtol `@rd`0]) + %log-2 ~(log-2 rd:math [rnd rtol `@rd`0]) + %sqrt ~(sqrt rd:math [rnd rtol `@rd`0]) + %cbrt ~(cbrt rd:math [rnd rtol `@rd`0]) == :: fun %5 ?- fun - %neg ~(neg rs:math [rnd rtol]) - %factorial ~(factorial rs:math [rnd rtol]) - %exp ~(exp rs:math [rnd rtol]) - %sin ~(sin rs:math [rnd rtol]) - %cos ~(cos rs:math [rnd rtol]) - %tan ~(tan rs:math [rnd rtol]) - %log ~(log rs:math [rnd rtol]) - %log-10 ~(log-10 rs:math [rnd rtol]) - %log-2 ~(log-2 rs:math [rnd rtol]) - %sqrt ~(sqrt rs:math [rnd rtol]) - %cbrt ~(cbrt rs:math [rnd rtol]) + %neg ~(neg rs:math [rnd rtol `@rs`0]) + %factorial ~(factorial rs:math [rnd rtol `@rs`0]) + %exp ~(exp rs:math [rnd rtol `@rs`0]) + %sin ~(sin rs:math [rnd rtol `@rs`0]) + %cos ~(cos rs:math [rnd rtol `@rs`0]) + %tan ~(tan rs:math [rnd rtol `@rs`0]) + %log ~(log rs:math [rnd rtol `@rs`0]) + %log-10 ~(log-10 rs:math [rnd rtol `@rs`0]) + %log-2 ~(log-2 rs:math [rnd rtol `@rs`0]) + %sqrt ~(sqrt rs:math [rnd rtol `@rs`0]) + %cbrt ~(cbrt rs:math [rnd rtol `@rs`0]) == :: fun %4 ?- fun - %neg ~(neg rh:math [rnd rtol]) - %factorial ~(factorial rh:math [rnd rtol]) - %exp ~(exp rh:math [rnd rtol]) - %sin ~(sin rh:math [rnd rtol]) - %cos ~(cos rh:math [rnd rtol]) - %tan ~(tan rh:math [rnd rtol]) - %log ~(log rh:math [rnd rtol]) - %log-10 ~(log-10 rh:math [rnd rtol]) - %log-2 ~(log-2 rh:math [rnd rtol]) - %sqrt ~(sqrt rh:math [rnd rtol]) - %cbrt ~(cbrt rh:math [rnd rtol]) + %neg ~(neg rh:math [rnd rtol `@rh`0]) + %factorial ~(factorial rh:math [rnd rtol `@rh`0]) + %exp ~(exp rh:math [rnd rtol `@rh`0]) + %sin ~(sin rh:math [rnd rtol `@rh`0]) + %cos ~(cos rh:math [rnd rtol `@rh`0]) + %tan ~(tan rh:math [rnd rtol `@rh`0]) + %log ~(log rh:math [rnd rtol `@rh`0]) + %log-10 ~(log-10 rh:math [rnd rtol `@rh`0]) + %log-2 ~(log-2 rh:math [rnd rtol `@rh`0]) + %sqrt ~(sqrt rh:math [rnd rtol `@rh`0]) + %cbrt ~(cbrt rh:math [rnd rtol `@rh`0]) == :: fun == :: bloq :: @@ -496,23 +498,23 @@ ?+ bloq.meta !! %7 ?- fun - %pow-n ~(pow-n rq:math [rnd rtol]) - %pow ~(pow rq:math [rnd rtol]) + %pow-n ~(pow-n rq:math [rnd rtol `@rq`0]) + %pow ~(pow rq:math [rnd rtol `@rq`0]) == :: fun %6 ?- fun - %pow-n ~(pow-n rd:math [rnd rtol]) - %pow ~(pow rd:math [rnd rtol]) + %pow-n ~(pow-n rd:math [rnd rtol `@rd`0]) + %pow ~(pow rd:math [rnd rtol `@rd`0]) == :: fun %5 ?- fun - %pow-n ~(pow-n rs:math [rnd rtol]) - %pow ~(pow rs:math [rnd rtol]) + %pow-n ~(pow-n rs:math [rnd rtol `@rs`0]) + %pow ~(pow rs:math [rnd rtol `@rs`0]) == :: fun %4 ?- fun - %pow-n ~(pow-n rh:math [rnd rtol]) - %pow ~(pow rh:math [rnd rtol]) + %pow-n ~(pow-n rh:math [rnd rtol `@rh`0]) + %pow ~(pow rh:math [rnd rtol `@rh`0]) == :: fun == :: bloq :: posits (/lib/unum): bloq 3/4/5/6. @@ -544,13 +546,13 @@ :: 6=@rd, 7=@rq). Each takes the bloq as its first argument and operates on :: raw component atoms. :: - ++ fadd |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(add rh:math [rnd rtol]) x y) ?:(=(5 b) (~(add rs:math [rnd rtol]) x y) ?:(=(6 b) (~(add rd:math [rnd rtol]) x y) (~(add rq:math [rnd rtol]) x y)))))) - ++ fsub |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(sub rh:math [rnd rtol]) x y) ?:(=(5 b) (~(sub rs:math [rnd rtol]) x y) ?:(=(6 b) (~(sub rd:math [rnd rtol]) x y) (~(sub rq:math [rnd rtol]) x y)))))) - ++ fmul |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(mul rh:math [rnd rtol]) x y) ?:(=(5 b) (~(mul rs:math [rnd rtol]) x y) ?:(=(6 b) (~(mul rd:math [rnd rtol]) x y) (~(mul rq:math [rnd rtol]) x y)))))) - ++ fdiv |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(div rh:math [rnd rtol]) x y) ?:(=(5 b) (~(div rs:math [rnd rtol]) x y) ?:(=(6 b) (~(div rd:math [rnd rtol]) x y) (~(div rq:math [rnd rtol]) x y)))))) - ++ fabs |=([b=@ x=@] ^-(@ ?:(=(4 b) (~(abs rh:math [rnd rtol]) x) ?:(=(5 b) (~(abs rs:math [rnd rtol]) x) ?:(=(6 b) (~(abs rd:math [rnd rtol]) x) (~(abs rq:math [rnd rtol]) x)))))) - ++ fgte |=([b=@ x=@ y=@] ^-(? ?:(=(4 b) (~(gte rh:math [rnd rtol]) x y) ?:(=(5 b) (~(gte rs:math [rnd rtol]) x y) ?:(=(6 b) (~(gte rd:math [rnd rtol]) x y) (~(gte rq:math [rnd rtol]) x y)))))) - ++ flte |=([b=@ x=@ y=@] ^-(? ?:(=(4 b) (~(lte rh:math [rnd rtol]) x y) ?:(=(5 b) (~(lte rs:math [rnd rtol]) x y) ?:(=(6 b) (~(lte rd:math [rnd rtol]) x y) (~(lte rq:math [rnd rtol]) x y)))))) + ++ fadd |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(add rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(add rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(add rd:math [rnd rtol `@rd`0]) x y) (~(add rq:math [rnd rtol `@rq`0]) x y)))))) + ++ fsub |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(sub rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(sub rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(sub rd:math [rnd rtol `@rd`0]) x y) (~(sub rq:math [rnd rtol `@rq`0]) x y)))))) + ++ fmul |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(mul rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(mul rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(mul rd:math [rnd rtol `@rd`0]) x y) (~(mul rq:math [rnd rtol `@rq`0]) x y)))))) + ++ fdiv |=([b=@ x=@ y=@] ^-(@ ?:(=(4 b) (~(div rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(div rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(div rd:math [rnd rtol `@rd`0]) x y) (~(div rq:math [rnd rtol `@rq`0]) x y)))))) + ++ fabs |=([b=@ x=@] ^-(@ ?:(=(4 b) (~(abs rh:math [rnd rtol `@rh`0]) x) ?:(=(5 b) (~(abs rs:math [rnd rtol `@rs`0]) x) ?:(=(6 b) (~(abs rd:math [rnd rtol `@rd`0]) x) (~(abs rq:math [rnd rtol `@rq`0]) x)))))) + ++ fgte |=([b=@ x=@ y=@] ^-(? ?:(=(4 b) (~(gte rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(gte rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(gte rd:math [rnd rtol `@rd`0]) x y) (~(gte rq:math [rnd rtol `@rq`0]) x y)))))) + ++ flte |=([b=@ x=@ y=@] ^-(? ?:(=(4 b) (~(lte rh:math [rnd rtol `@rh`0]) x y) ?:(=(5 b) (~(lte rs:math [rnd rtol `@rs`0]) x y) ?:(=(6 b) (~(lte rd:math [rnd rtol `@rd`0]) x y) (~(lte rq:math [rnd rtol `@rq`0]) x y)))))) ++ f0 |=(b=@ ^-(@ ?:(=(4 b) .~~0 ?:(=(5 b) .0 ?:(=(6 b) .~0 .~~~0))))) ++ f1 |=(b=@ ^-(@ ?:(=(4 b) .~~1 ?:(=(5 b) .1 ?:(=(6 b) .~1 .~~~1))))) ++ f2 |=(b=@ ^-(@ ?:(=(4 b) .~~2 ?:(=(5 b) .2 ?:(=(6 b) .~2 .~~~2))))) @@ -949,5 +951,200 @@ :: ~[2 2] :: Source ++ eigvecs |=(a=ray:ls ^-(ray:ls +:(eig a))) + :: + +| %rand + :: + :: +rand-ray (rand-spec.md section 8): fills a Lagoon $ray with random + :: values from /lib/rand's engines. Element order is ROW-MAJOR (C + :: order) -- fixed so a future jet parallelizing across elements lands + :: identical bits regardless of thread scheduling. + :: + :: NAMING FOOTGUN: this core (+sa) already shadows the stdlib + :: add/sub/mul/div/lth/lte/gth/gte (see the +| %uno arms above, e.g. + :: +lth wraps ray comparison), so every arithmetic/compare op below + :: uses the ^-escaped stdlib form -- a bare `lth`/`add`/`mul` here + :: would silently call the RAY-valued version instead of the integer + :: one. `con`/`lsh`/`bex`/`roll` are stdlib and safe bare. `+zeros` is + :: NOT one of `*lagoon`'s star-imported top-level arms -- it lives + :: inside `+la`'s own inner door -- so it needs `zeros:la` (matching + :: how the rest of this file already reaches it, e.g. `zeros:(lake + :: rnd)` above), not a bare call. + :: + :: %phil (Philox4x32-10) gets special per-element counter treatment, + :: since it's the engine meant to be jet-parallelized (rand-spec.md + :: section 3.1): +fill-uniform (a single, non-rejecting draw per + :: element) assigns element i counter ctr0+i directly, so the whole + :: fill decomposes into n independent, order-free draws. +fill-normal/ + :: +fill-expon/+fill-below are rejection-based (variable-consumption) + :: transforms, so they get a WIDER per-element counter WINDOW instead + :: (ctr0 + i*2^32): the rejection loop walks freely within its own + :: window, and the returned rng's counter is forced to ctr0+n*2^32 + :: regardless of how many sub-draws each element actually used -- so + :: the post-state is a pure function of n, not of how lucky/unlucky + :: each element's rejection loop was (needed for replay and for + :: composing a following +fill-* call cleanly). Window exhaustion + :: (walking past 2^32 sub-draws for one element) is astronomically + :: improbable but crashes rather than silently overflowing into the + :: next element's window. + :: + :: Non-%phil engines (%sm64, %pcg) have no cheap "jump to counter N" + :: operation generic enough to reuse here, and aren't the ones a jet + :: would parallelize anyway, so they just thread the rng sequentially + :: through each element in row-major order -- no window, no per- + :: element counter surgery, ordinary sequential composition. + :: + :: +fill-simple: the +fill-uniform case (no window; ctr0+i per + :: element for %phil). Shared so the packing loop exists once, not + :: once per draw kind -- `draw` closes over anything the specific + :: distribution needs (e.g. +fill-expon's lambda) before being passed + :: in. + ++ fill-simple + |= [meta=meta:ls r=rng:rand draw=$-(rng:rand [@ rng:rand])] + ^- [ray=ray:ls r=rng:rand] + =/ n (roll shape.meta ^mul) + =/ acc +:(zeros:la meta) + ?: ?=(%phil -.r) + :: Fresh [%phil key0 ctrN] literals throughout, rather than + :: mutating .r (`r(ctr.p ...)`), so nothing inside the recursive + :: trap depends on this ?= narrowing persisting across $ calls -- + :: a known Hoon footgun class (mint-vain/nest-fail from narrowing + :: not surviving a recursive rebind). + =/ key0 key.p.r + =/ ctr0 ctr.p.r + =/ i 0 + |- ^- [ray=ray:ls r=rng:rand] + ?: =(i n) + [[meta acc] [%phil key0 (^add ctr0 n)]] + =/ v -:(draw [%phil key0 (^add ctr0 i)]) + $(i +(i), acc (con acc (lsh [bloq.meta i] v))) + :: Explicit widen: the ?: ?=(%phil -.r) above narrows .r to exclude + :: %phil for this fall-through, but +draw's OWN return type is the + :: full rng:rand union (any engine can come back out) -- so the + :: trap's first entry (narrow) and its recursive re-entries (wide, + :: after the first =^ v r (draw r)) would otherwise disagree on .r's + :: type. Widening up front once makes every entry consistent. + =/ r0 `rng:rand`r + =/ i 0 + |- ^- [ray=ray:ls r=rng:rand] + ?: =(i n) + [[meta acc] r0] + =^ v r0 (draw r0) + $(i +(i), acc (con acc (lsh [bloq.meta i] v))) + :: +fill-windowed: the +fill-normal/+fill-expon/+fill-below case -- + :: rejection-based, so %phil elements get the ctr0+i*2^32 window + :: instead of a bare ctr0+i, with an explicit crash if a single + :: element's rejection loop ever walks past its own window. + ++ fill-windowed + |= [meta=meta:ls r=rng:rand draw=$-(rng:rand [@ rng:rand])] + ^- [ray=ray:ls r=rng:rand] + =/ n (roll shape.meta ^mul) + =/ acc +:(zeros:la meta) + ?: ?=(%phil -.r) + :: Same fresh-literal approach as +fill-simple, see its comment. + =/ key0 key.p.r + =/ ctr0 ctr.p.r + =/ win (bex 32) + =/ i 0 + |- ^- [ray=ray:ls r=rng:rand] + ?: =(i n) + [[meta acc] [%phil key0 (^add ctr0 (^mul n win))]] + =/ ri `rng:rand`[%phil key0 (^add ctr0 (^mul i win))] + =^ v ri (draw ri) + ~| %rand-ray-window-exhausted + ?> ?=(%phil -.ri) + ?> (^lth ctr.p.ri (^add ctr0 (^mul +(i) win))) + $(i +(i), acc (con acc (lsh [bloq.meta i] v))) + :: Explicit widen: the ?: ?=(%phil -.r) above narrows .r to exclude + :: %phil for this fall-through, but +draw's OWN return type is the + :: full rng:rand union (any engine can come back out) -- so the + :: trap's first entry (narrow) and its recursive re-entries (wide, + :: after the first =^ v r (draw r)) would otherwise disagree on .r's + :: type. Widening up front once makes every entry consistent. + =/ r0 `rng:rand`r + =/ i 0 + |- ^- [ray=ray:ls r=rng:rand] + ?: =(i n) + [[meta acc] r0] + =^ v r0 (draw r0) + $(i +(i), acc (con acc (lsh [bloq.meta i] v))) + :: +fill-uniform: [meta:ls rng:rand] -> [ray:ls rng:rand] + :: + :: Uniform [0,1) at the meta's float bloq (5=@rs, 6=@rd only -- v1 + :: scope per rand-spec.md section 8). + :: Examples + :: > (fill-uniform:sa [~[2 2] 5 %i754 ~] (from-atom:seed:rand %sm64 0)) + :: [ [meta=[shape=~[2 2] bloq=5 kind=%i754 tail=0] + :: data=0x1.3e99.03d8.3d14.54f0.3f39.65f4.3dee.6d78] + :: r=[%sm64 s=0x78dd.e6e5.fd29.f054] ] + :: Source + ++ fill-uniform + |= [meta=meta:ls r=rng:rand] + ^- [ray=ray:ls r=rng:rand] + ~| %rand-ray-bad-kind ?> =(%i754 kind.meta) + ~| %rand-ray-bad-bloq ?> |(=(5 bloq.meta) =(6 bloq.meta)) + %^ fill-simple meta r + ?: =(5 bloq.meta) + |=(rr=rng:rand (rs:uni:i754rand rr)) + |=(rr=rng:rand (rd:uni:i754rand rr)) + :: +fill-normal: [meta:ls rng:rand] -> [ray:ls rng:rand] + :: + :: Standard normal N(0,1) at the meta's float bloq (5/6 only). + :: Rejection-based (Marsaglia polar, via /lib/i754rand's +normal), so + :: %phil elements use the counter-window scheme (see header). + :: Examples + :: > (fill-normal:sa [~[2 2] 5 %i754 ~] (from-atom:seed:rand %sm64 0)) + :: [ [meta=[shape=~[2 2] bloq=5 kind=%i754 tail=0] + :: data=0x1.3fb9.4bba.be73.0e84.4009.430e.bf17.e80f] + :: r=[%sm64 s=0x2e2a.c13e.f8e8.d8d2] ] + :: Source + ++ fill-normal + |= [meta=meta:ls r=rng:rand] + ^- [ray=ray:ls r=rng:rand] + ~| %rand-ray-bad-kind ?> =(%i754 kind.meta) + ~| %rand-ray-bad-bloq ?> |(=(5 bloq.meta) =(6 bloq.meta)) + %^ fill-windowed meta r + ?: =(5 bloq.meta) + |=(rr=rng:rand (normal:rs:dist:i754rand rr)) + |=(rr=rng:rand (normal:rd:dist:i754rand rr)) + :: +fill-expon: [meta:ls rng:rand lambda=@] -> [ray:ls rng:rand] + :: + :: Exponential(lambda) at the meta's float bloq (5/6 only); .lambda is + :: a raw @ reinterpreted at that bloq's aura. Rejection-based (inverse + :: CDF via an open-open uniform draw), so %phil elements use the + :: counter-window scheme (see header). + :: Examples + :: > (fill-expon:sa [~[2 2] 6 %i754 ~] (from-atom:seed:rand %sm64 0) .~1) + :: [ [meta=[shape=~[2 2] bloq=6 kind=%i754 tail=0] + :: data=0x1.3ff0.11b8.a292.424b.3fff.e0f9.f3ad.2c6a.3fd0.c84b.1505 + :: .9506.400f.15bc.a87f.d682] + :: r=[%sm64 s=0x78dd.e6e5.fd29.f054] ] + :: Source + ++ fill-expon + |= [meta=meta:ls r=rng:rand lambda=@] + ^- [ray=ray:ls r=rng:rand] + ~| %rand-ray-bad-kind ?> =(%i754 kind.meta) + ~| %rand-ray-bad-bloq ?> |(=(5 bloq.meta) =(6 bloq.meta)) + %^ fill-windowed meta r + ?: =(5 bloq.meta) + |=(rr=rng:rand (expon:rs:dist:i754rand rr `@rs`lambda)) + |=(rr=rng:rand (expon:rd:dist:i754rand rr `@rd`lambda)) + :: +fill-below: [meta:ls n=@ rng:rand] -> [ray:ls rng:rand] + :: + :: Uniform in [0,n) via Lemire (/lib/rand's +below), for %uint rays at + :: any bloq -- caller's responsibility that .n fits the meta's bloq + :: width. Rejection-based, so %phil elements use the counter-window + :: scheme (see header). + :: Examples + :: > (fill-below:sa [~[2 2] 5 %uint ~] 100 (from-atom:seed:rand %sm64 0)) + :: [ [meta=[shape=~[2 2] bloq=5 kind=%uint tail=0] + :: data=0x1.0000.0061.0000.0002.0000.002b.0000.0058] + :: r=[%sm64 s=0x78dd.e6e5.fd29.f054] ] + :: Source + ++ fill-below + |= [meta=meta:ls n=@ r=rng:rand] + ^- [ray=ray:ls r=rng:rand] + ~| %rand-ray-bad-kind ?> =(%uint kind.meta) + %^ fill-windowed meta r + |=(rr=rng:rand (below:uni:rand rr n)) -- -- diff --git a/saloon/desk/tests/lib/saloon-rand-ray.hoon b/saloon/desk/tests/lib/saloon-rand-ray.hoon new file mode 100644 index 0000000..c99e4e7 --- /dev/null +++ b/saloon/desk/tests/lib/saloon-rand-ray.hoon @@ -0,0 +1,141 @@ +:: Tests for Saloon's +rand-ray layer (rand-spec.md section 8): +fill- +:: uniform/+fill-normal/+fill-expon/+fill-below in the +sa core. +:: +:: Every value below is cross-checked against a DIRECT call to the +:: underlying /lib/rand or /lib/i754rand primitive at the expected +:: counter/state, not hand-derived -- proving the per-element counter +:: math (rand-spec.md section 8's %phil "ctr0+i" / "ctr0+i*2^32 window" +:: design), not just that the ray-filling loop runs without crashing. +:: +/- ls=lagoon +/+ *test, *saloon, *lagoon, rand, i754rand +|% +:: +fill-uniform, %phil: element i must equal a plain rs:uni:i754rand +:: draw at counter i (rand-spec.md: "element i uses counter ctr0+i"), +:: and the post-state counter must be exactly ctr0+n (here 0+3=3). +++ test-fill-uniform-phil ^- tang + =/ res (fill-uniform:sa [~[3] 5 %i754 ~] (from-atom:seed:rand %phil 0)) + ?> ?=(%phil -.r.res) + =/ key0 key.p.r.res + ;: weld + %+ expect-eq !>(`@`3) !>(ctr.p.r.res) + %+ expect-eq + !>(out:(rs:uni:i754rand [%phil p=[key0 0]])) + !>((end [0 32] data.ray.res)) + %+ expect-eq + !>(out:(rs:uni:i754rand [%phil p=[key0 1]])) + !>((cut 0 [32 32] data.ray.res)) + %+ expect-eq + !>(out:(rs:uni:i754rand [%phil p=[key0 2]])) + !>((cut 0 [64 32] data.ray.res)) + == +:: +fill-uniform, %sm64 (sequential): elements thread the rng state +:: exactly the way two direct, sequential rs:uni:i754rand calls would. +++ test-fill-uniform-sm64 ^- tang + =/ res (fill-uniform:sa [~[2] 5 %i754 ~] (from-atom:seed:rand %sm64 0)) + =/ a0 (from-atom:seed:rand %sm64 0) + =^ e0 a0 (rs:uni:i754rand a0) + =^ e1 a0 (rs:uni:i754rand a0) + ;: weld + %+ expect-eq !>(e0) !>((end [0 32] data.ray.res)) + %+ expect-eq !>(e1) !>((cut 0 [32 32] data.ray.res)) + %+ expect-eq !>(a0) !>(r.res) + == +:: +fill-normal, %phil: element i must equal a normal:rd:dist:i754rand +:: draw at counter i*2^32 (rand-spec.md's counter-WINDOW design for +:: rejection-based transforms), and the post-state must be exactly +:: ctr0 + n*2^32 (here 0 + 2*2^32), regardless of how many sub-draws +:: each element's rejection loop actually used. +++ test-fill-normal-phil ^- tang + =/ res (fill-normal:sa [~[2] 6 %i754 ~] (from-atom:seed:rand %phil 0)) + ?> ?=(%phil -.r.res) + =/ key0 key.p.r.res + =/ win (bex 32) + ;: weld + %+ expect-eq !>(`@`(mul 2 win)) !>(ctr.p.r.res) + %+ expect-eq + !>(out:(normal:rd:dist:i754rand [%phil p=[key0 0]])) + !>((end [0 64] data.ray.res)) + %+ expect-eq + !>(out:(normal:rd:dist:i754rand [%phil p=[key0 win]])) + !>((cut 0 [64 64] data.ray.res)) + == +:: +fill-expon, %phil: same window scheme as +fill-normal, with the +:: distribution's own lambda parameter threaded through. +++ test-fill-expon-phil ^- tang + =/ res (fill-expon:sa [~[2] 6 %i754 ~] (from-atom:seed:rand %phil 0) .~2) + ?> ?=(%phil -.r.res) + =/ key0 key.p.r.res + =/ win (bex 32) + ;: weld + %+ expect-eq !>(`@`(mul 2 win)) !>(ctr.p.r.res) + %+ expect-eq + !>(out:(expon:rd:dist:i754rand [%phil p=[key0 0]] .~2)) + !>((end [0 64] data.ray.res)) + == +:: +fill-normal, %sm64 (sequential): same cross-check style as +:: +fill-uniform's %sm64 test. +++ test-fill-normal-sm64 ^- tang + =/ res (fill-normal:sa [~[2] 6 %i754 ~] (from-atom:seed:rand %sm64 0)) + =/ a0 (from-atom:seed:rand %sm64 0) + =^ e0 a0 (normal:rd:dist:i754rand a0) + =^ e1 a0 (normal:rd:dist:i754rand a0) + ;: weld + %+ expect-eq !>(e0) !>((end [0 64] data.ray.res)) + %+ expect-eq !>(e1) !>((cut 0 [64 64] data.ray.res)) + %+ expect-eq !>(a0) !>(r.res) + == +:: +fill-expon, %sm64 (sequential). +++ test-fill-expon-sm64 ^- tang + =/ res (fill-expon:sa [~[2] 6 %i754 ~] (from-atom:seed:rand %sm64 0) .~2) + =/ a0 (from-atom:seed:rand %sm64 0) + =^ e0 a0 (expon:rd:dist:i754rand a0 .~2) + =^ e1 a0 (expon:rd:dist:i754rand a0 .~2) + ;: weld + %+ expect-eq !>(e0) !>((end [0 64] data.ray.res)) + %+ expect-eq !>(e1) !>((cut 0 [64 64] data.ray.res)) + %+ expect-eq !>(a0) !>(r.res) + == +:: +fill-below, %sm64: elements are plain sequential Lemire draws. +++ test-fill-below-sm64 ^- tang + =/ res (fill-below:sa [~[1] 5 %uint ~] 100 (from-atom:seed:rand %sm64 0)) + %+ expect-eq + !>(out:(below:uni:rand (from-atom:seed:rand %sm64 0) 100)) + !>((end [0 32] data.ray.res)) +:: Domain violations crash, tagged (rand-spec.md section 9 policy). +++ test-fill-uniform-bad-kind-crashes ^- tang + %- expect-fail + |.((fill-uniform:sa [~[2] 5 %uint ~] (from-atom:seed:rand %sm64 0))) +++ test-fill-uniform-bad-bloq-crashes ^- tang + %- expect-fail + |.((fill-uniform:sa [~[2] 7 %i754 ~] (from-atom:seed:rand %sm64 0))) +++ test-fill-normal-bad-kind-crashes ^- tang + %- expect-fail + |.((fill-normal:sa [~[2] 5 %uint ~] (from-atom:seed:rand %sm64 0))) +++ test-fill-normal-bad-bloq-crashes ^- tang + %- expect-fail + |.((fill-normal:sa [~[2] 7 %i754 ~] (from-atom:seed:rand %sm64 0))) +++ test-fill-expon-bad-kind-crashes ^- tang + %- expect-fail + |.((fill-expon:sa [~[2] 5 %uint ~] (from-atom:seed:rand %sm64 0) .~2)) +++ test-fill-expon-bad-bloq-crashes ^- tang + %- expect-fail + |.((fill-expon:sa [~[2] 7 %i754 ~] (from-atom:seed:rand %sm64 0) .~2)) +++ test-fill-below-bad-kind-crashes ^- tang + %- expect-fail + |.((fill-below:sa [~[2] 5 %i754 ~] 100 (from-atom:seed:rand %sm64 0))) +:: +fill-windowed's own window-exhaustion crash (rand-spec.md section 8: +:: "window exhaustion... crash if it happens"), exercised directly since +:: no real /lib/i754rand rejection loop will ever plausibly walk past a +:: 2^32-draw window -- a hand-built .draw gate that always claims to +:: have consumed 2^40 sub-draws stands in for that astronomically +:: improbable case. +++ test-fill-window-exhausted-crashes ^- tang + %- expect-fail + |. + %^ fill-windowed:sa [~[1] 6 %i754 ~] (from-atom:seed:rand %phil 0) + |= rr=rng:rand + ^- [@ rng:rand] + ?> ?=(%phil -.rr) + [0 rr(ctr.p (add ctr.p.rr (bex 40)))] +--