diff --git a/.github/workflows/oracle.yml b/.github/workflows/oracle.yml
new file mode 100644
index 0000000..a08533e
--- /dev/null
+++ b/.github/workflows/oracle.yml
@@ -0,0 +1,28 @@
+name: Regenerate numerical oracle
+
+on:
+ schedule:
+ # The full supply-scale Arb corpus takes roughly two minutes. Running it
+ # nightly catches certificate/generator drift without burdening every PR.
+ - cron: "17 5 * * *"
+ workflow_dispatch:
+
+permissions:
+ contents: read
+
+jobs:
+ arb-certificates:
+ name: Verify Arb certificates
+ runs-on: ubuntu-latest
+ timeout-minutes: 15
+ steps:
+ - uses: actions/checkout@v7
+ - uses: actions/setup-python@v6
+ with:
+ python-version: "3.11"
+ cache: pip
+ cache-dependency-path: tools/requirements-oracle.txt
+ - name: Install python-flint
+ run: python -m pip install -r tools/requirements-oracle.txt
+ - name: Regenerate and compare certified quantiles
+ run: python tools/generate_arb_oracle.py --check testdata/f128_arb_certificates.json
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index f278759..7d48bca 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -6,17 +6,120 @@ on:
- main
pull_request:
+concurrency:
+ group: test-${{ github.event.pull_request.number || github.ref }}
+ cancel-in-progress: true
+
jobs:
build:
+ name: Tests (${{ matrix.os }})
+ runs-on: ${{ matrix.os }}
+ strategy:
+ matrix:
+ # Both architectures on each OS. Exact-output tests pin both paths, so
+ # a floating-point divergence in the cgo/Boost CDF between these
+ # targets shows up as a failure rather than a silent consensus split.
+ os:
+ - ubuntu-latest # linux/amd64
+ - ubuntu-24.04-arm # linux/arm64
+ - macos-latest # darwin/arm64
+ - macos-15-intel # darwin/amd64
+ - windows-latest # windows/amd64
+ steps:
+ - uses: actions/checkout@v7
+
+ - name: Set up Go
+ uses: actions/setup-go@v7
+ with:
+ go-version-file: go.mod
+
+ - name: Build
+ run: go build -v ./...
+
+ - name: Test
+ run: go test -v ./...
+
+ # Ordinary `go test` runs the rapid property tests at their default 100
+ # cases; run them again at a depth that reliably samples every generator
+ # band (the divU shallow-quotient bug reproduced within ~30 banded cases,
+ # so 20000 leaves ample margin while staying under a few seconds).
+ - name: Rapid property tests
+ shell: bash
+ run: go test -v -run TestRapid ./... -args -rapid.checks=20000
+
+ # Cross-checks the "bit-identical on every platform" claim on a 32-bit-int
+ # target, where Go's int is 32 bits -- the reason f128.exp is explicitly
+ # int64. x86-64 runs i386 user code natively, so this needs no emulation;
+ # only a 32-bit C/C++ toolchain for the cgo/Boost build, which the Go
+ # toolchain invokes with -m32 automatically under GOARCH=386. armv7 (the
+ # Raspberry Pi target) has the same 32-bit int but would need slow QEMU, so
+ # 386 buys the same integer width for free.
+ test-386:
+ name: Tests (linux/386, native 32-bit)
runs-on: ubuntu-latest
+ env:
+ GOARCH: "386"
+ CGO_ENABLED: "1"
steps:
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v7
- name: Set up Go
- uses: actions/setup-go@v3
+ uses: actions/setup-go@v7
+ with:
+ go-version-file: go.mod
+
+ - name: Install 32-bit C/C++ toolchain
+ run: sudo apt-get update && sudo apt-get install -y gcc-multilib g++-multilib
- name: Build
run: go build -v ./...
- name: Test
run: go test -v ./...
+
+ - name: Rapid property tests
+ shell: bash
+ run: go test -v -run TestRapid ./... -args -rapid.checks=20000
+
+ # Same 32-bit-int check as the 386 leg, but on the real Raspberry Pi target
+ # (armv7, GOARM=7). ARM has no i386-style native path on the amd64 runner, so
+ # it runs under QEMU in a foreign-arch container; the stock golang image
+ # carries gcc/g++ and the repo vendors Boost, so nothing extra is installed.
+ # Emulated, hence the slowest leg -- it runs the default rapid depth only,
+ # leaving the 20000-case sweep to the native legs.
+ test-armv7:
+ name: Tests (linux/arm/v7, emulated)
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v7
+
+ - name: Set up QEMU
+ uses: docker/setup-qemu-action@v3
+
+ - name: Test under QEMU (armv7 -- Raspberry Pi, 32-bit int)
+ run: |
+ docker run --rm --platform linux/arm/v7 \
+ -v "$PWD":/src -w /src -e CGO_ENABLED=1 \
+ golang:1.23 \
+ go test -v ./...
+
+ # Applies each curated single-site mutant of the f128 arithmetic and
+ # verifies the fast suite kills it (survivors marked equivalent carry
+ # written proofs and are expected). This keeps the kill-rate claim true
+ # over time: a weakened test lets a mutant survive, and a semantic change
+ # to the arithmetic makes an equivalent mutant die or a target string
+ # disappear -- both fail the job, deliberately coupling the mutant table to
+ # the consensus-critical lines it guards. Platform-independent, so one leg.
+ mutation:
+ name: Mutation campaign
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v7
+
+ - name: Set up Go
+ uses: actions/setup-go@v7
+ with:
+ go-version-file: go.mod
+
+ - name: Run mutation campaign
+ run: go run mutation_check.go
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..1cf67b5
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+.venv-oracle/
diff --git a/README.md b/README.md
index ada67b5..01a5b44 100644
--- a/README.md
+++ b/README.md
@@ -5,8 +5,11 @@ This package provides an implementation of cryptographic sortition used by [go-a
Please visit the [go-algorand README](https://github.com/algorand/go-algorand/blob/master/README.md) for more information on building and using this software.
+See [TESTING.md](TESTING.md) for the testing strategy, oracle layers, and
+developer runbook for the deterministic sortition implementation.
+
## License
Please see the [COPYING_FAQ](COPYING_FAQ) for details about how to apply our license.
-Copyright (C) 2019-2023, Algorand Inc.
+Copyright (C) 2019-2026, Algorand Foundation Ltd.
diff --git a/TESTING.md b/TESTING.md
new file mode 100644
index 0000000..89351a3
--- /dev/null
+++ b/TESTING.md
@@ -0,0 +1,95 @@
+# Testing sortition
+
+The test suite combines bit-identical differential testing, independent exact
+mathematics, metamorphic properties, certified offline vectors, and mutation
+testing. No single oracle is treated as sufficient for the consensus-critical
+`SelectF128` path.
+
+## Common commands
+
+```sh
+# Complete ordinary suite
+go test -count=1 ./...
+
+# Property tests at the depth used by CI
+go test -run TestRapid -rapid.checks=20000 -count=1 ./...
+
+# Curated production and test-oracle mutation campaign
+go run mutation_check.go
+
+# Longer fuzzing runs
+go test -run xxx -fuzz FuzzF128Ops -fuzztime 10m
+go test -run xxx -fuzz FuzzSelectF128 -fuzztime 10m
+```
+
+## Assurance layers
+
+| Approach | Purpose | Location |
+|---|---|---|
+| Ordinary and compatibility tests | Basic sortition behavior, benchmarks, agreement with the deployed C++/Boost path, defined frozen-tail behavior, and current-scale regression cases | `sortition_test.go`, `f128_test.go` |
+| Bit-identical 128-bit oracle | Mirrors the production CDF walk with 128-bit `math/big.Float`; detects implementation defects in the hand-written integer arithmetic | `selectBigOracle` and related tests in `f128_test.go` |
+| Primitive differential fuzzing | Checks f128 arithmetic against `big.Float`, with constructed tie/carry seeds and magnitude-banded rapid generators | `FuzzF128Ops` in `f128_test.go`, `TestRapidF128Ops` in `f128_rapid_test.go` |
+| Exact-integer RNE oracle | Computes normalization and round-to-nearest-even solely with `big.Int` quotient/remainder arithmetic, avoiding shared reliance on `big.Float` | `f128_integer_oracle_test.go` |
+| Division certificates | Proves sampled `divStep` digits satisfy `U = q*V + R` and `0 <= R < V`, including quotient-estimate corner paths | `TestDivStepExactCertificate` in `f128_integer_oracle_test.go` |
+| Exact binomial formula | Builds the CDF from binomial coefficients and integer cross-products rather than the production PMF recurrence | `f128_exact_test.go`, `cdf_exhaustive_test.go` |
+| Exhaustive boundary grids | Enumerates small parameter domains and checks digest integers below, at, and above every exact CDF boundary, with anti-vacuity counters | `cdf_exhaustive_test.go`, `oracle_hardening_test.go` |
+| Higher-precision convergence | Requires the 512- and 1024-bit CDF walks to converge and compares them with exact small-domain results and certified large-money results | `selectAtPrecision` in `f128_exact_test.go`, tests in `oracle_hardening_test.go` |
+| Internal trajectory and liveness | Checks PMF/CDF error envelopes and monotonicity, audits the maximum-domain exponent path, proves bounded freeze permanence, and pins CDF-evaluation limits | `f128_trajectory_test.go` |
+| Metamorphic properties | Checks digest monotonicity, power-of-two and arbitrary common-factor probability scaling, primitive identities, and arithmetic order properties without a numeric oracle | `f128_rapid_test.go` |
+| Distribution sanity | Checks aggregate selection weight against the expected binomial mean without reusing the CDF formula | `TestSelectF128Distribution` in `f128_exact_test.go` |
+| Arb-certified quantiles | Uses Arb's regularized incomplete beta implementation to certify large-money quantile inequalities with rigorous dyadic endpoints | `tools/generate_arb_oracle.py`, `testdata/f128_arb_certificates.json`, `f128_arb_certificate_test.go` |
+| Mutation testing | Applies curated bugs to production code and in-process test oracles; every non-equivalent mutant must be killed | `mutation_check.go` |
+
+## Arb corpus
+
+The Arb generator is an offline tool and is not a dependency of ordinary
+`go test`. Set it up in an ignored virtual environment:
+
+```sh
+python3 -m venv .venv-oracle
+.venv-oracle/bin/pip install -r tools/requirements-oracle.txt
+```
+
+Regenerate the corpus or verify it without rewriting the checked-in file:
+
+```sh
+.venv-oracle/bin/python tools/generate_arb_oracle.py \
+ > testdata/f128_arb_certificates.json
+
+.venv-oracle/bin/python tools/generate_arb_oracle.py \
+ --check testdata/f128_arb_certificates.json
+```
+
+The Go certificate test exactly compares the stored dyadic endpoints with the
+digest ratio. Trust in the external tool is limited to Arb's assertion that
+those endpoints enclose the true incomplete-beta CDF.
+
+## Important test semantics
+
+- The near-one frozen-tail sliver is defined to return `money`. Exact-math
+ tolerance tests exclude that interval deliberately; dedicated tests pin its
+ behavior and liveness at current-scale values.
+- `money` must remain below `SelectF128MaxMoney`. The maximum-domain exponent
+ audit checks the worst representable `1-p` combination without performing a
+ supply-sized walk.
+- Large-money in-process oracles have a committee-scale step budget. A broken
+ reference must fail loudly instead of hanging while attempting trillions of
+ iterations.
+- Randomized tests use deterministic seeds. Measure-zero events such as exact
+ ties and inclusive boundaries are constructed explicitly rather than left
+ to sampling.
+- Tests with skip or clamp paths need anti-vacuity counts or mandatory spot
+ cases proving that meaningful assertions ran.
+
+## Continuous integration
+
+- `.github/workflows/test.yml` runs builds, ordinary tests, and 20,000-check
+ rapid properties across Linux amd64/arm64, macOS amd64/arm64, and Windows
+ amd64. It also runs the mutation campaign on Linux.
+- `.github/workflows/oracle.yml` regenerates and compares the complete Arb
+ corpus nightly and on manual dispatch.
+
+When arithmetic or recurrence semantics change, update the relevant
+independent oracle and exact/certified vectors deliberately, then rerun the
+full suite and mutation campaign. A changed knife-edge pin is evidence that
+the numerical trajectory changed and should be explained in the commit.
diff --git a/cdf_exhaustive_test.go b/cdf_exhaustive_test.go
new file mode 100644
index 0000000..fcac795
--- /dev/null
+++ b/cdf_exhaustive_test.go
@@ -0,0 +1,176 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "math/big"
+ "testing"
+)
+
+func exactSelectFromCDFNumerators(money uint64, cdfNums []*big.Int, cdfDen, digestDen, digest *big.Int) uint64 {
+ lhs := new(big.Int).Mul(digest, cdfDen)
+ for j, cdfNum := range cdfNums {
+ rhs := new(big.Int).Mul(cdfNum, digestDen)
+ if lhs.Cmp(rhs) <= 0 {
+ return uint64(j)
+ }
+ }
+ return money
+}
+
+func digestFromBigInt(t *testing.T, x, digestDen *big.Int) Digest {
+ t.Helper()
+ if x.Sign() < 0 || x.Cmp(digestDen) > 0 {
+ t.Fatalf("digest integer outside [0, 2^256-1]: %s", x)
+ }
+ var d Digest
+ x.FillBytes(d[:])
+ return d
+}
+
+func uint64Distance(a, b uint64) uint64 {
+ if a >= b {
+ return a - b
+ }
+ return b - a
+}
+
+// TestSelectF128ExhaustiveSmallDomain complements the sampled exact-rational
+// test with a complete finite grid. Every non-degenerate rational probability,
+// every CDF boundary, and the adjacent digest integers are constructed from
+// the first-principles binomial formula. The f128 walk is allowed to move one
+// boundary because its per-operation RNE semantics intentionally differ from
+// exact mathematics at knife edges.
+func TestSelectF128ExhaustiveSmallDomain(t *testing.T) {
+ digestDen := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ const maxMoney = uint64(16)
+ const maxTotal = uint64(32)
+
+ var parameters, boundaries, bracketed, candidates, asserted, skippedSliver uint64
+ for money := uint64(1); money <= maxMoney; money++ {
+ for total := uint64(2); total <= maxTotal; total++ {
+ cdfDen := new(big.Int).Exp(new(big.Int).SetUint64(total), new(big.Int).SetUint64(money), nil)
+ for expected := uint64(1); expected < total; expected++ {
+ parameters++
+ cdfNums := make([]*big.Int, money)
+ points := make(map[string]*big.Int, 3*money+2)
+ points["0"] = new(big.Int)
+ points[digestDen.Text(16)] = new(big.Int).Set(digestDen)
+
+ for j := uint64(0); j < money; j++ {
+ boundaries++
+ cdfNum := exactCDFNum(money, total, expected, j)
+ cdfNums[j] = cdfNum
+ scaled := new(big.Int).Mul(cdfNum, digestDen)
+ at, rem := new(big.Int), new(big.Int)
+ at.QuoRem(scaled, cdfDen, rem)
+
+ // at/digestDen <= CDF(j) < (at+1)/digestDen. This
+ // assertion makes every generated window prove that it
+ // really straddles the intended exact boundary.
+ atScaled := new(big.Int).Mul(new(big.Int).Set(at), cdfDen)
+ above := new(big.Int).Add(new(big.Int).Set(at), big.NewInt(1))
+ aboveScaled := new(big.Int).Mul(new(big.Int).Set(above), cdfDen)
+ if atScaled.Cmp(scaled) > 0 || aboveScaled.Cmp(scaled) <= 0 || at.Cmp(digestDen) >= 0 {
+ t.Fatalf("invalid exact boundary bracket (money=%d total=%d expected=%d j=%d at=%s remainder=%s)",
+ money, total, expected, j, at, rem)
+ }
+ bracketed++
+
+ if at.Sign() > 0 {
+ below := new(big.Int).Sub(new(big.Int).Set(at), big.NewInt(1))
+ points[below.Text(16)] = below
+ }
+ points[at.Text(16)] = new(big.Int).Set(at)
+ points[above.Text(16)] = above
+ }
+
+ for _, point := range points {
+ candidates++
+ d := digestFromBigInt(t, point, digestDen)
+ got := SelectF128(money, total, expected, d)
+ if bitExact := selectBigOracle(money, total, expected, d); got != bitExact {
+ t.Fatalf("SelectF128=%d vs 128-bit walk oracle=%d (money=%d total=%d expected=%d digest=%x)",
+ got, bitExact, money, total, expected, d)
+ }
+ if inFrozenSliver(money, digestRatioBig(d, 512)) {
+ skippedSliver++
+ continue
+ }
+ want := exactSelectFromCDFNumerators(money, cdfNums, cdfDen, digestDen, point)
+ asserted++
+ if uint64Distance(got, want) > 1 {
+ t.Fatalf("SelectF128=%d vs exact=%d (money=%d total=%d expected=%d digest=%x)",
+ got, want, money, total, expected, d)
+ }
+ }
+ }
+ }
+ }
+
+ if boundaries == 0 || bracketed != boundaries {
+ t.Fatalf("anti-vacuity: bracketed %d of %d exact boundaries", bracketed, boundaries)
+ }
+ if parameters < 7_000 || asserted < 100_000 || skippedSliver >= candidates/2 {
+ t.Fatalf("anti-vacuity: parameters=%d boundaries=%d candidates=%d asserted=%d skippedSliver=%d",
+ parameters, boundaries, candidates, asserted, skippedSliver)
+ }
+ t.Logf("exhaustive exact grid: parameters=%d boundaries=%d candidates=%d asserted=%d skippedSliver=%d",
+ parameters, boundaries, candidates, asserted, skippedSliver)
+}
+
+// TestSelectF128ExhaustiveSmallDegenerateDomain keeps p=0, p>=1,
+// totalMoney=0, and money=0 out of the regular-grid special cases above. It
+// exhausts them separately with endpoint digests so none is lost to a skip.
+func TestSelectF128ExhaustiveSmallDegenerateDomain(t *testing.T) {
+ var maximum Digest
+ for i := range maximum {
+ maximum[i] = 0xff
+ }
+ var one Digest
+ one[len(one)-1] = 1
+ digests := []Digest{{}, one, maximum}
+
+ for money := uint64(0); money <= 16; money++ {
+ for total := uint64(0); total <= 32; total++ {
+ expectedValues := map[uint64]struct{}{0: {}, total: {}}
+ if total != ^uint64(0) {
+ expectedValues[total+1] = struct{}{}
+ }
+ for expected := range expectedValues {
+ for _, d := range digests {
+ got := SelectF128(money, total, expected, d)
+ var want uint64
+ switch {
+ case money == 0:
+ want = 0
+ case expected < total: // only expected==0: all mass at zero
+ want = 0
+ case d == (Digest{}): // p>=1 and ratio==0 crosses cdf(0)==0
+ want = 0
+ default:
+ want = money
+ }
+ if got != want {
+ t.Fatalf("degenerate SelectF128=%d, want %d (money=%d total=%d expected=%d digest=%x)",
+ got, want, money, total, expected, d)
+ }
+ }
+ }
+ }
+ }
+}
diff --git a/f128.go b/f128.go
new file mode 100644
index 0000000..8d87a1e
--- /dev/null
+++ b/f128.go
@@ -0,0 +1,568 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "encoding/binary"
+ "math/bits"
+)
+
+// f128 is a minimal NON-NEGATIVE binary floating-point value with a 128-bit
+// mantissa, used for a fast, allocation-free, deterministic binomial CDF (see
+// SelectF128). It is a value type -- arithmetic returns new f128s on the stack,
+// never the heap -- and uses only integer ops (math/bits), so it is bit-identical
+// on every platform (no hardware FP, no FMA, no libm).
+//
+// value = (hi<<64 | lo) * 2^exp
+//
+// The 128-bit mantissa is normalized so bit 127 (the MSB of hi) is set, or the
+// value is zero (hi==lo==0). All sortition quantities (p, 1-p, ratio, pmf, cdf,
+// factors) are >= 0, so there is no sign bit. Arithmetic ROUNDS TO NEAREST, TIES
+// TO EVEN (matching math/big.Float), so every result is the correctly-rounded
+// 128-bit value -- this is what makes SelectF128 match the big.Float oracle.
+// It also shapes the ratio == 1.0 edge (all-0xff digest, or any digest with
+// >= 129 leading one bits): for some distributions the accumulated cdf rounds
+// up to exactly 1.0 at an early j and the walk stops there; in others it stays
+// below 1.0 for all j < money and the walk runs to money (see
+// TestSelectF128RatioExactlyOne and the SelectF128 doc comment).
+type f128 struct {
+ hi, lo uint64
+ // exp is explicitly 64-bit: int is 32 bits on 386/arm, and a
+ // platform-sized exponent would make cross-GOARCH bit-identity depend on
+ // exponents staying small rather than on the type.
+ exp int64
+}
+
+// f128MantBits is the f128 mantissa width. The big.Float oracle in the test
+// forms its constants at this same precision so they round to f128 exactly.
+const f128MantBits = 128
+
+func (a f128) isZero() bool { return a.hi == 0 && a.lo == 0 }
+
+func shl128(hi, lo uint64, n uint) (uint64, uint64) {
+ switch {
+ case n == 0:
+ return hi, lo
+ case n < 64:
+ return hi<>(64-n), lo << n
+ case n < 128:
+ return lo << (n - 64), 0
+ default:
+ return 0, 0
+ }
+}
+
+func shr128(hi, lo uint64, n uint) (uint64, uint64) {
+ switch {
+ case n == 0:
+ return hi, lo
+ case n < 64:
+ return hi >> n, lo>>n | hi<<(64-n)
+ case n < 128:
+ return 0, hi >> (n - 64)
+ default:
+ return 0, 0
+ }
+}
+
+// shr128gs shifts hi:lo right by n (1..128), returning the result plus the round
+// bit (the most-significant shifted-out bit, at position n-1) and sticky (any
+// lower shifted-out bit). Used for round-to-nearest on exponent alignment.
+func shr128gs(hi, lo uint64, n uint) (rhi, rlo uint64, round, sticky bool) {
+ rhi, rlo = shr128(hi, lo, n)
+ switch {
+ case n <= 64:
+ round = (lo>>(n-1))&1 != 0
+ if n >= 2 {
+ sticky = lo&((uint64(1)<<(n-1))-1) != 0
+ }
+ case n < 128:
+ m := n - 64
+ round = (hi>>(m-1))&1 != 0
+ if m >= 2 {
+ sticky = hi&((uint64(1)<<(m-1))-1) != 0
+ }
+ sticky = sticky || lo != 0
+ default: // n == 128
+ round = hi&(1<<63) != 0
+ sticky = (hi&^(uint64(1)<<63) != 0) || lo != 0
+ }
+ return rhi, rlo, round, sticky
+}
+
+// shl256 shifts a 256-bit value (a3:a2:a1:a0) left by n < 256.
+func shl256(a3, a2, a1, a0 uint64, n uint) (uint64, uint64, uint64, uint64) {
+ words := n / 64
+ shift := n % 64
+ in := [4]uint64{a3, a2, a1, a0}
+ var out [4]uint64
+ for i := 0; i < 4; i++ {
+ src := i + int(words)
+ if src >= len(in) {
+ break
+ }
+ out[i] = in[src] << shift
+ if shift != 0 && src+1 < len(in) {
+ out[i] |= in[src+1] >> (64 - shift)
+ }
+ }
+ return out[0], out[1], out[2], out[3]
+}
+
+// norm128 normalizes a 128-bit mantissa (shifts MSB to bit 127). No bits are
+// dropped, so no rounding is needed; used for exact conversions.
+func norm128(hi, lo uint64, exp int64) f128 {
+ if hi == 0 && lo == 0 {
+ return f128{}
+ }
+ var s int
+ if hi != 0 {
+ s = bits.LeadingZeros64(hi)
+ } else {
+ s = 64 + bits.LeadingZeros64(lo)
+ }
+ if s != 0 {
+ hi, lo = shl128(hi, lo, uint(s))
+ exp -= int64(s)
+ }
+ return f128{hi, lo, exp}
+}
+
+// roundNE rounds the normalized 128-bit mantissa hi:lo to nearest, ties to even,
+// given the round bit and sticky of the discarded tail, and renormalizes on
+// carry-out. hi:lo must already be normalized (bit 127 set).
+func roundNE(hi, lo uint64, roundBit, sticky bool, exp int64) f128 {
+ if roundBit && (sticky || lo&1 != 0) {
+ var c uint64
+ lo, c = bits.Add64(lo, 1, 0)
+ hi, c = bits.Add64(hi, 0, c)
+ if c != 0 { // mantissa overflowed to 2^128 -> renormalize to 2^127
+ return f128{1 << 63, 0, exp + 1}
+ }
+ }
+ return f128{hi, lo, exp}
+}
+
+func f128FromUint64(u uint64) f128 {
+ if u == 0 {
+ return f128{}
+ }
+ s := bits.LeadingZeros64(u)
+ return f128{u << uint(s), 0, -int64(s) - 64}
+}
+
+// f128FromDigestRatio returns digest/(2^256-1), rounded to nearest-even at
+// f128 precision, without reducing the digest to float64 first.
+func f128FromDigestRatio(d Digest) f128 {
+ w3 := binary.BigEndian.Uint64(d[0:8])
+ w2 := binary.BigEndian.Uint64(d[8:16])
+ w1 := binary.BigEndian.Uint64(d[16:24])
+ w0 := binary.BigEndian.Uint64(d[24:32])
+
+ var leading int
+ switch {
+ case w3 != 0:
+ leading = bits.LeadingZeros64(w3)
+ case w2 != 0:
+ leading = 64 + bits.LeadingZeros64(w2)
+ case w1 != 0:
+ leading = 128 + bits.LeadingZeros64(w1)
+ case w0 != 0:
+ leading = 192 + bits.LeadingZeros64(w0)
+ default:
+ return f128{}
+ }
+
+ n3, n2, n1, n0 := shl256(w3, w2, w1, w0, uint(leading))
+ roundBit := n1&(uint64(1)<<63) != 0
+ sticky := n1&^(uint64(1)<<63) != 0 || n0 != 0
+
+ // Dividing by 2^256 places the binary point directly after the digest.
+ // The real denominator is one smaller, so the exact ratio is slightly
+ // larger. That correction is at most one discarded-bit unit and changes
+ // rounding only when the 2^256 quotient is exactly halfway.
+ if roundBit && !sticky {
+ sticky = true
+ }
+ return roundNE(n3, n2, roundBit, sticky, -128-int64(leading))
+}
+
+// mul returns a*b rounded to nearest even.
+func (a f128) mul(b f128) f128 {
+ if a.isZero() || b.isZero() {
+ return f128{}
+ }
+ hhHi, hhLo := bits.Mul64(a.hi, b.hi)
+ hlHi, hlLo := bits.Mul64(a.hi, b.lo)
+ lhHi, lhLo := bits.Mul64(a.lo, b.hi)
+ llHi, llLo := bits.Mul64(a.lo, b.lo)
+ p0 := llLo
+ p1, cA := bits.Add64(llHi, hlLo, 0)
+ p1, cB := bits.Add64(p1, lhLo, 0)
+ cp1 := cA + cB
+ p2, cC := bits.Add64(hhLo, hlHi, 0)
+ p2, cD := bits.Add64(p2, lhHi, 0)
+ p2, cE := bits.Add64(p2, cp1, 0)
+ p3 := hhHi + cC + cD + cE
+ if p3&(1<<63) != 0 { // product >= 2^255: mantissa p3:p2, tail p1:p0
+ roundBit := p1&(1<<63) != 0
+ sticky := (p1&^(uint64(1)<<63) != 0) || p0 != 0
+ return roundNE(p3, p2, roundBit, sticky, a.exp+b.exp+128)
+ }
+ // product in [2^254, 2^255): shift left 1 to normalize
+ hi := p3<<1 | p2>>63
+ lo := p2<<1 | p1>>63
+ roundBit := p1&(1<<62) != 0
+ sticky := (p1&^(uint64(3)<<62) != 0) || p0 != 0
+ return roundNE(hi, lo, roundBit, sticky, a.exp+b.exp+127)
+}
+
+// divU returns a/u for an unsigned integer u (the walk's per-step denominator
+// j), rounded to nearest even.
+func (a f128) divU(u uint64) f128 {
+ if a.isZero() || u == 0 {
+ return f128{}
+ }
+ // 256-bit quotient (M/u)*2^128 = q3:q2:q1:q0 by long division. Two full
+ // fraction digits below the mantissa keep the round bit inside the
+ // computed digits even for the shallowest quotient (u > a.hi, where the
+ // quotient has only 128 significant bits), so the remainder only ever
+ // contributes to sticky, strictly below the round bit. (A 192-bit
+ // quotient with the remainder OR'd into its last digit is NOT enough:
+ // for u > ~2^62 that marker lands in the mantissa or at the round bit
+ // and breaks round-to-nearest-even about half the time.)
+ q3, r := bits.Div64(0, a.hi, u)
+ q2, r := bits.Div64(r, a.lo, u)
+ q1, r := bits.Div64(r, 0, u)
+ q0, rem := bits.Div64(r, 0, u)
+ // M >= 2^127 and u < 2^64 make the quotient >= 2^191, so at most 64
+ // leading zeros: the mantissa always comes from o3:o2 below.
+ var lz int
+ if q3 != 0 {
+ lz = bits.LeadingZeros64(q3)
+ } else {
+ lz = 64
+ }
+ o3, o2, o1, o0 := shl256(q3, q2, q1, q0, uint(lz))
+ round := o1&(uint64(1)<<63) != 0
+ sticky := o1&^(uint64(1)<<63) != 0 || o0 != 0 || rem != 0
+ return roundNE(o3, o2, round, sticky, a.exp-int64(lz))
+}
+
+// divStep is one digit step of Knuth's Algorithm D for a normalized (top bit
+// set) 128-bit divisor v1:v0: it divides the 192-bit value uHi:uMid:uLo by
+// v1:v0, where the running-remainder prefix uHi:uMid is already < v1:v0, and
+// returns the 64-bit quotient digit q and the new 128-bit remainder rHi:rLo.
+func divStep(uHi, uMid, uLo, v1, v0 uint64) (q, rHi, rLo uint64) {
+ // qhat = min((uHi:uMid)/v1, 2^64-1); Div64 requires uHi < v1, so cap when
+ // uHi == v1 (the running remainder guarantees uHi never exceeds v1).
+ var qhat, rhat uint64
+ refine := true
+ if uHi >= v1 {
+ qhat = ^uint64(0)
+ var c uint64
+ rhat, c = bits.Add64(uMid, v1, 0) // rhat = (v1:uMid) - qhat*v1 = uMid + v1
+ if c != 0 {
+ refine = false // rhat >= 2^64: the refine test is already false
+ }
+ } else {
+ qhat, rhat = bits.Div64(uHi, uMid, v1)
+ }
+ // Lower qhat (over-estimated by at most 2) until qhat*v0 <= rhat:uLo.
+ for refine {
+ hi, lo := bits.Mul64(qhat, v0)
+ if hi > rhat || (hi == rhat && lo > uLo) {
+ qhat--
+ var c uint64
+ rhat, c = bits.Add64(rhat, v1, 0)
+ if c != 0 {
+ break
+ }
+ continue
+ }
+ break
+ }
+ // u - qhat*(v1:v0), a 192-bit subtraction.
+ p1hi, p1lo := bits.Mul64(qhat, v1)
+ p0hi, p0lo := bits.Mul64(qhat, v0)
+ prodMid, c := bits.Add64(p1lo, p0hi, 0)
+ prodHi := p1hi + c
+ sLo, br := bits.Sub64(uLo, p0lo, 0)
+ sMid, br := bits.Sub64(uMid, prodMid, br)
+ _, br = bits.Sub64(uHi, prodHi, br)
+ q = qhat
+ // Defense in depth: unlike Knuth's Algorithm D, which bounds the D3
+ // adjustment at two rounds and relies on this D6 add-back, the refine loop
+ // above runs to fixpoint with an exact 128-bit test (qhat*v0 <= rhat:uLo
+ // is equivalent to U - qhat*V >= 0 given rhat's bookkeeping), so qhat
+ // should already be exact and br always 0. Kept in case that analysis is
+ // wrong; an instrumented 3M-case search never fired it.
+ if br != 0 { // qhat was 1 too large: add the divisor back
+ q--
+ sLo, c = bits.Add64(sLo, v0, 0)
+ sMid, _ = bits.Add64(sMid, v1, c)
+ }
+ return q, sMid, sLo
+}
+
+// div returns a/b rounded to nearest even, for NON-NEGATIVE operands (b != 0).
+// It divides the 256-bit a.hi:a.lo:0:0 by the normalized 128-bit b.hi:b.lo via
+// Knuth long division. Since both mantissas are in [2^127,2^128), the ratio is
+// in (0.5,2), so the 129-bit integer quotient Q is either already normalized
+// (Q < 2^128) or one bit wide (Q >= 2^128); the division remainder supplies the
+// bits below Q for round-to-nearest-even. Unlike divU (small integer divisor),
+// this handles a full f128 divisor; it is used once per setup, not in the walk.
+func (a f128) div(b f128) f128 {
+ if a.isZero() || b.isZero() {
+ return f128{}
+ }
+ v1, v0 := b.hi, b.lo
+ // Long-divide [a.hi, a.lo, 0, 0] by v1:v0, most-significant limb first.
+ var remHi, remLo uint64
+ var q1, q0 uint64 // the two low quotient digits; the top two are 0 and {0,1}
+ var q2 uint64
+ digit, remHi, remLo := divStep(remHi, remLo, a.hi, v1, v0) // = 0
+ _ = digit
+ q2, remHi, remLo = divStep(remHi, remLo, a.lo, v1, v0) // in {0,1}
+ q1, remHi, remLo = divStep(remHi, remLo, 0, v1, v0)
+ q0, remHi, remLo = divStep(remHi, remLo, 0, v1, v0)
+
+ expq := a.exp - b.exp - 128
+ if q2 != 0 { // Q in [2^128,2^129): mantissa = Q>>1, dropped low bit is round
+ mantHi := q2<<63 | q1>>1
+ mantLo := q1<<63 | q0>>1
+ round := q0&1 != 0
+ sticky := remHi != 0 || remLo != 0
+ return roundNE(mantHi, mantLo, round, sticky, expq+1)
+ }
+ // Q in [2^127,2^128): already normalized; round/sticky come from rem/b, i.e.
+ // round = (2*rem >= b), sticky = the leftover after that comparison.
+ dblLo := remLo << 1
+ dblHi := remHi<<1 | remLo>>63
+ carry := remHi >> 63
+ var round, sticky bool
+ if carry != 0 || dblHi > v1 || (dblHi == v1 && dblLo >= v0) {
+ round = true
+ sLo, br := bits.Sub64(dblLo, v0, 0)
+ sHi, _ := bits.Sub64(dblHi, v1, br)
+ sticky = sHi != 0 || sLo != 0
+ } else {
+ sticky = remHi != 0 || remLo != 0
+ }
+ return roundNE(q1, q0, round, sticky, expq)
+}
+
+// add returns a+b (both non-negative), rounded to nearest even.
+func (a f128) add(b f128) f128 {
+ if a.isZero() {
+ return b
+ }
+ if b.isZero() {
+ return a
+ }
+ if a.exp < b.exp {
+ a, b = b, a
+ }
+ // Compare in int64 before converting to a shift count: uint is 32-bit on
+ // 386/arm, where a large exponent difference would otherwise truncate.
+ if a.exp-b.exp > 128 {
+ return a // b is below the round bit
+ }
+ diff := uint(a.exp - b.exp)
+ bhi, blo, round, sticky := shr128gs(b.hi, b.lo, diff)
+ slo, c := bits.Add64(a.lo, blo, 0)
+ shi, c2 := bits.Add64(a.hi, bhi, c)
+ exp := a.exp
+ if c2 != 0 { // carry into bit 128: shift right 1, recompute round/sticky
+ sticky = sticky || round
+ round = slo&1 != 0
+ slo = slo>>1 | shi<<63
+ shi = shi>>1 | 1<<63
+ exp++
+ }
+ return roundNE(shi, slo, round, sticky, exp)
+}
+
+func (a f128) cmp(b f128) int {
+ az, bz := a.isZero(), b.isZero()
+ switch {
+ case az && bz:
+ return 0
+ case az:
+ return -1
+ case bz:
+ return 1
+ case a.exp != b.exp:
+ if a.exp < b.exp {
+ return -1
+ }
+ return 1
+ case a.hi != b.hi:
+ if a.hi < b.hi {
+ return -1
+ }
+ return 1
+ case a.lo != b.lo:
+ if a.lo < b.lo {
+ return -1
+ }
+ return 1
+ }
+ return 0
+}
+
+// intPow returns base^e by exponentiation by squaring (integer exponent).
+func (base f128) intPow(e uint64) f128 {
+ result := f128FromUint64(1)
+ b := base
+ for e > 0 {
+ if e&1 == 1 {
+ result = result.mul(b)
+ }
+ e >>= 1
+ if e > 0 {
+ b = b.mul(b)
+ }
+ }
+ return result
+}
+
+// binomialF128 evaluates the CDF of Binomial(money trials, success probability p)
+// in software f128 -- the counterpart of boost::math::binomial_distribution.
+// cdf(j) returns P(X <= j). Boost computes that as ibetac(j+1, n-j, p); this
+// accumulates the IDENTICAL value as the running sum of the binomial PMF:
+//
+// pmf(0) = (1-p)^money
+// pmf(j) = pmf(j-1) * (money-j+1)/j * p/(1-p)
+// cdf(j) = pmf(0) + pmf(1) + ... + pmf(j)
+//
+// cdf MUST be called with j = 0, 1, 2, ... in increasing order (as the walk
+// below does); each call advances the running PMF/CDF using only f128 software
+// arithmetic, so cdf(j) is bit-reproducible on every platform.
+type binomialF128 struct {
+ money uint64
+ pq f128 // p/(1-p), the per-step PMF multiplier
+ pmf f128 // pmf(at): the current term
+ cum f128 // cdf(at) = P(X <= at)
+ at uint64 // index that pmf/cum currently hold
+
+ frozen bool // cum can never change again: cdf(k) == cum for every k >= at
+}
+
+// newBinomialF128 constructs the CDF evaluator for Binomial(money trials,
+// p = expectedSize/totalMoney) -- the analogue of constructing
+// binomial_distribution(n=money, p). Taking p as its exact integer
+// numerator and denominator lets each PMF-recurrence constant be a SINGLE
+// round-to-nearest-even f128 divide of exact values:
+//
+// 1-p = (totalMoney-expectedSize) / totalMoney
+// p/(1-p) = expectedSize / (totalMoney-expectedSize)
+//
+// rather than stacking roundings through an intermediate float64 p (whose
+// float64(totalMoney) conversion is itself inexact above 2^53). Note that
+// intPow then amplifies qf's single rounding by up to the trial count, so
+// pmf(0) carries ~money*2^-129 of relative error; the SelectF128 tail-edge
+// note documents the resulting CDF plateau and how the walk's freeze
+// detection handles it. Returns nil for the degenerate p >= 1 --
+// expectedSize >= totalMoney, an exact integer comparison; all probability
+// mass at j == money -- which the caller handles. totalMoney == 0 also lands
+// on the nil path.
+func newBinomialF128(expectedSize, totalMoney, money uint64) *binomialF128 {
+ if expectedSize >= totalMoney { // p >= 1
+ return nil
+ }
+ qf := f128FromUint64(totalMoney - expectedSize).div(f128FromUint64(totalMoney)) // 1-p
+ pq := f128FromUint64(expectedSize).div(f128FromUint64(totalMoney - expectedSize)) // p/(1-p)
+ pmf0 := qf.intPow(money) // (1-p)^money
+ return &binomialF128{money: money, pq: pq, pmf: pmf0, cum: pmf0, at: 0}
+}
+
+func (b *binomialF128) cdf(j uint64) f128 {
+ for b.at < j && !b.frozen {
+ b.at++
+ // pmf(at) = pmf(at-1) * (money-at+1)/at * p/(1-p)
+ pmfPrev := b.pmf
+ b.pmf = f128FromUint64(b.money - b.at + 1).divU(b.at).mul(b.pq).mul(b.pmf)
+ cumPrev := b.cum
+ b.cum = b.cum.add(b.pmf)
+ // cum is frozen once an add is a no-op while pmf strictly shrank: a
+ // shrinking pmf proves the rounded step factor is < 1, and the factor
+ // only decreases with at (round-to-nearest is monotone), so every later
+ // pmf is <= this one; and if adding THIS pmf could not move cum, no
+ // later, smaller pmf can either. From here cdf(k) == cum for all k.
+ if b.cum.cmp(cumPrev) == 0 && b.pmf.cmp(pmfPrev) < 0 {
+ b.frozen = true
+ }
+ }
+ return b.cum
+}
+
+// binomialCDFWalkF128 is the pure-Go, deterministic counterpart of the C++
+// sortition_binomial_cdf_walk in sortition.cpp. It performs the same boundary
+// walk, using an f128 ratio and f128 CDF values:
+//
+// C++ sortition.cpp Go this function
+// ------------------------------------------ -------------------------------------------
+// uint64_t sortition_binomial_cdf_walk( func binomialCDFWalkF128(
+// double n, double p, double ratio, expectedSize, totalMoney uint64,
+// uint64_t money) { ratio f128, money uint64) uint64 {
+// binomial_distribution dist(n, p); dist := newBinomialF128(expectedSize, totalMoney, money)
+// for (uint64_t j = 0; j < money; j++) { for j := uint64(0); j < money; j++ {
+// double boundary = cdf(dist, j); boundary := dist.cdf(j)
+// if (ratio <= boundary) { if ratio.cmp(boundary) <= 0 {
+// return j; return j
+// } }
+// } }
+// return money; return money
+// } }
+//
+// Boost computes cdf(dist, j) = ibetac(j+1, n-j, p) afresh each step in hardware
+// double, whereas dist.cdf(j) returns the same mathematical value as a running
+// PMF sum in software f128 (see binomialF128). The f128 path also receives the
+// digest ratio directly at f128 precision, and the success probability as its
+// exact integer numerator and denominator rather than a float64 quotient.
+//
+// Precondition: money < SelectF128MaxMoney (2^56). Below that bound no int64
+// exponent arithmetic in the walk can wrap, even at the most extreme
+// representable probability -- see the constant's doc for the accounting,
+// which must include the mantissa normalization offset (stored exp is
+// log2(v) - 127) and mul's intermediate exponent sums. Behavior beyond the
+// bound is undefined (Boost's Select cannot evaluate such money either).
+func binomialCDFWalkF128(expectedSize, totalMoney uint64, ratio f128, money uint64) uint64 {
+ dist := newBinomialF128(expectedSize, totalMoney, money)
+ if dist == nil { // p >= 1: cdf(j)==0 for j.
+
+package sortition
+
+import (
+ "encoding/hex"
+ "encoding/json"
+ "math/big"
+ "os"
+ "testing"
+)
+
+type arbDyadicCertificate struct {
+ Mantissa string `json:"mantissa"`
+ Exponent int64 `json:"exponent"`
+}
+
+type arbCDFVector struct {
+ Label string `json:"label"`
+ Money uint64 `json:"money"`
+ TotalMoney uint64 `json:"total_money"`
+ ExpectedSize uint64 `json:"expected_size"`
+ Digest string `json:"digest"`
+ Want uint64 `json:"want"`
+ PrecisionBits uint64 `json:"precision_bits"`
+ CDFPreviousUpper arbDyadicCertificate `json:"cdf_previous_upper"`
+ CDFSelectedLower arbDyadicCertificate `json:"cdf_selected_lower"`
+}
+
+type arbCertificateFile struct {
+ Format uint64 `json:"format"`
+ Generator string `json:"generator"`
+ Semantics string `json:"semantics"`
+ Vectors []arbCDFVector `json:"vectors"`
+}
+
+func parseCertificateMantissa(t *testing.T, d arbDyadicCertificate) *big.Int {
+ t.Helper()
+ m, ok := new(big.Int).SetString(d.Mantissa, 10)
+ if !ok {
+ t.Fatalf("invalid certificate mantissa %q", d.Mantissa)
+ }
+ return m
+}
+
+// compareDyadicToRat compares mantissa*2^exponent with num/den using exact
+// integer arithmetic. The certificate endpoints themselves come from Arb;
+// this check makes the inequalities carried by the JSON independently
+// inspectable by the Go test.
+func compareDyadicToRat(t *testing.T, d arbDyadicCertificate, num, den *big.Int) int {
+ t.Helper()
+ m := parseCertificateMantissa(t, d)
+ left, right := new(big.Int).Set(m), new(big.Int).Set(num)
+ if d.Exponent >= 0 {
+ left.Lsh(left, uint(d.Exponent))
+ left.Mul(left, den)
+ } else {
+ left.Mul(left, den)
+ right.Lsh(right, uint(-d.Exponent))
+ }
+ return left.Cmp(right)
+}
+
+// TestSelectF128ArbCertificates consumes offline certificates generated by
+// python-flint/Arb's regularized incomplete beta implementation. Unlike all
+// in-process floating references, Arb neither initializes at q^money nor uses
+// SelectF128's forward PMF recurrence. Every vector first proves the exact
+// binomial quantile inequality, then pins SelectF128 to that certified count.
+func TestSelectF128ArbCertificates(t *testing.T) {
+ raw, err := os.ReadFile("testdata/f128_arb_certificates.json")
+ if err != nil {
+ t.Fatal(err)
+ }
+ var certificates arbCertificateFile
+ if err := json.Unmarshal(raw, &certificates); err != nil {
+ t.Fatal(err)
+ }
+ if certificates.Format != 1 || certificates.Generator == "" || certificates.Semantics == "" || len(certificates.Vectors) < 10 {
+ t.Fatalf("unexpected Arb certificate corpus: format=%d generator=%q semantics=%q vectors=%d",
+ certificates.Format, certificates.Generator, certificates.Semantics, len(certificates.Vectors))
+ }
+
+ digestDen := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ labels := make(map[string]struct{}, len(certificates.Vectors))
+ for _, vector := range certificates.Vectors {
+ t.Run(vector.Label, func(t *testing.T) {
+ if _, exists := labels[vector.Label]; exists || vector.Label == "" {
+ t.Fatalf("empty or duplicate certificate label %q", vector.Label)
+ }
+ labels[vector.Label] = struct{}{}
+ if vector.Money >= SelectF128MaxMoney || vector.TotalMoney == 0 || vector.ExpectedSize == 0 ||
+ vector.ExpectedSize >= vector.TotalMoney || vector.Want > vector.Money || vector.PrecisionBits < 128 {
+ t.Fatalf("invalid certified parameters or result: money=%d total=%d expected=%d want=%d precision=%d",
+ vector.Money, vector.TotalMoney, vector.ExpectedSize, vector.Want, vector.PrecisionBits)
+ }
+ digestBytes, err := hex.DecodeString(vector.Digest)
+ if err != nil || len(digestBytes) != DigestSize {
+ t.Fatalf("invalid digest %q", vector.Digest)
+ }
+ var d Digest
+ copy(d[:], digestBytes)
+ digestInt := new(big.Int).SetBytes(d[:])
+ if digestInt.Sign() <= 0 {
+ t.Fatal("certificate corpus deliberately excludes the ratio==0 degenerate boundary")
+ }
+ zero, one := new(big.Int), big.NewInt(1)
+ for name, endpoint := range map[string]arbDyadicCertificate{
+ "CDF(want-1) upper": vector.CDFPreviousUpper,
+ "CDF(want) lower": vector.CDFSelectedLower,
+ } {
+ if compareDyadicToRat(t, endpoint, zero, one) < 0 ||
+ compareDyadicToRat(t, endpoint, one, one) > 0 {
+ t.Fatalf("%s endpoint is outside [0,1]", name)
+ }
+ }
+
+ // Arb supplied an upper enclosure for CDF(want-1) and a lower
+ // enclosure for CDF(want). Exact integer comparison with the
+ // digest ratio proves the selected quantile.
+ if compareDyadicToRat(t, vector.CDFPreviousUpper, digestInt, digestDen) >= 0 {
+ t.Fatal("certificate does not prove CDF(want-1) < ratio")
+ }
+ if compareDyadicToRat(t, vector.CDFSelectedLower, digestInt, digestDen) < 0 {
+ t.Fatal("certificate does not prove ratio <= CDF(want)")
+ }
+ if inFrozenSliver(vector.Money, digestRatioBig(d, 512)) {
+ t.Fatal("exact-math certificate unexpectedly lies in the separately defined frozen sliver")
+ }
+ if got := SelectF128(vector.Money, vector.TotalMoney, vector.ExpectedSize, d); got != vector.Want {
+ t.Fatalf("SelectF128=%d, Arb-certified exact binomial quantile=%d", got, vector.Want)
+ }
+ })
+ }
+}
diff --git a/f128_exact_test.go b/f128_exact_test.go
new file mode 100644
index 0000000..d8fa376
--- /dev/null
+++ b/f128_exact_test.go
@@ -0,0 +1,404 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "math/big"
+ "math/rand"
+ "testing"
+)
+
+// The differential tests compare SelectF128 against an oracle that computes
+// at the SAME 128-bit precision, so a defect mirrored into the oracle (as the
+// pmf(0) plateau was) is invisible to them. The tests in this file compare
+// against a 512-bit reference instead: the 128-bit walk must land within one
+// boundary of the near-exact answer everywhere outside the documented
+// frozen-tail sliver, and must agree with the oracle bit-for-bit at
+// deliberately constructed knife-edge digests.
+
+// selectAtPrecision runs the binomial-CDF walk at the requested big.Float
+// precision. Multiple precisions are cross-checked below so the 512-bit
+// reference is not trusted merely because it has a larger mantissa.
+func selectAtPrecision(money, totalMoney, expectedSize uint64, vrfOutput Digest, prec uint) uint64 {
+ ratio := digestRatioBig(vrfOutput, prec)
+ if expectedSize >= totalMoney { // p >= 1
+ if ratio.Sign() <= 0 {
+ return 0
+ }
+ return money
+ }
+ q := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney))
+ pq := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize))
+ pmf := bigIntPow(q, money, prec)
+ cdf := new(big.Float).SetPrec(prec).Set(pmf)
+ if cdf.Cmp(ratio) >= 0 {
+ return 0
+ }
+ for j := uint64(1); j < money && j <= testOracleMaxCDFSteps; j++ {
+ factor := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(money-j+1),
+ new(big.Float).SetPrec(prec).SetUint64(j))
+ pmf = new(big.Float).SetPrec(prec).Mul(pmf, new(big.Float).SetPrec(prec).Mul(factor, pq))
+ cdf = new(big.Float).SetPrec(prec).Add(cdf, pmf)
+ if cdf.Cmp(ratio) >= 0 {
+ return j
+ }
+ }
+ if money > testOracleMaxCDFSteps {
+ panic("selectAtPrecision exceeded the test oracle step budget")
+ }
+ return money
+}
+
+// selectHighPrec runs the tolerance oracle at its ordinary 512-bit precision.
+func selectHighPrec(money, totalMoney, expectedSize uint64, vrfOutput Digest) uint64 {
+ return selectAtPrecision(money, totalMoney, expectedSize, vrfOutput, 512)
+}
+
+// inFrozenSliver reports whether the digest ratio is within the carved-out
+// near-1.0 region where the 128-bit walk's answer is DEFINED as money (see
+// the SelectF128 doc comment) and tolerance against exact math does not
+// apply. The bound (money+2)*2^-122 covers the plateau's ~money*2^-129 with
+// two orders of margin, including the boundary-crowding zone just above it.
+func inFrozenSliver(money uint64, ratio *big.Float) bool {
+ gap := new(big.Float).SetPrec(512).Sub(new(big.Float).SetPrec(512).SetInt64(1), ratio)
+ bound := new(big.Float).SetMantExp(new(big.Float).SetPrec(64).SetUint64(money+2), -122)
+ return gap.Cmp(bound) < 0
+}
+
+// TestSelectF128NearExactMath asserts |SelectF128 - selectHighPrec| <= 1 for
+// every input outside the frozen sliver. A mirrored magnitude error like the
+// pmf(0) plateau cannot hide from this: with the pre-freeze-fix code the
+// supply-scale spots below would have diverged by ~2e15 (or hung without the
+// short-circuit). Nothing else on this branch compares supply-scale results
+// against any reference outside the sliver.
+func TestSelectF128NearExactMath(t *testing.T) {
+ rng := rand.New(rand.NewSource(6))
+ // mustAssert guards against the comparison being skipped by the sliver
+ // carve-out: the deliberately chosen spot cases must actually assert (an
+ // earlier version of the carve-out bound was accidentally so wide that
+ // the supply-scale spots passed vacuously).
+ check := func(money, total, expected uint64, d Digest, mustAssert bool) {
+ t.Helper()
+ if inFrozenSliver(money, digestRatioBig(d, 512)) {
+ if mustAssert {
+ t.Fatalf("spot case unexpectedly in the frozen sliver (money=%d vrf=%x)", money, d)
+ }
+ return
+ }
+ j128 := SelectF128(money, total, expected, d)
+ jhp := selectHighPrec(money, total, expected, d)
+ diff := j128 - jhp
+ if jhp > j128 {
+ diff = jhp - j128
+ }
+ if diff > 1 {
+ t.Fatalf("SelectF128=%d vs 512-bit reference=%d (money=%d total=%d expected=%d vrf=%x)",
+ j128, jhp, money, total, expected, d)
+ }
+ }
+
+ for i := 0; i < 400; i++ {
+ money := rng.Uint64() % 2000
+ total := rng.Uint64() >> uint(rng.Intn(40))
+ if total == 0 {
+ total = 1
+ }
+ expected := rng.Uint64() >> uint(rng.Intn(50))
+ var d Digest
+ if i%3 == 0 {
+ rng.Read(d[:])
+ } else {
+ // near-maximum digests stress the tail, where boundaries crowd
+ for j := range d {
+ d[j] = 0xff
+ }
+ d[rng.Intn(14)] = byte(rng.Uint64())
+ }
+ check(money, total, expected, d, false)
+ }
+
+ // supply-scale spots, chosen outside the sliver (gaps 2^-56..2^-66 vs a
+ // sliver bound near 2^-68 at the supply ceiling)
+ const online = uint64(2_000_000_000_000_000)
+ const supply = uint64(10_000_000_000_000_000)
+ check(online, online, 1500, maxDigestMinusPowerOfTwo(196), true) // ratio ~1-2^-60
+ check(online, online, 6000, maxDigestMinusPowerOfTwo(200), true) // ratio ~1-2^-56
+ check(supply, supply, 5000, maxDigestMinusPowerOfTwo(190), true) // ratio ~1-2^-66
+ check(supply/3, supply, 2990, maxDigestMinusPowerOfTwo(198), true) // ratio ~1-2^-58
+}
+
+// exactCDFNum returns the numerator of the EXACT binomial CDF at j over the
+// denominator totalMoney^money: sum_{i<=j} C(money,i) * E^i * (T-E)^(money-i),
+// built from first principles with stdlib binomial coefficients. It shares no
+// formula with the walk or the differential oracle, which both use the PMF
+// recurrence pmf(j) = pmf(j-1)*(money-j+1)/j * pq -- a shared algebra error
+// there would pass every differential test but not this.
+func exactCDFNum(money, totalMoney, expectedSize, j uint64) *big.Int {
+ e := new(big.Int).SetUint64(expectedSize)
+ q := new(big.Int).SetUint64(totalMoney - expectedSize)
+ sum := new(big.Int)
+ for i := uint64(0); i <= j; i++ {
+ term := new(big.Int).Binomial(int64(money), int64(i))
+ term.Mul(term, new(big.Int).Exp(e, new(big.Int).SetUint64(i), nil))
+ term.Mul(term, new(big.Int).Exp(q, new(big.Int).SetUint64(money-i), nil))
+ sum.Add(sum, term)
+ }
+ return sum
+}
+
+// selectExactRat runs the walk against the exact CDF with pure integer
+// arithmetic: ratio <= cdf(j) becomes t * T^money <= cdfNum(j) * (2^256-1).
+// No rounding anywhere. Cost grows as T^money-sized integers, so callers keep
+// money small.
+func selectExactRat(money, totalMoney, expectedSize uint64, vrfOutput Digest) uint64 {
+ tDig := new(big.Int).SetBytes(vrfOutput[:])
+ if expectedSize >= totalMoney { // p >= 1
+ if tDig.Sign() == 0 {
+ return 0
+ }
+ return money
+ }
+ den256 := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ lhs := new(big.Int).Mul(tDig,
+ new(big.Int).Exp(new(big.Int).SetUint64(totalMoney), new(big.Int).SetUint64(money), nil))
+ for j := uint64(0); j < money; j++ {
+ rhs := new(big.Int).Mul(exactCDFNum(money, totalMoney, expectedSize, j), den256)
+ if lhs.Cmp(rhs) <= 0 {
+ return j
+ }
+ }
+ return money
+}
+
+// TestSelectF128VsExactRat compares the 128-bit walk against exact
+// mathematics for small money, including digests constructed to sit exactly
+// on true CDF boundaries. Agreement must be within one boundary outside the
+// frozen sliver; a formula error shared by the implementation and the
+// differential oracle cannot hide here.
+func TestSelectF128VsExactRat(t *testing.T) {
+ rng := rand.New(rand.NewSource(9))
+ den256 := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ check := func(money, total, expected uint64, d Digest) {
+ t.Helper()
+ if inFrozenSliver(money, digestRatioBig(d, 512)) {
+ return
+ }
+ got := SelectF128(money, total, expected, d)
+ want := selectExactRat(money, total, expected, d)
+ diff := got - want
+ if want > got {
+ diff = want - got
+ }
+ if diff > 1 {
+ t.Fatalf("SelectF128=%d vs exact=%d (money=%d total=%d expected=%d vrf=%x)",
+ got, want, money, total, expected, d)
+ }
+ }
+
+ for i := 0; i < 300; i++ {
+ money := rng.Uint64() % 51
+ total := 1 + rng.Uint64()>>uint(rng.Intn(40))
+ expected := rng.Uint64() >> uint(rng.Intn(50))
+ var d Digest
+ switch {
+ case i%3 == 0:
+ rng.Read(d[:])
+ case i%3 == 1: // near-maximum digests
+ for j := range d {
+ d[j] = 0xff
+ }
+ d[rng.Intn(10)] = byte(rng.Uint64())
+ default: // a digest exactly on (or one off) a true CDF boundary
+ if money == 0 || expected >= total {
+ rng.Read(d[:])
+ break
+ }
+ j := rng.Uint64() % money
+ tt := new(big.Int).Mul(exactCDFNum(money, total, expected, j), den256)
+ tt.Div(tt, new(big.Int).Exp(new(big.Int).SetUint64(total), new(big.Int).SetUint64(money), nil))
+ tt.Add(tt, big.NewInt(int64(rng.Intn(3)-1)))
+ if tt.Sign() < 0 || tt.Cmp(den256) > 0 {
+ rng.Read(d[:])
+ break
+ }
+ tt.FillBytes(d[:])
+ }
+ check(money, total, expected, d)
+ }
+}
+
+// TestSelectF128ExactRatInclusiveBoundary constructs a digest ratio exactly
+// equal to CDF(0): for Binomial(1, 1/3), CDF(0)=2/3 and 3 divides 2^256-1.
+// Sampling cannot distinguish <= from < at this measure-zero boundary.
+func TestSelectF128ExactRatInclusiveBoundary(t *testing.T) {
+ den := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ at := new(big.Int).Quo(new(big.Int).Mul(new(big.Int).Set(den), big.NewInt(2)), big.NewInt(3))
+ if new(big.Int).Mod(new(big.Int).Set(den), big.NewInt(3)).Sign() != 0 {
+ t.Fatal("test construction requires 3 to divide 2^256-1")
+ }
+ for offset, want := range map[int64]uint64{-1: 0, 0: 0, 1: 1} {
+ point := new(big.Int).Add(new(big.Int).Set(at), big.NewInt(offset))
+ d := digestFromBigInt(t, point, den)
+ if got := selectExactRat(1, 3, 1, d); got != want {
+ t.Fatalf("offset %d: exact selector=%d, want %d", offset, got, want)
+ }
+ }
+}
+
+// TestSelectF128Distribution mirrors TestSortitionBasic for the pure-Go path:
+// summed selection weight over many digests must track N * money * p. This is
+// oracle-free AND formula-free -- a gross semantic error (wrong probability,
+// shifted distribution) fails it even if perfectly mirrored everywhere else.
+func TestSelectF128Distribution(t *testing.T) {
+ rng := rand.New(rand.NewSource(4))
+ cases := []struct {
+ name string
+ money, total, expected, rounds uint64
+ }{
+ {"half stake", 100, 200, 20, 1000},
+ {"realistic stake", 1_000_000_000_000_000, 2_000_000_000_000_000, 1500, 200},
+ }
+ for _, c := range cases {
+ var hits uint64
+ for i := uint64(0); i < c.rounds; i++ {
+ var d Digest
+ rng.Read(d[:])
+ hits += SelectF128(c.money, c.total, c.expected, d)
+ }
+ want := float64(c.rounds) * float64(c.expected) * float64(c.money) / float64(c.total)
+ if diff := float64(hits) - want; diff < -0.02*want || diff > 0.02*want {
+ t.Errorf("%s: %d selections over %d rounds, want %.0f +/- 2%%",
+ c.name, hits, c.rounds, want)
+ }
+ }
+}
+
+// oracleCDFAt returns the oracle's cdf(j), mirroring selectBigOracle's
+// 128-bit recurrence exactly.
+func oracleCDFAt(money, totalMoney, expectedSize, j uint64) *big.Float {
+ const prec = f128MantBits
+ q := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney))
+ pq := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize))
+ pmf := bigIntPow(q, money, prec)
+ cdf := new(big.Float).SetPrec(prec).Set(pmf)
+ for i := uint64(1); i <= j; i++ {
+ factor := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(money-i+1),
+ new(big.Float).SetPrec(prec).SetUint64(i))
+ step := new(big.Float).SetPrec(prec).Mul(factor, pq)
+ pmf = new(big.Float).SetPrec(prec).Mul(pmf, step)
+ cdf = new(big.Float).SetPrec(prec).Add(cdf, pmf)
+ }
+ return cdf
+}
+
+// TestSelectF128BoundaryStraddle constructs digests one ratio-ulp below, at,
+// and above oracle CDF boundaries -- the knife edges where rounding bugs
+// live, which uniform sampling hits with probability ~2^-128 -- and asserts
+// bit-level agreement with the oracle at each, plus monotone results across
+// the constructed digests.
+func TestSelectF128BoundaryStraddle(t *testing.T) {
+ rng := rand.New(rand.NewSource(8))
+ den := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ denF := new(big.Float).SetPrec(600).SetInt(den)
+ var windows, bracketed, clampedWindows, candidates uint64
+
+ for iter := 0; iter < 150; iter++ {
+ money := 2 + rng.Uint64()%250
+ total := 2 + rng.Uint64()>>uint(rng.Intn(40))
+ expected := rng.Uint64() % total // p < 1
+ boundaries := []uint64{0, money / 2, money - 1, rng.Uint64() % money}
+ for _, j := range boundaries {
+ windows++
+ c := oracleCDFAt(money, total, expected, j)
+ // One f128 ulp at c's own magnitude, in digest units: the ulp is
+ // 2^(exp-128) for c in [2^(exp-1), 2^exp), so it scales with the
+ // boundary -- a fixed 2^128 step would be one ulp only for c in
+ // [0.5, 1) and many ulps for small boundaries like a large-mean
+ // cdf(0). Clamp to one digest unit when the ulp is finer than the
+ // digest grid.
+ step := new(big.Int).Rsh(den, uint(128-c.MantExp(nil)))
+ if step.Sign() == 0 {
+ step.SetInt64(1)
+ }
+ ulp := new(big.Float).SetMantExp(new(big.Float).SetPrec(600).SetInt64(1), c.MantExp(nil)-128)
+ wantStep, _ := new(big.Float).SetPrec(600).Mul(ulp, denF).Int(nil)
+ if wantStep.Sign() == 0 {
+ wantStep.SetInt64(1)
+ }
+ if step.Cmp(wantStep) != 0 {
+ t.Fatalf("straddle step %s is not one f128 ulp (%s) at boundary j=%d", step, wantStep, j)
+ }
+ tt, _ := new(big.Float).SetPrec(600).Mul(c, denF).Int(nil)
+ prev := uint64(0)
+ first := true
+ clamped := false
+ var loRatio, hiRatio *big.Float
+ for k := int64(-2); k <= 2; k++ {
+ ti := new(big.Int).Add(tt, new(big.Int).Mul(big.NewInt(k), step))
+ if ti.Sign() < 0 || ti.Cmp(den) > 0 {
+ clamped = true
+ continue
+ }
+ candidates++
+ var d Digest
+ ti.FillBytes(d[:])
+ got := SelectF128(money, total, expected, d)
+ want := selectBigOracle(money, total, expected, d)
+ if got != want {
+ t.Fatalf("straddle mismatch at boundary j=%d offset %d ulp: SelectF128=%d oracle=%d (money=%d total=%d expected=%d vrf=%x)",
+ j, k, got, want, money, total, expected, d)
+ }
+ if !first && got < prev {
+ t.Fatalf("non-monotone across straddle at boundary j=%d: %d then %d (money=%d total=%d expected=%d)",
+ j, prev, got, money, total, expected)
+ }
+ prev, first = got, false
+ if loRatio == nil {
+ loRatio = digestRatioBig(d, f128MantBits)
+ }
+ hiRatio = digestRatioBig(d, f128MantBits)
+ }
+ // The candidates must genuinely bracket the boundary (this is
+ // what makes it a straddle rather than ordinary sampling); only
+ // checkable when no candidate was clamped away at 0 or the
+ // maximum digest.
+ if !clamped && (loRatio.Cmp(c) > 0 || hiRatio.Cmp(c) < 0) {
+ t.Fatalf("straddle window does not bracket boundary j=%d: [%v, %v] vs cdf=%v (money=%d total=%d expected=%d)",
+ j, loRatio, hiRatio, c, money, total, expected)
+ }
+ if clamped {
+ clampedWindows++
+ } else {
+ bracketed++
+ }
+ }
+ }
+ if windows != 600 || bracketed < 400 || candidates < 2_000 || bracketed+clampedWindows != windows {
+ t.Fatalf("anti-vacuity: windows=%d bracketed=%d clamped=%d candidates=%d",
+ windows, bracketed, clampedWindows, candidates)
+ }
+}
diff --git a/f128_integer_oracle_test.go b/f128_integer_oracle_test.go
new file mode 100644
index 0000000..6f42ec2
--- /dev/null
+++ b/f128_integer_oracle_test.go
@@ -0,0 +1,282 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "math/big"
+ "math/rand"
+ "testing"
+)
+
+var integerOracleOne = big.NewInt(1)
+
+// floorLog2Rat returns floor(log2(num/den)) for positive integers. It uses
+// integer comparison only; in particular, the exact-integer oracle below does
+// not ask big.Float to decide normalization or rounding.
+func floorLog2Rat(num, den *big.Int) int {
+ e := num.BitLen() - den.BitLen()
+ if e >= 0 {
+ if num.Cmp(new(big.Int).Lsh(new(big.Int).Set(den), uint(e))) < 0 {
+ e--
+ }
+ } else if new(big.Int).Lsh(new(big.Int).Set(num), uint(-e)).Cmp(den) < 0 {
+ e--
+ }
+ return e
+}
+
+// roundExactRatToF128 rounds (num/den)*2^exp2 to f128 using only big.Int
+// quotient/remainder arithmetic. num must be non-negative and den positive.
+func roundExactRatToF128(num, den *big.Int, exp2 int64) f128 {
+ if num.Sign() == 0 {
+ return f128{}
+ }
+ if num.Sign() < 0 || den.Sign() <= 0 {
+ panic("roundExactRatToF128 requires num >= 0 and den > 0")
+ }
+
+ e := floorLog2Rat(num, den)
+ n, d := new(big.Int).Set(num), new(big.Int).Set(den)
+ if shift := 127 - e; shift >= 0 {
+ n.Lsh(n, uint(shift))
+ } else {
+ d.Lsh(d, uint(-shift))
+ }
+
+ mant, rem := new(big.Int), new(big.Int)
+ mant.QuoRem(n, d, rem)
+ twiceRem := new(big.Int).Lsh(new(big.Int).Set(rem), 1)
+ if cmp := twiceRem.Cmp(d); cmp > 0 || (cmp == 0 && mant.Bit(0) != 0) {
+ mant.Add(mant, integerOracleOne)
+ }
+
+ outExp := exp2 + int64(e) - 127
+ if mant.BitLen() == 129 { // rounded 2^128: renormalize to 2^127
+ mant.Rsh(mant, 1)
+ outExp++
+ }
+ if mant.BitLen() != 128 {
+ panic("integer oracle produced a non-normalized mantissa")
+ }
+ lo := mant.Uint64()
+ hi := new(big.Int).Rsh(new(big.Int).Set(mant), 64).Uint64()
+ return f128{hi: hi, lo: lo, exp: outExp}
+}
+
+func f128MantissaInt(x f128) *big.Int {
+ m := new(big.Int).SetUint64(x.hi)
+ m.Lsh(m, 64)
+ return m.Add(m, new(big.Int).SetUint64(x.lo))
+}
+
+func exactMulF128(a, b f128) f128 {
+ if a.isZero() || b.isZero() {
+ return f128{}
+ }
+ n := new(big.Int).Mul(f128MantissaInt(a), f128MantissaInt(b))
+ return roundExactRatToF128(n, integerOracleOne, a.exp+b.exp)
+}
+
+func exactAddF128(a, b f128) f128 {
+ if a.isZero() {
+ return b
+ }
+ if b.isZero() {
+ return a
+ }
+ baseExp := a.exp
+ if b.exp < baseExp {
+ baseExp = b.exp
+ }
+ an := f128MantissaInt(a)
+ bn := f128MantissaInt(b)
+ an.Lsh(an, uint(a.exp-baseExp))
+ bn.Lsh(bn, uint(b.exp-baseExp))
+ return roundExactRatToF128(new(big.Int).Add(an, bn), integerOracleOne, baseExp)
+}
+
+func exactDivF128(a, b f128) f128 {
+ if a.isZero() || b.isZero() {
+ return f128{}
+ }
+ return roundExactRatToF128(f128MantissaInt(a), f128MantissaInt(b), a.exp-b.exp)
+}
+
+func exactDivUF128(a f128, u uint64) f128 {
+ if a.isZero() || u == 0 {
+ return f128{}
+ }
+ return roundExactRatToF128(f128MantissaInt(a), new(big.Int).SetUint64(u), a.exp)
+}
+
+func assertCanonicalF128(t *testing.T, name string, x f128) {
+ t.Helper()
+ if x.isZero() {
+ if x.exp != 0 {
+ t.Fatalf("%s: zero has nonzero exponent: %+v", name, x)
+ }
+ return
+ }
+ if x.hi&(uint64(1)<<63) == 0 {
+ t.Fatalf("%s: mantissa is not normalized: %+v", name, x)
+ }
+}
+
+func assertExactF128(t *testing.T, name string, got, want f128) {
+ t.Helper()
+ assertCanonicalF128(t, name, got)
+ if got != want {
+ t.Fatalf("%s: got %+v, exact-integer RNE oracle wants %+v", name, got, want)
+ }
+}
+
+// TestF128OpsExactIntegerOracle validates the public arithmetic contract with
+// an oracle structurally independent of math/big.Float. Exact values are
+// represented as integer ratios; normalization and ties-to-even are decided
+// from quotient and remainder.
+func TestF128OpsExactIntegerOracle(t *testing.T) {
+ rng := rand.New(rand.NewSource(10))
+ gaps := [...]int64{0, 1, -1, 63, -63, 64, -64, 65, -65, 127, -127, 128, -128, 129, -129, 511, -511, 2000, -2000}
+ for i := 0; i < 100_000; i++ {
+ aexp := int64(rng.Intn(4001) - 2000)
+ bexp := aexp - gaps[i%len(gaps)]
+ a := f128{hi: rng.Uint64() | 1<<63, lo: rng.Uint64(), exp: aexp}
+ b := f128{hi: rng.Uint64() | 1<<63, lo: rng.Uint64(), exp: bexp}
+ if i%997 == 0 {
+ a = f128{}
+ }
+ if i%991 == 0 {
+ b = f128{}
+ }
+
+ assertExactF128(t, "mul", a.mul(b), exactMulF128(a, b))
+ assertExactF128(t, "add", a.add(b), exactAddF128(a, b))
+ assertExactF128(t, "div", a.div(b), exactDivF128(a, b))
+
+ var u uint64
+ switch i % 8 {
+ case 0:
+ u = 0
+ case 1:
+ u = 1
+ case 2:
+ u = SelectF128MaxMoney - 1
+ case 3:
+ u = 1<<62 | rng.Uint64()&((1<<62)-1)
+ case 4:
+ u = 1<<63 | rng.Uint64()&((1<<63)-1)
+ case 5:
+ u = ^uint64(0)
+ default:
+ u = rng.Uint64()
+ }
+ assertExactF128(t, "divU", a.divU(u), exactDivUF128(a, u))
+ }
+}
+
+// TestF128ConversionsExactIntegerOracle independently checks the two exact
+// input conversions. The existing big.Float matrix remains useful; this test
+// makes its rounding oracle non-circular.
+func TestF128ConversionsExactIntegerOracle(t *testing.T) {
+ rng := rand.New(rand.NewSource(11))
+ den := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), integerOracleOne)
+ for i := 0; i < 20_000; i++ {
+ u := rng.Uint64()
+ assertExactF128(t, "fromUint64", f128FromUint64(u),
+ roundExactRatToF128(new(big.Int).SetUint64(u), integerOracleOne, 0))
+
+ var d Digest
+ _, _ = rng.Read(d[:])
+ if i == 0 {
+ d = Digest{}
+ } else if i == 1 {
+ for j := range d {
+ d[j] = 0xff
+ }
+ }
+ assertExactF128(t, "digestRatio", f128FromDigestRatio(d),
+ roundExactRatToF128(new(big.Int).SetBytes(d[:]), den, 0))
+ }
+}
+
+func words128(x *big.Int) (uint64, uint64) {
+ lo := x.Uint64()
+ hi := new(big.Int).Rsh(new(big.Int).Set(x), 64).Uint64()
+ return hi, lo
+}
+
+// TestDivStepExactCertificate checks the long-division digit directly rather
+// than only through the final rounded f128 quotient. Under divStep's prefix
+// precondition, the returned digit and remainder must be the Euclidean
+// quotient and remainder of the 192-by-128-bit division.
+func TestDivStepExactCertificate(t *testing.T) {
+ rng := rand.New(rand.NewSource(12))
+ check := func(v1, v0 uint64, prefix *big.Int, uLo uint64) {
+ t.Helper()
+ v := new(big.Int).SetUint64(v1)
+ v.Lsh(v, 64).Add(v, new(big.Int).SetUint64(v0))
+ if prefix.Sign() < 0 || prefix.Cmp(v) >= 0 {
+ t.Fatal("invalid divStep test prefix")
+ }
+ uHi, uMid := words128(prefix)
+ q, rHi, rLo := divStep(uHi, uMid, uLo, v1, v0)
+
+ u := new(big.Int).Lsh(new(big.Int).Set(prefix), 64)
+ u.Add(u, new(big.Int).SetUint64(uLo))
+ wantQ, wantR := new(big.Int), new(big.Int)
+ wantQ.QuoRem(u, v, wantR)
+ if wantQ.BitLen() > 64 || q != wantQ.Uint64() {
+ t.Fatalf("divStep quotient: got %#x, want %s", q, wantQ.Text(16))
+ }
+ wantRHi, wantRLo := words128(wantR)
+ if rHi != wantRHi || rLo != wantRLo {
+ t.Fatalf("divStep remainder: got %#x:%#x, want %#x:%#x", rHi, rLo, wantRHi, wantRLo)
+ }
+
+ // Check the certificate separately from comparison with QuoRem.
+ gotR := new(big.Int).SetUint64(rHi)
+ gotR.Lsh(gotR, 64).Add(gotR, new(big.Int).SetUint64(rLo))
+ reconstructed := new(big.Int).Mul(new(big.Int).SetUint64(q), v)
+ reconstructed.Add(reconstructed, gotR)
+ if reconstructed.Cmp(u) != 0 || gotR.Sign() < 0 || gotR.Cmp(v) >= 0 {
+ t.Fatalf("invalid division certificate: U=%x q=%x V=%x R=%x", u, q, v, gotR)
+ }
+ }
+
+ for i := 0; i < 200_000; i++ {
+ v1 := rng.Uint64() | 1<<63
+ v0 := rng.Uint64()
+ if i%16 == 0 {
+ v0 |= 1 // V-1 then has uHi == v1 and exercises the qhat cap
+ }
+ v := new(big.Int).SetUint64(v1)
+ v.Lsh(v, 64).Add(v, new(big.Int).SetUint64(v0))
+ prefix := new(big.Int).SetUint64(rng.Uint64())
+ prefix.Lsh(prefix, 64).Add(prefix, new(big.Int).SetUint64(rng.Uint64()))
+ prefix.Mod(prefix, v)
+ if i%16 == 0 {
+ prefix.Sub(v, integerOracleOne)
+ }
+ check(v1, v0, prefix, rng.Uint64())
+ }
+
+ // qhat cap with rhat carry: uHi == v1, uMid+v1 overflows, and uMid < v0.
+ v1, v0 := uint64(1)<<63, ^uint64(0)
+ prefix := new(big.Int).SetUint64(v1)
+ prefix.Lsh(prefix, 64).Add(prefix, new(big.Int).SetUint64(1<<63))
+ check(v1, v0, prefix, 0)
+}
diff --git a/f128_rapid_test.go b/f128_rapid_test.go
new file mode 100644
index 0000000..4652d29
--- /dev/null
+++ b/f128_rapid_test.go
@@ -0,0 +1,296 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "bytes"
+ "math/big"
+ "testing"
+
+ "pgregory.net/rapid"
+)
+
+// These rapid property tests complement the go fuzz targets rather than
+// replacing them. Go's fuzzer mutates locally around its corpus, so integer
+// arguments dwell near seed values -- divU's shallow-quotient misrounding
+// (u > ~2^62) survived ~30M fuzz execs with u=1 and u=7 seeds, while rapid
+// found it within its first hundred generated cases. The generators below add
+// explicit magnitude-band guidance on top of rapid's own boundary bias so
+// every regime is sampled every run. Ordinary go test runs 100 cases per
+// test; CI runs them again at -rapid.checks=20000, and a long run can go
+// deeper:
+//
+// go test -run TestRapid -rapid.checks=100000
+//
+// A failure writes a minimized reproducer under testdata/rapid/.
+
+// bandedUint64 draws a uint64 spread across magnitude bands: the sortition
+// walk's divisor domain, the beyond-domain deep-quotient region, and the two
+// shallow-quotient bands where a divU rounding bug previously hid.
+func bandedUint64(t *rapid.T, label string) uint64 {
+ return rapid.OneOf(
+ rapid.Uint64Range(0, SelectF128MaxMoney), // walk domain (u = step index)
+ rapid.Uint64Range(SelectF128MaxMoney, 1<<62), // deep quotients beyond the domain
+ rapid.Uint64Range(1<<62, 1<<63), // quotient round bit in the last digit
+ rapid.Uint64Range(1<<63, ^uint64(0)), // 128-bit shallow quotients
+ ).Draw(t, label)
+}
+
+// bandedMantissa draws mantissa words by structure, not just value: uniform
+// bits, sparse (a few set bits), or dense (all ones with a few holes). Sparse
+// and dense mantissas raise the density of exact ties and of carry-out
+// renormalization, which sit near ~2^-64 measure under uniform draws.
+func bandedMantissa(t *rapid.T, label string) (uint64, uint64) {
+ switch rapid.IntRange(0, 2).Draw(t, label+"Kind") {
+ case 0:
+ return rapid.Uint64().Draw(t, label+"Hi"), rapid.Uint64().Draw(t, label+"Lo")
+ case 1: // sparse
+ var hi, lo uint64
+ for range rapid.IntRange(1, 3).Draw(t, label+"Bits") {
+ b := rapid.IntRange(0, 127).Draw(t, label+"Bit")
+ if b >= 64 {
+ hi |= 1 << (b - 64)
+ } else {
+ lo |= 1 << b
+ }
+ }
+ return hi, lo
+ default: // dense
+ hi, lo := ^uint64(0), ^uint64(0)
+ for range rapid.IntRange(0, 3).Draw(t, label+"Holes") {
+ b := rapid.IntRange(0, 127).Draw(t, label+"Hole")
+ if b >= 64 {
+ hi &^= 1 << (b - 64)
+ } else {
+ lo &^= 1 << b
+ }
+ }
+ return hi, lo
+ }
+}
+
+// TestRapidF128Ops property-tests every f128 primitive against 128-bit
+// big.Float, mirroring FuzzF128Ops, plus the oracle-independent commutativity
+// of mul and add (asymmetric partial-product assembly cannot hide from these
+// even if a bug were somehow mirrored into the reference computation).
+func TestRapidF128Ops(t *testing.T) {
+ rapid.Check(t, func(t *rapid.T) {
+ ahi, alo := bandedMantissa(t, "a")
+ aexp := rapid.Int64Range(-2000, 2000).Draw(t, "aexp")
+ a := norm128(ahi, alo, aexp)
+ // Band the exponent gap: add's alignment boundaries sit at gaps of
+ // exactly 64, 127, 128, and 129, which a uniform pair of exponents
+ // rarely produces.
+ gap := rapid.OneOf(
+ rapid.Int64Range(-4, 4),
+ rapid.Int64Range(-68, -60), rapid.Int64Range(60, 68),
+ rapid.Int64Range(-132, -124), rapid.Int64Range(124, 132),
+ rapid.Int64Range(-2000, 2000),
+ ).Draw(t, "gap")
+ bhi, blo := bandedMantissa(t, "b")
+ b := norm128(bhi, blo, aexp-gap)
+ u := bandedUint64(t, "u")
+ ab, bb := f128ToBig(a), f128ToBig(b)
+ check := func(name string, got f128, want *big.Float) {
+ if g := f128ToBig(got); g.Cmp(want) != 0 {
+ t.Fatalf("%s: f128=%v big.Float=%v (a=%v b=%v u=%d)", name, g, want, ab, bb, u)
+ }
+ }
+ check("mul", a.mul(b), new(big.Float).SetPrec(f128MantBits).Mul(ab, bb))
+ check("add", a.add(b), new(big.Float).SetPrec(f128MantBits).Add(ab, bb))
+ if u != 0 {
+ check("divU", a.divU(u),
+ new(big.Float).SetPrec(f128MantBits).Quo(ab, new(big.Float).SetPrec(f128MantBits).SetUint64(u)))
+ }
+ if !b.isZero() {
+ check("div", a.div(b), new(big.Float).SetPrec(f128MantBits).Quo(ab, bb))
+ }
+ if x, y := a.mul(b), b.mul(a); x != y {
+ t.Fatalf("mul not commutative: %+v vs %+v (a=%v b=%v)", x, y, ab, bb)
+ }
+ if x, y := a.add(b), b.add(a); x != y {
+ t.Fatalf("add not commutative: %+v vs %+v (a=%v b=%v)", x, y, ab, bb)
+ }
+ })
+}
+
+// TestRapidSelectF128VsOracle property-tests the full walk against the
+// big.Float oracle, mirroring FuzzSelectF128. Digests are drawn both uniformly
+// and from the near-maximum regime (mostly-0xff), where the ratio ~1 tail
+// edges live.
+func TestRapidSelectF128VsOracle(t *testing.T) {
+ nearMaxDigest := rapid.Custom(func(t *rapid.T) []byte {
+ b := bytes.Repeat([]byte{0xff}, DigestSize)
+ i := rapid.IntRange(0, DigestSize-1).Draw(t, "hole")
+ b[i] = rapid.Byte().Draw(t, "holeval")
+ return b
+ })
+ rapid.Check(t, func(t *rapid.T) {
+ money := rapid.Uint64Range(0, 3000).Draw(t, "money") // bound the walk, as in FuzzSelectF128
+ total := rapid.OneOf(
+ rapid.Uint64Range(0, 1_000_000),
+ rapid.Uint64Range(1_000_000, 10_000_000_000_000_000), // through the supply ceiling
+ rapid.Uint64Range(10_000_000_000_000_000, ^uint64(0)),
+ ).Draw(t, "total")
+ expected := rapid.OneOf(
+ rapid.Uint64Range(0, 10_000), // committee sizes
+ rapid.Uint64Range(10_000, ^uint64(0)), // through and beyond p >= 1
+ ).Draw(t, "expected")
+ var d Digest
+ copy(d[:], rapid.OneOf(
+ rapid.SliceOfN(rapid.Byte(), DigestSize, DigestSize),
+ nearMaxDigest,
+ ).Draw(t, "vrf"))
+ got := SelectF128(money, total, expected, d)
+ want := selectBigOracle(money, total, expected, d)
+ if got != want {
+ t.Fatalf("SelectF128=%d != oracle=%d (money=%d total=%d expected=%d vrf=%x)",
+ got, want, money, total, expected, d)
+ }
+ })
+}
+
+// TestRapidSelectF128ScaleInvariance checks an ORACLE-INDEPENDENT exact
+// property: scaling totalMoney and expectedSize by a common power of two
+// leaves every rational in the walk identical -- the same quotients round to
+// the same 128-bit values -- so the result must be BIT-IDENTICAL. Unlike
+// tolerance-based metamorphic properties (monotonicity in p or money), this
+// cannot flake at knife edges, and it catches numerator/denominator
+// mishandling without consulting the oracle.
+func TestRapidSelectF128ScaleInvariance(t *testing.T) {
+ rapid.Check(t, func(t *rapid.T) {
+ money := rapid.Uint64Range(0, 3000).Draw(t, "money")
+ total := rapid.Uint64Range(0, 1<<40).Draw(t, "total")
+ expected := rapid.Uint64Range(0, 1<<40).Draw(t, "expected")
+ k := rapid.IntRange(1, 23).Draw(t, "k")
+ var d Digest
+ copy(d[:], rapid.SliceOfN(rapid.Byte(), DigestSize, DigestSize).Draw(t, "vrf"))
+ base := SelectF128(money, total, expected, d)
+ scaled := SelectF128(money, total< 1_000_000_000_000_000 {
+ maxBase = 1_000_000_000_000_000
+ }
+ total := rapid.Uint64Range(1, maxBase).Draw(t, "total")
+ expected := rapid.Uint64Range(0, total).Draw(t, "expected")
+ var d Digest
+ copy(d[:], rapid.SliceOfN(rapid.Byte(), DigestSize, DigestSize).Draw(t, "vrf"))
+
+ baseDist := newBinomialF128(expected, total, money)
+ scaledDist := newBinomialF128(expected*factor, total*factor, money)
+ if (baseDist == nil) != (scaledDist == nil) {
+ t.Fatalf("common factor changed degenerate classification: T=%d E=%d k=%d", total, expected, factor)
+ }
+ if baseDist != nil && (baseDist.pq != scaledDist.pq || baseDist.pmf != scaledDist.pmf || baseDist.cum != scaledDist.cum) {
+ t.Fatalf("common factor changed initialized distribution: T=%d E=%d k=%d", total, expected, factor)
+ }
+ base := SelectF128(money, total, expected, d)
+ scaled := SelectF128(money, total*factor, expected*factor, d)
+ if base != scaled {
+ t.Fatalf("common-factor variance: SelectF128=%d but factor %d gives %d (money=%d total=%d expected=%d vrf=%x)",
+ base, factor, scaled, money, total, expected, d)
+ }
+ })
+}
+
+// TestRapidF128Metamorphic checks algebraic and order contracts without any
+// floating-point oracle. These overlap the exact-integer differential test on
+// purpose: corruption of either oracle cannot make an identity or monotonicity
+// failure disappear.
+func TestRapidF128Metamorphic(t *testing.T) {
+ rapid.Check(t, func(t *rapid.T) {
+ makeValue := func(label string) f128 {
+ return f128{
+ hi: rapid.Uint64().Draw(t, label+"Hi") | 1<<63,
+ lo: rapid.Uint64().Draw(t, label+"Lo"),
+ exp: rapid.Int64Range(-2000, 2000).Draw(t, label+"Exp"),
+ }
+ }
+ a, b, c, d := makeValue("a"), makeValue("b"), makeValue("c"), makeValue("d")
+ zero, one := f128{}, f128FromUint64(1)
+ if a.add(zero) != a || a.mul(one) != a || a.divU(1) != a || a.div(one) != a {
+ t.Fatalf("primitive identity failed for a=%+v", a)
+ }
+ if got := a.div(a); got != one {
+ t.Fatalf("self division=%+v, want one=%+v", got, one)
+ }
+ shift := rapid.Int64Range(-1000, 1000).Draw(t, "powerShift")
+ power := one
+ power.exp += shift
+ wantShift := a
+ wantShift.exp += shift
+ if got := a.mul(power); got != wantShift {
+ t.Fatalf("power-of-two shift=%+v, want %+v (a=%+v shift=%d)", got, wantShift, a, shift)
+ }
+
+ if a.cmp(b) > 0 {
+ a, b = b, a
+ }
+ if a.add(c).cmp(b.add(c)) > 0 || a.mul(c).cmp(b.mul(c)) > 0 || a.div(c).cmp(b.div(c)) > 0 {
+ t.Fatalf("numerator monotonicity failed: a=%+v b=%+v c=%+v", a, b, c)
+ }
+ if c.cmp(d) > 0 {
+ c, d = d, c
+ }
+ if a.div(c).cmp(a.div(d)) < 0 {
+ t.Fatalf("division denominator antitonicity failed: a=%+v c=%+v d=%+v", a, c, d)
+ }
+ })
+}
+
+// TestRapidSelectF128DigestMonotonic checks an ORACLE-INDEPENDENT property:
+// for fixed (money, total, expected), the selection count is non-decreasing in
+// the digest. This holds exactly -- the digest-to-ratio conversion is monotone
+// (round-to-nearest of a monotone quotient, and the halfway correction only
+// ever rounds up), and a larger ratio can only cross the same CDF boundaries
+// later or freeze to money. The differential tests share one structural blind
+// spot: a defect mirrored into the big.Float oracle (as the pmf(0) plateau
+// was) is invisible to them; a property test against mathematics is not.
+func TestRapidSelectF128DigestMonotonic(t *testing.T) {
+ rapid.Check(t, func(t *rapid.T) {
+ money := rapid.Uint64Range(0, 3000).Draw(t, "money")
+ total := rapid.Uint64Range(0, 10_000_000_000_000_000).Draw(t, "total")
+ expected := rapid.Uint64Range(0, 10_000).Draw(t, "expected")
+ var d1, d2 Digest
+ copy(d1[:], rapid.SliceOfN(rapid.Byte(), DigestSize, DigestSize).Draw(t, "vrf"))
+ // OR-ing random bits into d1 yields d2 >= d1 as a 256-bit integer.
+ d2 = d1
+ for _, i := range rapid.SliceOfN(rapid.IntRange(0, DigestSize-1), 1, 8).Draw(t, "orBytes") {
+ d2[i] |= rapid.Byte().Draw(t, "orVal")
+ }
+ low, high := SelectF128(money, total, expected, d1), SelectF128(money, total, expected, d2)
+ if low > high {
+ t.Fatalf("monotonicity violated: SelectF128(d1)=%d > SelectF128(d2)=%d with d1 <= d2 (money=%d total=%d expected=%d d1=%x d2=%x)",
+ low, high, money, total, expected, d1, d2)
+ }
+ })
+}
diff --git a/f128_test.go b/f128_test.go
new file mode 100644
index 0000000..2b00ee2
--- /dev/null
+++ b/f128_test.go
@@ -0,0 +1,611 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "bytes"
+ "math/big"
+ "math/rand"
+ "testing"
+)
+
+// Test oracles only exercise committee-scale quantiles at large money. A
+// broken oracle must fail instead of attempting a supply-sized fall-through;
+// small-money exhaustive and fuzz cases remain allowed to walk to money.
+const testOracleMaxCDFSteps = uint64(20_000)
+
+// selectBigOracle is an independent, "obviously correct" reference for SelectF128:
+// the SAME binomial-CDF walk, but every arithmetic step uses math/big.Float (Go's
+// standard arbitrary-precision float, round-to-nearest-even) at the f128 mantissa
+// width instead of the hand-rolled f128. SelectF128 is fuzz-checked to be
+// bit-identical to this (FuzzSelectF128), which is what justifies trusting the
+// hand-rolled integer arithmetic in f128.go. Mirrors binomialCDFWalkF128 exactly.
+//
+// Valid range: big.Float's exponent is an int32, so pmf0 = (1-p)^money underflows
+// to exactly 0 once it drops below ~2^(-2^31). That requires the mean money*p to
+// exceed ~1.5e9 (a committee of ~1.5 billion) -- far outside any reachable
+// sortition input, where the mean is a committee size (<= a few thousand). So the
+// oracle is exact for every reachable input: FuzzSelectF128 stays well inside it
+// (money %= 3001), and TestF128MatchesOracleLargeMoney exercises realistic large
+// money (up to ~2^51) where it is likewise exact. (A truly unbounded oracle would
+// need big.Rat, which is infeasible at large money -- (1-p)^money has a
+// total^money denominator -- so this range, covering all reachable inputs, is the
+// practical maximum.)
+func selectBigOracle(money uint64, totalMoney uint64, expectedSize uint64, vrfOutput Digest) uint64 {
+ const prec = f128MantBits
+ ratio := digestRatioBig(vrfOutput, prec)
+ if expectedSize >= totalMoney { // p >= 1
+ if ratio.Sign() <= 0 {
+ return 0
+ }
+ return money
+ }
+ q := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney))
+ pq := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(expectedSize),
+ new(big.Float).SetPrec(prec).SetUint64(totalMoney-expectedSize))
+ pmf := bigIntPow(q, money, prec) // (1-p)^money
+ cdf := new(big.Float).SetPrec(prec).Set(pmf)
+ if cdf.Cmp(ratio) >= 0 {
+ return 0
+ }
+ for j := uint64(1); j < money && j <= testOracleMaxCDFSteps; j++ {
+ factor := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(money-j+1),
+ new(big.Float).SetPrec(prec).SetUint64(j))
+ step := new(big.Float).SetPrec(prec).Mul(factor, pq)
+ pmf = new(big.Float).SetPrec(prec).Mul(pmf, step)
+ cdf = new(big.Float).SetPrec(prec).Add(cdf, pmf)
+ if cdf.Cmp(ratio) >= 0 {
+ return j
+ }
+ }
+ if money > testOracleMaxCDFSteps {
+ panic("selectBigOracle exceeded the test oracle step budget")
+ }
+ return money
+}
+
+func digestRatioBig(d Digest, prec uint) *big.Float {
+ numerator := new(big.Int).SetBytes(d[:])
+ denominator := new(big.Int).Lsh(big.NewInt(1), DigestSize*8)
+ denominator.Sub(denominator, big.NewInt(1))
+ return new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetInt(numerator),
+ new(big.Float).SetInt(denominator),
+ )
+}
+
+func bigIntPow(base *big.Float, e uint64, prec uint) *big.Float {
+ result := new(big.Float).SetPrec(prec).SetInt64(1)
+ b := new(big.Float).SetPrec(prec).Set(base)
+ for e > 0 {
+ if e&1 == 1 {
+ result = new(big.Float).SetPrec(prec).Mul(result, b)
+ }
+ e >>= 1
+ if e > 0 {
+ b = new(big.Float).SetPrec(prec).Mul(b, b)
+ }
+ }
+ return result
+}
+
+// FuzzSelectF128 differentially fuzzes the hand-rolled f128 implementation
+// against the independent math/big.Float oracle.
+func FuzzSelectF128(f *testing.F) {
+ seedVRF := append(bytes.Repeat([]byte{0xff}, 7), make([]byte, 25)...)
+ f.Add(uint64(1954), uint64(1_999_999_999_999_964), uint64(1500), seedVRF)
+ f.Add(uint64(1141), uint64(1000), uint64(250), seedVRF)
+ f.Add(uint64(0), uint64(2_000_000_000_000_000), uint64(20), make([]byte, 32))
+ f.Add(uint64(1000), uint64(1000), uint64(1000), make([]byte, 32))
+ f.Add(uint64(100), uint64(1000), uint64(2000), make([]byte, 32)) // expectedSize > totalMoney
+ // all-0xff digest: ratio is exactly 1.0, the cdf-reaches-1.0 regime
+ f.Add(uint64(1954), uint64(1_999_999_999_999_964), uint64(1500), bytes.Repeat([]byte{0xff}, 32))
+ // fall-through to money by full walk (p=1/2: the pmf never drops below
+ // cum's half-ulp, so the CDF never freezes) and by the freeze short-circuit
+ // (tiny p: pmf underflows within a few steps); the second must equal the
+ // oracle's unshortened walk
+ f.Add(uint64(100), uint64(200), uint64(100), bytes.Repeat([]byte{0xff}, 32))
+ f.Add(uint64(1954), uint64(1_999_999_999_999_960), uint64(1500), bytes.Repeat([]byte{0xff}, 32))
+ // expectedSize > totalMoney with a nonzero digest: the degenerate path's money return
+ f.Add(uint64(100), uint64(1000), uint64(2000), bytes.Repeat([]byte{0xff}, 32))
+ // digest exactly halfway between ratio ulps: the denominator-correction
+ // round-up branch in f128FromDigestRatio, ~2^-128 density under mutation
+ halfway := make([]byte, DigestSize)
+ halfway[0], halfway[16] = 0x80, 0x80
+ f.Add(uint64(2000), uint64(4000), uint64(2000), halfway)
+ // tiny digest: the low-word normalization branch of the ratio conversion
+ tiny := make([]byte, DigestSize)
+ tiny[DigestSize-1] = 1
+ f.Add(uint64(1500), uint64(3000), uint64(1500), tiny)
+ mid := make([]byte, DigestSize)
+ mid[13] = 0x40 // leading zeros into the digest's third word
+ f.Add(uint64(1500), uint64(3000), uint64(1500), mid)
+ mid2 := make([]byte, DigestSize)
+ mid2[20] = 0x10 // leading zeros into the digest's second word
+ f.Add(uint64(1500), uint64(3000), uint64(1500), mid2)
+ // p just below 1: (1-p)^money exercises deep exponents; and extremes of
+ // total/expected magnitude the mutator will not reach from mid-range seeds
+ f.Add(uint64(2500), uint64(10_000_000_000_000_000), uint64(9_999_999_999_999_999), make([]byte, 32))
+ f.Add(uint64(3000), ^uint64(0), uint64(1), make([]byte, 32))
+ f.Add(uint64(7), uint64(0), uint64(0), make([]byte, 32)) // totalMoney == 0: degenerate nil path, zero ratio
+
+ // The uint64 expectedSize needs no input filtering: NaN/Inf/negative/
+ // fractional sizes are unrepresentable, and expectedSize >= totalMoney
+ // (including totalMoney == 0) takes the exact-integer p >= 1 path in both
+ // implementations.
+ f.Fuzz(func(t *testing.T, money, total, expected uint64, vrf []byte) {
+ money %= 3001 // bound the walk so each fuzz exec stays fast
+ var d Digest
+ copy(d[:], vrf)
+ got := SelectF128(money, total, expected, d)
+ want := selectBigOracle(money, total, expected, d)
+ if got != want {
+ t.Fatalf("SelectF128=%d != big.Float oracle=%d (money=%d total=%d expected=%d vrf=%x)",
+ got, want, money, total, expected, d)
+ }
+ })
+}
+
+// TestF128PrimitiveEdges pins the defensive arms of the shift and comparison
+// primitives that no production caller reaches (norm128 only shifts by
+// 1..127, and CDF boundaries are never zero, so neither the walk nor the
+// differential harnesses can cover them). They are total functions with
+// defined answers; assert them directly. The one remaining uncoverable branch
+// is divStep's add-back, which is unreachable under any inputs (see the
+// comment there).
+func TestF128PrimitiveEdges(t *testing.T) {
+ if hi, lo := shl128(5, 7, 0); hi != 5 || lo != 7 {
+ t.Fatalf("shl128 by 0: got %d,%d", hi, lo)
+ }
+ if hi, lo := shl128(5, 7, 128); hi != 0 || lo != 0 {
+ t.Fatalf("shl128 by 128: got %d,%d", hi, lo)
+ }
+ if hi, lo := shr128(5, 7, 200); hi != 0 || lo != 0 {
+ t.Fatalf("shr128 by 200: got %d,%d", hi, lo)
+ }
+ if !f128FromUint64(0).isZero() {
+ t.Fatal("f128FromUint64(0) is not zero")
+ }
+ one := f128FromUint64(1)
+ if c := (f128{}).cmp(f128{}); c != 0 {
+ t.Fatalf("cmp(0,0) = %d, want 0", c)
+ }
+ if c := one.cmp(f128{}); c != 1 {
+ t.Fatalf("cmp(1,0) = %d, want 1", c)
+ }
+ if c := (f128{}).cmp(one); c != -1 {
+ t.Fatalf("cmp(0,1) = %d, want -1", c)
+ }
+ if c := one.cmp(one); c != 0 {
+ t.Fatalf("cmp(1,1) = %d, want 0", c)
+ }
+}
+
+// TestF128AgreesWithCurrent checks broad agreement with the deployed
+// Boost-double implementation. Knife-edge differences remain expected because
+// SelectF128 uses an f128 digest ratio and CDF.
+func TestF128AgreesWithCurrent(t *testing.T) {
+ rng := rand.New(rand.NewSource(1))
+ const total = uint64(2_000_000_000_000_000)
+ committees := []uint64{20, 1500, 2990, 6000}
+ const n = 200000
+ match, considered, diffs := 0, 0, 0
+ for i := 0; i < n; i++ {
+ exp := committees[rng.Intn(len(committees))]
+ mean := 0.05 + rng.Float64()*30
+ money := uint64(mean * float64(total) / float64(exp))
+ if money == 0 {
+ continue
+ }
+ var d Digest
+ rng.Read(d[:])
+ considered++
+ cpp := Select(money, total, float64(exp), d)
+ f := SelectF128(money, total, exp, d)
+ if cpp == f {
+ match++
+ continue
+ }
+ diffs++
+ if diffs <= 10 {
+ t.Logf("knife-edge diff: money=%d exp=%d cpp=%d f128=%d", money, exp, cpp, f)
+ }
+ }
+ t.Logf("SelectF128 vs C++ Select: %d/%d agree (%d differ)", match, considered, diffs)
+ if match*1000 < considered*999 {
+ t.Errorf("agreement %d/%d too low", match, considered)
+ }
+}
+
+func TestF128DigestRatioMatchesBigFloat(t *testing.T) {
+ cases := make([]Digest, 0, 1389)
+ cases = append(cases, Digest{})
+ setBit := func(d *Digest, bit int) {
+ d[len(d)-1-bit/8] |= byte(1) << uint(bit%8)
+ }
+
+ // Exercise every possible normalization shift.
+ for bit := 0; bit < DigestSize*8; bit++ {
+ var d Digest
+ setBit(&d, bit)
+ cases = append(cases, d)
+ }
+
+ // Exercise exact halfway tails at every shift where a tail remains. The
+ // denominator correction must make each of these round upward.
+ for leading := 0; leading < f128MantBits; leading++ {
+ var d Digest
+ setBit(&d, DigestSize*8-1-leading)
+ setBit(&d, f128MantBits-1-leading)
+ cases = append(cases, d)
+ }
+
+ // Halfway tails whose sticky comes only from the lowest word: rounding up
+ // must come from the digest's own low bits, not the denominator
+ // correction.
+ for leading := 0; leading < 63; leading++ {
+ var d Digest
+ setBit(&d, DigestSize*8-1-leading)
+ setBit(&d, f128MantBits-1-leading)
+ setBit(&d, 0)
+ cases = append(cases, d)
+ }
+
+ var one Digest
+ one[len(one)-1] = 1
+ cases = append(cases, one)
+
+ var halfway Digest
+ halfway[0] = 0x80
+ halfway[16] = 0x80
+ cases = append(cases, halfway)
+
+ var maximum Digest
+ for i := range maximum {
+ maximum[i] = 0xff
+ }
+ cases = append(cases, maximum)
+
+ var nearMaximum Digest
+ for i := 0; i < 7; i++ {
+ nearMaximum[i] = 0xff
+ }
+ cases = append(cases, nearMaximum)
+
+ rng := rand.New(rand.NewSource(3))
+ for i := 0; i < 1000; i++ {
+ var d Digest
+ rng.Read(d[:])
+ cases = append(cases, d)
+ }
+
+ for _, d := range cases {
+ got := f128ToBig(f128FromDigestRatio(d))
+ want := digestRatioBig(d, f128MantBits)
+ if got.Cmp(want) != 0 {
+ t.Fatalf("digest ratio mismatch for %x: f128=%v big.Float=%v", d, got, want)
+ }
+ }
+}
+
+func TestSelectF128NearMaximumDigest(t *testing.T) {
+ var d Digest
+ for i := 0; i < 7; i++ {
+ d[i] = 0xff
+ }
+ got := SelectF128(1954, 1_999_999_999_999_964, 1500, d)
+ if got != 1 {
+ t.Fatalf("SelectF128=%d, want 1 for exact near-maximum digest ratio", got)
+ }
+}
+
+// TestSelectF128RatioExactlyOne pins the walk when the f128 ratio is exactly
+// 1.0: mathematically for the all-0xff digest, and by 128-bit rounding for any
+// digest with at least 129 leading one bits. With the f128-rounded threshold
+// fixed at 1.0, money is the exact-CDF count; the walk returns an earlier j
+// only when the accumulated f128 CDF happens to round up to exactly 1.0 (see
+// the SelectF128 doc comment). Each case pins one branch:
+//
+// - money=1954 with total=1_999_999_999_999_964: stops early at j=3, while
+// the same distribution with total=2_000_000_000_000_000 (36 more) falls
+// through to money. A hair-trigger pair pinned together: if a rounding
+// change flips either, the cdf trajectory moved by an ulp -- the walk did
+// not break.
+// - money=100, p=1/2: provably falls through to money -- cdf(99) is
+// 1 - 2^-100, which sits 2^28 ulps below 1.0, a gap no rounding can
+// bridge.
+// - money=129, p=1/2: the exact boundary -- true cdf(128) = 1 - 2^-129 is
+// precisely the rounding midpoint, and ties-to-even rounds it up to
+// exactly 1.0, stopping at j=128.
+func TestSelectF128RatioExactlyOne(t *testing.T) {
+ one := f128FromUint64(1)
+
+ var maximum Digest
+ for i := range maximum {
+ maximum[i] = 0xff
+ }
+ // exactly 129 leading one bits: the minimal digest that rounds to 1.0
+ var minLeadingOnes Digest
+ for i := 0; i < 16; i++ {
+ minLeadingOnes[i] = 0xff
+ }
+ minLeadingOnes[16] = 0x80
+
+ cases := []struct {
+ money, total, expected uint64
+ want uint64
+ }{
+ {1954, 1_999_999_999_999_964, 1500, 3},
+ {1954, 2_000_000_000_000_000, 1500, 1954},
+ {100, 200, 100, 100},
+ {129, 258, 129, 128},
+ }
+ for _, d := range []Digest{maximum, minLeadingOnes} {
+ if f128FromDigestRatio(d).cmp(one) != 0 {
+ t.Fatalf("digest %x: ratio is not exactly 1.0", d)
+ }
+ for _, c := range cases {
+ got := SelectF128(c.money, c.total, c.expected, d)
+ if oracle := selectBigOracle(c.money, c.total, c.expected, d); got != oracle {
+ t.Fatalf("digest %x money=%d: SelectF128=%d != oracle=%d", d, c.money, got, oracle)
+ }
+ if got != c.want {
+ t.Fatalf("digest %x money=%d: SelectF128=%d, want %d", d, c.money, got, c.want)
+ }
+ }
+ }
+}
+
+// FuzzF128Ops validates the f128 arithmetic primitives directly against
+// math/big.Float at the same mantissa width.
+func FuzzF128Ops(f *testing.F) {
+ f.Add(uint64(1)<<63, uint64(0), 0, uint64(3)<<62, uint64(0), 0, uint64(7))
+ f.Add(uint64(0), uint64(0), 0, uint64(1)<<63, uint64(1), -5, uint64(1))
+ // Large divisors: go's mutator only explores integers near corpus values,
+ // so without these seeds ~30M execs never left the small-u neighborhood
+ // and missed divU's shallow-quotient misrounding (u > ~2^62).
+ f.Add(uint64(1)<<63|12345, uint64(67890), 3, uint64(1)<<63, uint64(1), 0, uint64(1)<<63|54321)
+ f.Add(uint64(1)<<63|999, uint64(777), -9, uint64(1)<<63, uint64(1), 0, uint64(5)<<61|33)
+ // For the same reason, each regime below gets its own seed. Exponents are
+ // encoded as E+2000 (the body maps aexp%4000-2000). Exact ties, carry-out
+ // renormalization, and the rare divStep branches have ~2^-64 density under
+ // uniform inputs; these vectors were constructed or mined by instrumented
+ // search and verified against big.Float:
+ f.Add(uint64(1)<<63, uint64(1), 2000, uint64(3)<<62, uint64(0), 2000, uint64(2)) // mul tail exactly half, odd mantissa: tie rounds up
+ f.Add(uint64(1)<<63, uint64(3), 2000, uint64(3)<<62, uint64(0), 2000, uint64(2)) // mul tie, even mantissa: ties to even
+ f.Add(^uint64(0), ^uint64(0), 2128, uint64(1)<<63, uint64(0), 2000, uint64(3)) // add tie at exp gap 128 into all-ones: carry renormalizes
+ f.Add(uint64(1)<<63, uint64(1), 2129, uint64(1)<<63, uint64(0), 2000, uint64(3)) // add exp gap 129: addend entirely below the round bit
+ f.Add(uint64(1)<<63, uint64(5), 2000, uint64(1)<<63, uint64(100), 2000, uint64(6)) // div: step remainder high word reaches v1 (qhat cap)
+ f.Add(uint64(1)<<63, uint64(1)<<63, 2000, uint64(1)<<63, uint64(1)<<63|2, 2000, uint64(6)) // div: qhat cap where rhat carries out
+ f.Add(uint64(0xff2432b605ae124e), uint64(0x7e0bcdcc481c2dbd), 2000,
+ uint64(0xa7c5c2e99ad4c9a0), uint64(0xd065b805fe1d2cf5), 2000, uint64(9)) // div: refine decrement hits the rhat overflow break
+ f.Add(uint64(0), uint64(12345), 2000, uint64(1), uint64(0), 2000, uint64(2)) // denormalized mantissas: norm128/shl128 paths
+ f.Add(uint64(1)<<63, uint64(0), 2000, uint64(0), uint64(0), 2000, uint64(5)) // zero second operand
+ // sticky carried ONLY by the term a mutation could drop (found by the
+ // mutation campaign in mutation_check.go): a mul tie broken only by p0,
+ // an add tie at gap 128 broken only by the addend's low word, and a divU
+ // tie broken only by the division remainder
+ f.Add(uint64(3)<<62|1, uint64(1), 2000, uint64(3)<<62-1, uint64(1), 2000, uint64(11))
+ f.Add(uint64(1)<<63, uint64(2), 2128, uint64(1)<<63, uint64(5), 2000, uint64(3))
+ f.Add(uint64(1)<<63, uint64(0), 2000, uint64(1)<<63, uint64(1), 2000, ^uint64(0))
+ f.Fuzz(func(t *testing.T, ahi, alo uint64, aexp int, bhi, blo uint64, bexp int, u uint64) {
+ a := norm128(ahi, alo, int64(aexp%4000-2000))
+ b := norm128(bhi, blo, int64(bexp%4000-2000))
+ ab, bb := f128ToBig(a), f128ToBig(b)
+ // norm128 is this harness's own input constructor, so nothing
+ // downstream would notice it corrupting the value (a mutation
+ // campaign caught exactly that); check it against the raw words.
+ raw := new(big.Float).SetPrec(300).SetUint64(ahi)
+ raw.SetMantExp(raw, 64)
+ raw.Add(raw, new(big.Float).SetPrec(300).SetUint64(alo))
+ raw.SetMantExp(raw, aexp%4000-2000)
+ if raw.Cmp(ab) != 0 {
+ t.Fatalf("norm128 changed the value: raw=%v normalized=%v (ahi=%#x alo=%#x)", raw, ab, ahi, alo)
+ }
+ check := func(name string, got f128, want *big.Float) {
+ gotBig := f128ToBig(got)
+ if gotBig.Cmp(want) != 0 {
+ t.Fatalf("%s: f128=%v big.Float=%v (a=%v b=%v u=%d)", name, gotBig, want, ab, bb, u)
+ }
+ }
+ check("mul", a.mul(b), new(big.Float).SetPrec(f128MantBits).Mul(ab, bb))
+ check("add", a.add(b), new(big.Float).SetPrec(f128MantBits).Add(ab, bb))
+ if u != 0 {
+ check("divU", a.divU(u),
+ new(big.Float).SetPrec(f128MantBits).Quo(ab, new(big.Float).SetPrec(f128MantBits).SetUint64(u)))
+ }
+ if !b.isZero() {
+ check("div", a.div(b), new(big.Float).SetPrec(f128MantBits).Quo(ab, bb))
+ }
+ })
+}
+
+// TestF128MatchesOracleLargeMoney extends the strict SelectF128 == oracle check
+// beyond FuzzSelectF128's money<=3000 cap to realistic LARGE money (up to ~2^51),
+// where the big.Float oracle is still exact (mean = money*p stays a committee
+// size). Each committee's mean is drawn in [0.05, size] so money stays <= total.
+func TestF128MatchesOracleLargeMoney(t *testing.T) {
+ rng := rand.New(rand.NewSource(2))
+ const total = uint64(2_000_000_000_000_000)
+ committees := []uint64{20, 1500, 2990, 6000}
+ for i := 0; i < 3000; i++ {
+ size := committees[rng.Intn(len(committees))]
+ mean := 0.05 + rng.Float64()*float64(size) // mean <= size => money <= total
+ money := uint64(mean * float64(total) / float64(size)) // up to ~2^51
+ if money == 0 {
+ continue
+ }
+ var d Digest
+ rng.Read(d[:])
+ if got, want := SelectF128(money, total, size, d), selectBigOracle(money, total, size, d); got != want {
+ t.Fatalf("SelectF128=%d != oracle=%d (money=%d size=%d vrf=%x)", got, want, money, size, d)
+ }
+ }
+}
+
+// BenchmarkSelectF128 mirrors BenchmarkSortition (same parameters) so the pure-Go
+// deterministic path can be compared directly against the cgo/Boost Select.
+func BenchmarkSelectF128(b *testing.B) {
+ b.StopTimer()
+ keys := make([]Digest, b.N)
+ for i := 0; i < b.N; i++ {
+ rand.Read(keys[i][:])
+ }
+ b.StartTimer()
+ for i := 0; i < b.N; i++ {
+ SelectF128(1000000, 1000000000000, 2500, keys[i])
+ }
+}
+
+// f128ToBig returns the exact value of an f128 as a big.Float (test helper).
+func f128ToBig(x f128) *big.Float {
+ hi := new(big.Float).SetPrec(300).SetUint64(x.hi)
+ hi.SetMantExp(hi, 64)
+ m := new(big.Float).SetPrec(300).Add(hi, new(big.Float).SetPrec(300).SetUint64(x.lo))
+ return m.SetMantExp(m, int(x.exp))
+}
+
+// TestDivVsBig checks f128.div against a 128-bit round-nearest-even big.Float
+// divide on broad random operands. The SelectF128 fuzz only exercises div with
+// the p/(1-p) shape, so this independently validates the hand-rolled Knuth
+// long division across the full input space (spanning exponents and the
+// Q<2^128 vs Q>=2^128 normalization cases).
+func TestDivVsBig(t *testing.T) {
+ rng := rand.New(rand.NewSource(2))
+ randF128 := func() f128 {
+ return f128{rng.Uint64() | 1<<63, rng.Uint64(), int64(rng.Intn(4000) - 2000)} // normalized
+ }
+ for i := 0; i < 5_000_000; i++ {
+ a, b := randF128(), randF128()
+ got := f128ToBig(a.div(b))
+ want := new(big.Float).SetPrec(f128MantBits).Quo(f128ToBig(a), f128ToBig(b))
+ if got.Cmp(want) != 0 {
+ t.Fatalf("div mismatch: a=%+v b=%+v got=%s want=%s", a, b, got.Text('p', 0), want.Text('p', 0))
+ }
+ }
+}
+
+// maxDigestMinusPowerOfTwo returns the digest integer 2^256-1-2^bit.
+func maxDigestMinusPowerOfTwo(bit uint) Digest {
+ var d Digest
+ for i := range d {
+ d[i] = 0xff
+ }
+ d[len(d)-1-int(bit/8)] &^= byte(1) << (bit % 8)
+ return d
+}
+
+// TestDivUVsBig checks divU against a 128-bit round-nearest-even big.Float
+// divide across divisor magnitudes. The log-spread divisor matters: the
+// shallow-quotient region (u > ~2^62, where the quotient has at most 129
+// significant bits) is unreachable from the sortition walk, whose divisor is
+// the step index bounded by money, but the round-to-nearest-even contract
+// covers it, and a previous divU broke there by folding the remainder's
+// sticky marker into a digit that landed in or at the rounding position.
+func TestDivUVsBig(t *testing.T) {
+ rng := rand.New(rand.NewSource(5))
+ for i := 0; i < 2_000_000; i++ {
+ a := f128{rng.Uint64() | 1<<63, rng.Uint64(), int64(rng.Intn(4000) - 2000)}
+ u := rng.Uint64() >> uint(rng.Intn(64)) // log-spread magnitudes
+ if u == 0 {
+ continue
+ }
+ got := f128ToBig(a.divU(u))
+ want := new(big.Float).SetPrec(f128MantBits).Quo(
+ f128ToBig(a), new(big.Float).SetPrec(f128MantBits).SetUint64(u))
+ if got.Cmp(want) != 0 {
+ t.Fatalf("divU mismatch: a={%#x,%#x,%d} u=%d got=%s want=%s",
+ a.hi, a.lo, a.exp, u, got.Text('p', 0), want.Text('p', 0))
+ }
+ }
+}
+
+// TestSelectF128CurrentConsensusFrozenTail pins the accepted frozen-tail
+// behavior at values admitted by current go-algorand consensus parameters.
+// Consensus v41 inherits NumProposers=20, NextCommitteeSize=5000, and
+// MinBalance=100,000 microalgos. Its payout-eligibility interval is 30,000
+// through 70,000,000 Algos, and the mainnet genesis supply is 10,000,000,000
+// Algos. The payout maximum is not a voting-stake cap--online accounts above
+// it can still take part in consensus--so the cases cover a proposer plus the
+// base account minimum, both payout landmarks, and the supply ceiling.
+//
+// In every case q=(1-p) rounds downward. Raising q to money scales every PMF
+// term down enough that the accumulated f128 CDF freezes below the chosen
+// digest ratio. SelectF128 defines this interval to return money; completion
+// is also the liveness assertion, since the unshortened walk would perform up
+// to money no-op iterations. The deployed Boost walk does not share the
+// plateau: the digest rounds to binary64 1.0, and its independently evaluated
+// CDF reaches 1.0 at the finite values pinned in boostWant.
+func TestSelectF128CurrentConsensusFrozenTail(t *testing.T) {
+ const (
+ mainnetSupply = uint64(10_000_000_000_000_000)
+ )
+ tests := []struct {
+ name string
+ money uint64
+ total uint64
+ expected uint64
+ clearBit uint
+ boostWant uint64
+ }{
+ {"proposer committee", 1_999_999_999_999_999, 1_999_999_999_999_999, 20, 175, 67},
+ {"base minimum balance", 100_000, mainnetSupply, 5_000, 141, 2},
+ {"payout minimum balance", 30_000_000_000, mainnetSupply, 5_000, 159, 6},
+ {"payout maximum balance", 70_000_000_000_000, mainnetSupply, 5_000, 170, 94},
+ {"mainnet supply ceiling", mainnetSupply, mainnetSupply, 5_000, 178, 5_598},
+ }
+ for _, test := range tests {
+ t.Run(test.name, func(t *testing.T) {
+ d := maxDigestMinusPowerOfTwo(test.clearBit)
+ if got := SelectF128(test.money, test.total, test.expected, d); got != test.money {
+ t.Fatalf("SelectF128=%d, want money=%d for a ratio above the CDF plateau", got, test.money)
+ }
+ if got := Select(test.money, test.total, float64(test.expected), d); got != test.boostWant {
+ t.Fatalf("Boost Select=%d, want finite tail count %d", got, test.boostWant)
+ }
+ })
+ }
+}
+
+// TestSelectF128FrozenTailReportedCase retains the original supply-sized
+// reproducer: pmf(0)'s trial-count-amplified rounding leaves the accumulated
+// CDF around 2^-78 below 1 when money == totalMoney == 2e15 and committee size
+// is the current certification size 1500. The ratio 1-2^-80 is above that
+// plateau, while Boost terminates at its binary64 tail boundary.
+func TestSelectF128FrozenTailReportedCase(t *testing.T) {
+ const onlineStake = uint64(2_000_000_000_000_000)
+ d := maxDigestMinusPowerOfTwo(176) // ratio ~= 1 - 2^-80
+ if got := SelectF128(onlineStake, onlineStake, 1500, d); got != onlineStake {
+ t.Fatalf("SelectF128=%d, want money=%d for a ratio above the CDF plateau", got, onlineStake)
+ }
+ if got := Select(onlineStake, onlineStake, 1500, d); got != 1832 {
+ t.Fatalf("Boost Select=%d, want finite tail count 1832", got)
+ }
+
+ // The all-0xff digest (ratio exactly 1.0) is also above this
+ // distribution's plateau and must take the same frozen path.
+ for i := range d {
+ d[i] = 0xff
+ }
+ if got := SelectF128(onlineStake, onlineStake, 1500, d); got != onlineStake {
+ t.Fatalf("SelectF128=%d, want money=%d for ratio exactly 1.0", got, onlineStake)
+ }
+}
diff --git a/f128_trajectory_test.go b/f128_trajectory_test.go
new file mode 100644
index 0000000..019bf91
--- /dev/null
+++ b/f128_trajectory_test.go
@@ -0,0 +1,353 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "math/big"
+ "testing"
+)
+
+// auditedMul performs the same multiplication used by intPow while proving
+// with unbounded integers that both possible normalization exponents fit in
+// int64. Checking both +127 and +128 is conservative and avoids duplicating
+// mul's partial-product branch decision in the test.
+func auditedMul(t *testing.T, a, b f128, minExp, maxExp *int64) f128 {
+ t.Helper()
+ if !a.isZero() && !b.isZero() {
+ for _, normalization := range []int64{127, 128} {
+ e := new(big.Int).SetInt64(a.exp)
+ e.Add(e, new(big.Int).SetInt64(b.exp))
+ e.Add(e, new(big.Int).SetInt64(normalization))
+ if !e.IsInt64() {
+ t.Fatalf("mul exponent overflows int64: %d + %d + %d = %s", a.exp, b.exp, normalization, e)
+ }
+ }
+ }
+ out := a.mul(b)
+ assertCanonicalF128(t, "audited mul", out)
+ if !out.isZero() {
+ if out.exp < *minExp {
+ *minExp = out.exp
+ }
+ if out.exp > *maxExp {
+ *maxExp = out.exp
+ }
+ }
+ return out
+}
+
+func auditedIntPow(t *testing.T, base f128, e uint64) (f128, int64, int64) {
+ t.Helper()
+ minExp, maxExp := base.exp, base.exp
+ result, b := f128FromUint64(1), base
+ for e > 0 {
+ if e&1 == 1 {
+ result = auditedMul(t, result, b, &minExp, &maxExp)
+ }
+ e >>= 1
+ if e > 0 {
+ b = auditedMul(t, b, b, &minExp, &maxExp)
+ }
+ }
+ return result, minExp, maxExp
+}
+
+// TestSelectF128MaxDomainExponentAudit exercises the documented worst case:
+// money is one below the exported limit and 1-p is the smallest positive
+// rational representable by the uint64 API. The test audits every exponent
+// addition in exponentiation-by-squaring with big.Int, so a future bound or
+// normalization change fails before platform-independent int64 arithmetic can
+// wrap.
+func TestSelectF128MaxDomainExponentAudit(t *testing.T) {
+ const total = ^uint64(0)
+ money := SelectF128MaxMoney - 1
+ expected := total - 1 // 1-p = 1/(2^64-1)
+ q := f128FromUint64(1).div(f128FromUint64(total))
+ pmf, minExp, maxExp := auditedIntPow(t, q, money)
+ dist := newBinomialF128(expected, total, money)
+ if dist == nil || pmf != dist.pmf {
+ t.Fatalf("audited pmf(0)=%+v, constructor produced %+v", pmf, dist)
+ }
+
+ // q is just above 2^-64. Across money < 2^56 its logarithmic correction
+ // is less than one, so floor(log2(q^money)) is exactly -64*money.
+ wantExp := -int64(money<<6) - 127
+ if pmf.exp != wantExp {
+ t.Fatalf("domain-edge pmf exponent=%d, want %d", pmf.exp, wantExp)
+ }
+ if minExp < wantExp || maxExp > 0 {
+ t.Fatalf("unexpected audited exponent range [%d,%d], final lower bound %d", minExp, maxExp, wantExp)
+ }
+ if got := SelectF128(money, total, expected, Digest{}); got != 0 {
+ t.Fatalf("domain-edge ratio zero selected %d, want 0", got)
+ }
+}
+
+// TestSelectF128DegenerateDistributionShortCircuit pins construction, not
+// only the eventual return value. If p==1 accidentally enters the ordinary
+// walk, its all-zero PMF produces the right answer for nonzero digests only
+// after money iterations, recreating the liveness class that freeze prevents.
+func TestSelectF128DegenerateDistributionShortCircuit(t *testing.T) {
+ const money = SelectF128MaxMoney - 1
+ for _, total := range []uint64{0, 1, 6000, ^uint64(0)} {
+ if dist := newBinomialF128(total, total, money); dist != nil {
+ t.Fatalf("p==1 constructed an ordinary distribution (total=%d): %+v", total, dist)
+ }
+ if total != ^uint64(0) {
+ if dist := newBinomialF128(total+1, total, money); dist != nil {
+ t.Fatalf("p>1 constructed an ordinary distribution (total=%d): %+v", total, dist)
+ }
+ }
+ }
+}
+
+type highPrecisionBinomial struct {
+ money uint64
+ prec uint
+ at uint64
+ pq *big.Float
+ pmf *big.Float
+ cum *big.Float
+}
+
+func newHighPrecisionBinomial(money, total, expected uint64, prec uint) *highPrecisionBinomial {
+ q := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(total-expected),
+ new(big.Float).SetPrec(prec).SetUint64(total))
+ pq := new(big.Float).SetPrec(prec).Quo(
+ new(big.Float).SetPrec(prec).SetUint64(expected),
+ new(big.Float).SetPrec(prec).SetUint64(total-expected))
+ pmf := bigIntPow(q, money, prec)
+ return &highPrecisionBinomial{
+ money: money,
+ prec: prec,
+ pq: pq,
+ pmf: pmf,
+ cum: new(big.Float).SetPrec(prec).Set(pmf),
+ }
+}
+
+func (b *highPrecisionBinomial) advance() {
+ b.at++
+ factor := new(big.Float).SetPrec(b.prec).Quo(
+ new(big.Float).SetPrec(b.prec).SetUint64(b.money-b.at+1),
+ new(big.Float).SetPrec(b.prec).SetUint64(b.at))
+ step := new(big.Float).SetPrec(b.prec).Mul(factor, b.pq)
+ b.pmf = new(big.Float).SetPrec(b.prec).Mul(b.pmf, step)
+ b.cum = new(big.Float).SetPrec(b.prec).Add(b.cum, b.pmf)
+}
+
+func binomialMode(money, total, expected uint64) uint64 {
+ n := new(big.Int).SetUint64(money)
+ n.Add(n, big.NewInt(1))
+ n.Mul(n, new(big.Int).SetUint64(expected))
+ n.Quo(n, new(big.Int).SetUint64(total))
+ return n.Uint64()
+}
+
+func f128AsBig(x f128, prec uint) *big.Float {
+ return new(big.Float).SetPrec(prec).Set(f128ToBig(x))
+}
+
+func absFloatDifference(a, b *big.Float, prec uint) *big.Float {
+ d := new(big.Float).SetPrec(prec).Sub(a, b)
+ return d.Abs(d)
+}
+
+func trajectoryErrorBound(money, at uint64, prec uint) *big.Float {
+ // RNE unit roundoff is 2^-128. Sixty-four ulps per initialization trial
+ // and recurrence index is deliberately conservative: it covers q^money's
+ // amplified input rounding, four rounded recurrence operations per term,
+ // and CDF summation without turning this into an empirical pinned maximum.
+ units := new(big.Int).SetUint64(money)
+ units.Add(units, new(big.Int).SetUint64(32*(at+1)+128))
+ return new(big.Float).SetMantExp(new(big.Float).SetPrec(prec).SetInt(units), -122)
+}
+
+func assertTrajectoryError(t *testing.T, name string, at uint64, got f128, reference, bound *big.Float, relative bool) {
+ t.Helper()
+ const prec = uint(512)
+ diff := absFloatDifference(f128AsBig(got, prec), reference, prec)
+ if relative {
+ if reference.Sign() == 0 {
+ if diff.Sign() != 0 {
+ t.Fatalf("%s(%d): nonzero value against zero reference", name, at)
+ }
+ return
+ }
+ diff.Quo(diff, reference)
+ }
+ if diff.Cmp(bound) > 0 {
+ t.Fatalf("%s(%d) error %s exceeds bound %s", name, at, diff.Text('p', 8), bound.Text('p', 8))
+ }
+}
+
+// TestSelectF128TrajectoryErrorBudget compares every production PMF/CDF state
+// through and beyond the mode with a 512-bit recurrence. Exact/Arb tests cover
+// shared-formula risk; this test instead detects substantial internal drift
+// that happens not to move the sampled final quantile.
+func TestSelectF128TrajectoryErrorBudget(t *testing.T) {
+ tests := []struct {
+ name string
+ money, total, expected uint64
+ maxAt uint64
+ }{
+ {"toy balanced", 256, 512, 256, 255},
+ {"online 1500", 2_000_000_000_000_000, 2_000_000_000_000_000, 1500, 2500},
+ {"one third supply", 10_000_000_000_000_000 / 3, 10_000_000_000_000_000, 2990, 1800},
+ {"payout minimum", 30_000_000_000, 10_000_000_000_000_000, 5000, 120},
+ }
+ for _, test := range tests {
+ t.Run(test.name, func(t *testing.T) {
+ const prec = uint(512)
+ got := newBinomialF128(test.expected, test.total, test.money)
+ want := newHighPrecisionBinomial(test.money, test.total, test.expected, prec)
+ mode := binomialMode(test.money, test.total, test.expected)
+ observed := uint64(0)
+ for at := uint64(0); at <= test.maxAt; at++ {
+ if at > 0 {
+ pmfPrev, cumPrev := got.pmf, got.cum
+ got.cdf(at)
+ want.advance()
+ if got.cum.cmp(cumPrev) < 0 {
+ t.Fatalf("CDF decreased at %d", at)
+ }
+ if at <= mode && got.pmf.cmp(pmfPrev) < 0 {
+ t.Fatalf("PMF decreased before mode %d at %d", mode, at)
+ }
+ if at > mode && got.pmf.cmp(pmfPrev) > 0 {
+ t.Fatalf("PMF increased after mode %d at %d", mode, at)
+ }
+ }
+ bound := trajectoryErrorBound(test.money, at, prec)
+ assertTrajectoryError(t, "pmf", at, got.pmf, want.pmf, bound, true)
+ assertTrajectoryError(t, "cdf", at, got.cum, want.cum, bound, false)
+ observed++
+ if got.frozen {
+ break
+ }
+ }
+ if observed <= mode {
+ t.Fatalf("trajectory stopped after %d states before crossing mode %d", observed, mode)
+ }
+ })
+ }
+}
+
+// TestSelectF128FreezePermanence ignores the production short-circuit after a
+// bounded walk freezes and explicitly advances the recurrence to money-1.
+// Every later PMF must remain non-increasing and every rounded CDF add must be
+// a no-op, proving on the exercised grid that the optimization returns the
+// same answer as the otherwise impractical unshortened walk.
+func TestSelectF128FreezePermanence(t *testing.T) {
+ tests := []struct {
+ money, total, expected uint64
+ }{
+ {256, 512, 256},
+ {512, 1000, 1},
+ {1000, 10_000, 20},
+ {3000, 2_000_000_000_000_000, 1500},
+ }
+ frozenCases, checkedTailSteps := 0, uint64(0)
+ for _, test := range tests {
+ b := newBinomialF128(test.expected, test.total, test.money)
+ for j := uint64(1); j < test.money && !b.frozen; j++ {
+ b.cdf(j)
+ }
+ if !b.frozen {
+ continue
+ }
+ frozenCases++
+ pmf, cum, at := b.pmf, b.cum, b.at
+ for at+1 < test.money {
+ at++
+ pmfNext := f128FromUint64(test.money - at + 1).divU(at).mul(b.pq).mul(pmf)
+ if pmfNext.cmp(pmf) > 0 {
+ t.Fatalf("PMF increased after freeze at %d (money=%d total=%d expected=%d)", at, test.money, test.total, test.expected)
+ }
+ if cumNext := cum.add(pmfNext); cumNext != cum {
+ t.Fatalf("CDF changed after freeze at %d (money=%d total=%d expected=%d)", at, test.money, test.total, test.expected)
+ }
+ pmf = pmfNext
+ checkedTailSteps++
+ }
+ }
+ if frozenCases < 3 || checkedTailSteps < 100 {
+ t.Fatalf("anti-vacuity: frozenCases=%d checkedTailSteps=%d", frozenCases, checkedTailSteps)
+ }
+}
+
+func selectF128WithStepCount(money, total, expected uint64, d Digest) (selected, evaluations uint64, froze bool) {
+ ratio := f128FromDigestRatio(d)
+ dist := newBinomialF128(expected, total, money)
+ if dist == nil {
+ if ratio.isZero() {
+ return 0, 0, false
+ }
+ return money, 0, false
+ }
+ for j := uint64(0); j < money; j++ {
+ evaluations++
+ boundary := dist.cdf(j)
+ if ratio.cmp(boundary) <= 0 {
+ return j, evaluations, false
+ }
+ if dist.frozen {
+ return money, evaluations, true
+ }
+ }
+ return money, evaluations, false
+}
+
+// TestSelectF128ConsensusStepBounds makes liveness deterministic by counting
+// CDF evaluations rather than timing them. The cases cover certified tails
+// outside the frozen sliver and defined money-returning cases inside it.
+func TestSelectF128ConsensusStepBounds(t *testing.T) {
+ const (
+ online = uint64(2_000_000_000_000_000)
+ supply = uint64(10_000_000_000_000_000)
+ )
+ tests := []struct {
+ name string
+ money, total, expected uint64
+ digest Digest
+ want uint64
+ wantFreeze bool
+ }{
+ {"online 1500 certified tail", online, online, 1500, maxDigestMinusPowerOfTwo(196), 1852, false},
+ {"online 6000 certified tail", online, online, 6000, maxDigestMinusPowerOfTwo(200), 6667, false},
+ {"supply 5000 certified tail", supply, supply, 5000, maxDigestMinusPowerOfTwo(190), 5667, false},
+ {"proposer frozen tail", online - 1, online - 1, 20, maxDigestMinusPowerOfTwo(175), online - 1, true},
+ {"base minimum frozen tail", 100_000, supply, 5000, maxDigestMinusPowerOfTwo(141), 100_000, true},
+ {"supply frozen tail", supply, supply, 5000, maxDigestMinusPowerOfTwo(178), supply, true},
+ }
+ for _, test := range tests {
+ t.Run(test.name, func(t *testing.T) {
+ got, evaluations, froze := selectF128WithStepCount(test.money, test.total, test.expected, test.digest)
+ if got != test.want || got != SelectF128(test.money, test.total, test.expected, test.digest) {
+ t.Fatalf("instrumented selection=%d, want %d", got, test.want)
+ }
+ if froze != test.wantFreeze {
+ t.Fatalf("freeze=%v, want %v after %d evaluations", froze, test.wantFreeze, evaluations)
+ }
+ limit := test.expected + 4096
+ if evaluations == 0 || evaluations > limit {
+ t.Fatalf("walk used %d CDF evaluations, limit %d", evaluations, limit)
+ }
+ })
+ }
+}
diff --git a/go.mod b/go.mod
index 96fec71..ad11e8a 100644
--- a/go.mod
+++ b/go.mod
@@ -1,3 +1,5 @@
module github.com/algorand/sortition
-go 1.17
+go 1.23
+
+require pgregory.net/rapid v1.3.0
diff --git a/go.sum b/go.sum
new file mode 100644
index 0000000..2ecce80
--- /dev/null
+++ b/go.sum
@@ -0,0 +1,2 @@
+pgregory.net/rapid v1.3.0 h1:vBvO0VSqti75J1jjYqpgPNBLKMd1+gxa9fYo7vk/Exc=
+pgregory.net/rapid v1.3.0/go.mod h1:dPlE4OBBxgXPqkP79flB6sJL1dx5azpI7HQ9MY9Z7uk=
diff --git a/mutation_check.go b/mutation_check.go
new file mode 100644
index 0000000..4c9f22c
--- /dev/null
+++ b/mutation_check.go
@@ -0,0 +1,253 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+//go:build ignore
+
+// mutation_check applies curated single-site mutants to f128.go and its test
+// oracles, then verifies the fast suite kills each one. Every mutant models a
+// plausible implementation or harness bug (flipped rounding masks, dropped
+// carries, wrong initialization, corrupt recurrence/boundary logic); an
+// UNEXPECTED survivor is a test gap and must be answered with a new test.
+// Mutants marked equivalent encode a proof of why the mutation cannot change
+// behavior; they are expected to survive and document that analysis.
+//
+// Run from the repo root:
+//
+// go run mutation_check.go
+package main
+
+import (
+ "fmt"
+ "os"
+ "os/exec"
+ "strings"
+)
+
+type mutant struct {
+ name string
+ target string // defaults to f128.go; test-oracle mutants name their file
+ old, new string
+ equivalent string // non-empty: why this mutant is expected to survive
+}
+
+var mutants = []mutant{
+ {name: "mul-roundbit-hi", old: "roundBit := p1&(1<<63) != 0", new: "roundBit := p1&(1<<62) != 0"},
+ {name: "mul-sticky-hi-drop-p0", old: "sticky := (p1&^(uint64(1)<<63) != 0) || p0 != 0", new: "sticky := p1&^(uint64(1)<<63) != 0"},
+ {name: "mul-carry-or", old: "cp1 := cA + cB", new: "cp1 := cA | cB"},
+ {name: "mul-exp-hi", old: "roundNE(p3, p2, roundBit, sticky, a.exp+b.exp+128)", new: "roundNE(p3, p2, roundBit, sticky, a.exp+b.exp+129)"},
+ {name: "mul-roundbit-lo", old: "roundBit := p1&(1<<62) != 0", new: "roundBit := p1&(1<<61) != 0"},
+ {
+ name: "mul-sticky-lo-mask", old: "sticky := (p1&^(uint64(3)<<62) != 0) || p0 != 0", new: "sticky := (p1&^(uint64(1)<<62) != 0) || p0 != 0",
+ equivalent: "p1 bit 63 is also the mantissa's low bit in this branch, so whenever the leaked sticky could matter (round bit set) the mantissa is odd and ties round up anyway; clearing bit 63 is necessary only cosmetically",
+ },
+ {name: "roundNE-invert-tie", old: "if roundBit && (sticky || lo&1 != 0) {", new: "if roundBit && (sticky || lo&1 == 0) {"},
+ {name: "roundNE-carry-exp", old: "return f128{1 << 63, 0, exp + 1}", new: "return f128{1 << 63, 0, exp}"},
+ {name: "gs-roundbit-pos", old: "round = (lo>>(n-1))&1 != 0", new: "round = (lo>>n)&1 != 0"},
+ {name: "gs-sticky-guard", old: "if n >= 2 {", new: "if n >= 3 {"},
+ {name: "gs-128-drop-lo", old: "sticky = (hi&^(uint64(1)<<63) != 0) || lo != 0", new: "sticky = hi&^(uint64(1)<<63) != 0"},
+ {name: "add-gap-boundary", old: "if a.exp-b.exp > 128 {", new: "if a.exp-b.exp > 127 {"},
+ {name: "add-carry-drop-sticky", old: "sticky = sticky || round", new: "sticky = round"},
+ {name: "add-carry-drop-round", old: "round = slo&1 != 0", new: "round = false"},
+ {name: "divU-drop-rem-sticky", old: "sticky := o1&^(uint64(1)<<63) != 0 || o0 != 0 || rem != 0", new: "sticky := o1&^(uint64(1)<<63) != 0 || o0 != 0 || rem > rem"},
+ {name: "divU-roundbit", old: "round := o1&(uint64(1)<<63) != 0", new: "round := o1&(uint64(1)<<62) != 0"},
+ {name: "divU-exp", old: "roundNE(o3, o2, round, sticky, a.exp-int64(lz))", new: "roundNE(o3, o2, round, sticky, a.exp-int64(lz)-1)"},
+ {name: "ratio-drop-halfway", old: "if roundBit && !sticky {", new: "if roundBit && sticky {"},
+ {name: "ratio-exp", old: "roundNE(n3, n2, roundBit, sticky, -128-int64(leading))", new: "roundNE(n3, n2, roundBit, sticky, -129-int64(leading))"},
+ {
+ name: "ratio-sticky-drop-n0", old: "sticky := n1&^(uint64(1)<<63) != 0 || n0 != 0", new: "sticky := n1&^(uint64(1)<<63) != 0 || n0 > n0",
+ equivalent: "the halfway denominator-correction turns any roundBit-without-sticky into sticky, so n0's contribution is recreated exactly when it could matter, and sticky is irrelevant when roundBit is clear",
+ },
+ {name: "digest-little-endian", old: "w3 := binary.BigEndian.Uint64(d[0:8])", new: "w3 := binary.LittleEndian.Uint64(d[0:8])"},
+ {name: "intpow-invert-bit", old: "if e&1 == 1 {", new: "if e&1 == 0 {"},
+ {name: "intpow-wrong-square", old: "b = b.mul(b)", new: "b = b.mul(base)"},
+ {name: "distribution-degenerate-strict", old: "if expectedSize >= totalMoney { // p >= 1", new: "if expectedSize > totalMoney { // p >= 1"},
+ {
+ name: "distribution-q-numerator",
+ old: "qf := f128FromUint64(totalMoney - expectedSize).div(f128FromUint64(totalMoney)) // 1-p",
+ new: "qf := f128FromUint64(expectedSize).div(f128FromUint64(totalMoney)) // 1-p",
+ },
+ {
+ name: "distribution-pq-denominator",
+ old: "pq := f128FromUint64(expectedSize).div(f128FromUint64(totalMoney - expectedSize)) // p/(1-p)",
+ new: "pq := f128FromUint64(expectedSize).div(f128FromUint64(totalMoney)) // p/(1-p)",
+ },
+ {name: "distribution-pmf0-power", old: "pmf0 := qf.intPow(money)", new: "pmf0 := qf.intPow(money + 1)"},
+ {name: "distribution-cum-zero", old: "pmf: pmf0, cum: pmf0, at: 0", new: "pmf: pmf0, cum: f128{}, at: 0"},
+ {name: "cdf-skip-index", old: "b.at++", new: "b.at += 2"},
+ {name: "cdf-recurrence-off-by-one", old: "f128FromUint64(b.money - b.at + 1)", new: "f128FromUint64(b.money - b.at)"},
+ {
+ name: "walk-start-at-one",
+ old: "for j := uint64(0); j < money; j++ {\n\t\tboundary := dist.cdf(j)",
+ new: "for j := uint64(1); j < money; j++ {\n\t\tboundary := dist.cdf(j)",
+ },
+ {name: "walk-strict-compare", old: "P(X <= j)\n\t\tif ratio.cmp(boundary) <= 0 {", new: "P(X <= j)\n\t\tif ratio.cmp(boundary) < 0 {"},
+ {name: "freeze-fire-on-change", old: "if b.cum.cmp(cumPrev) == 0 && b.pmf.cmp(pmfPrev) < 0 {", new: "if b.cum.cmp(cumPrev) != 0 && b.pmf.cmp(pmfPrev) < 0 {"},
+ {
+ name: "freeze-nonstrict-pmf", old: "b.pmf.cmp(pmfPrev) < 0 {", new: "b.pmf.cmp(pmfPrev) <= 0 {",
+ equivalent: "an equal pmf whose add was a no-op still freezes cum forever: the step factor is non-increasing, so later pmfs stay <= this one and later adds stay no-ops",
+ },
+ {name: "divstep-cap-boundary", old: "if uHi >= v1 {", new: "if uHi > v1 {"},
+ {
+ name: "divstep-refine-nonstrict", old: "if hi > rhat || (hi == rhat && lo > uLo) {", new: "if hi > rhat || (hi == rhat && lo >= uLo) {",
+ equivalent: "the extra decrement fires only on exact quotient digits, whose zero-padded tail makes the remaining walk produce capped all-ones digits ending with remainder exactly V, and the final round-up carries the representation back to Q; the tie-sensitive case would need an odd exact 129-bit quotient, impossible since the dividend's 2^128 factor cannot come from M_b",
+ },
+ {name: "div-mantissa-shift", old: "mantHi := q2<<63 | q1>>1", new: "mantHi := q2<<62 | q1>>1"},
+ {name: "div-roundbit-q2", old: "round := q0&1 != 0", new: "round := false"},
+ {
+ name: "div-tie-nonstrict", old: "if carry != 0 || dblHi > v1 || (dblHi == v1 && dblLo >= v0) {", new: "if carry != 0 || dblHi > v1 || (dblHi == v1 && dblLo > v0) {",
+ equivalent: "exact ties are impossible in div: 2*rem == M_b requires 2^129 to divide M_b*(2q+1), and M_b < 2^128 with 2q+1 odd cannot supply it, so the >= vs > distinction is unreachable",
+ },
+ {name: "norm128-direction", old: "exp -= int64(s)", new: "exp += int64(s)"},
+ {name: "fromuint64-exp", old: "return f128{u << uint(s), 0, -int64(s) - 64}", new: "return f128{u << uint(s), 0, -int64(s) - 63}"},
+ {name: "shl256-fill-shift", old: "out[i] |= in[src+1] >> (64 - shift)", new: "out[i] |= in[src+1] >> (63 - shift)"},
+ {name: "cmp-lo-invert", old: "if a.lo < b.lo {", new: "if a.lo > b.lo {"},
+
+ // Test the tests: corrupt each major in-process oracle independently. The
+ // exact-integer, exact-CDF, high-precision, Arb, metamorphic, and production
+ // layers should make these harness defects observable rather than allowing
+ // all references to agree with one another accidentally.
+ {
+ name: "big-oracle-pmf0-power",
+ target: "f128_test.go",
+ old: "pmf := bigIntPow(q, money, prec) // (1-p)^money",
+ new: "pmf := bigIntPow(q, money+1, prec) // (1-p)^money",
+ },
+ {
+ name: "big-oracle-recurrence-off-by-one",
+ target: "f128_test.go",
+ old: "new(big.Float).SetPrec(prec).SetUint64(money-j+1),",
+ new: "new(big.Float).SetPrec(prec).SetUint64(money-j),",
+ },
+ {
+ name: "digest-oracle-denominator",
+ target: "f128_test.go",
+ old: "denominator.Sub(denominator, big.NewInt(1))",
+ new: "denominator.Sub(denominator, big.NewInt(2))",
+ },
+ {
+ name: "highprec-oracle-recurrence-off-by-one",
+ target: "f128_exact_test.go",
+ old: "new(big.Float).SetPrec(prec).SetUint64(money-j+1),",
+ new: "new(big.Float).SetPrec(prec).SetUint64(money-j),",
+ },
+ {
+ name: "frozen-sliver-invert",
+ target: "f128_exact_test.go",
+ old: "return gap.Cmp(bound) < 0",
+ new: "return gap.Cmp(bound) > 0",
+ },
+ {
+ name: "exact-cdf-binomial-off-by-one",
+ target: "f128_exact_test.go",
+ old: "term := new(big.Int).Binomial(int64(money), int64(i))",
+ new: "term := new(big.Int).Binomial(int64(money-1), int64(i))",
+ },
+ {
+ name: "exact-select-strict-boundary",
+ target: "f128_exact_test.go",
+ old: "if lhs.Cmp(rhs) <= 0 {",
+ new: "if lhs.Cmp(rhs) < 0 {",
+ },
+ {
+ name: "straddle-wrong-ulp",
+ target: "f128_exact_test.go",
+ old: "step := new(big.Int).Rsh(den, uint(128-c.MantExp(nil)))",
+ new: "step := new(big.Int).Rsh(den, uint(127-c.MantExp(nil)))",
+ },
+}
+
+func main() {
+ const defaultTarget = "f128.go"
+ originals := make(map[string][]byte)
+ for i := range mutants {
+ if mutants[i].target == "" {
+ mutants[i].target = defaultTarget
+ }
+ if _, ok := originals[mutants[i].target]; ok {
+ continue
+ }
+ orig, err := os.ReadFile(mutants[i].target)
+ if err != nil {
+ fmt.Fprintln(os.Stderr, err)
+ os.Exit(1)
+ }
+ originals[mutants[i].target] = orig
+ }
+ restore := func(target string) {
+ if err := os.WriteFile(target, originals[target], 0o644); err != nil {
+ fmt.Fprintf(os.Stderr, "restore %s: %v\n", target, err)
+ os.Exit(1)
+ }
+ }
+ restoreAll := func() {
+ for target := range originals {
+ restore(target)
+ }
+ }
+ defer restoreAll()
+
+ killArgs := []string{
+ "test", "-count=1",
+ "-run", "TestF128Digest|TestF128Primitive|TestF128OpsExact|TestF128ConversionsExact|TestDivStepExact|TestDivUVsBig|TestDivVsBig|TestSelectF128|TestSelectHighPrecision|TestRapid|Fuzz",
+ "-rapid.checks=2000",
+ // killed mutants fail the rapid tests by design; don't litter
+ // testdata/rapid with failfiles for them
+ "-rapid.nofailfile",
+ }
+
+ killed, expectedSurvivors := 0, 0
+ var unexpected []string
+ for _, m := range mutants {
+ src := string(originals[m.target])
+ if c := strings.Count(src, m.old); c != 1 {
+ restoreAll()
+ fmt.Fprintf(os.Stderr, "mutant %s: target string occurs %d times in %s, need exactly 1\n", m.name, c, m.target)
+ os.Exit(1)
+ }
+ if err := os.WriteFile(m.target, []byte(strings.Replace(src, m.old, m.new, 1)), 0o644); err != nil {
+ fmt.Fprintln(os.Stderr, err)
+ os.Exit(1)
+ }
+ out, testErr := exec.Command("go", killArgs...).CombinedOutput()
+ restore(m.target)
+ if strings.Contains(string(out), "[build failed]") {
+ fmt.Fprintf(os.Stderr, "mutant %s: does not compile, fix the table\n%s\n", m.name, out)
+ os.Exit(1)
+ }
+ switch {
+ case testErr != nil && m.equivalent == "":
+ killed++
+ fmt.Printf("KILLED %s\n", m.name)
+ case testErr != nil && m.equivalent != "":
+ fmt.Printf("UNEXPECTED KILL of equivalent mutant %s: the equivalence proof is wrong or the code changed\n", m.name)
+ unexpected = append(unexpected, m.name+" (equivalent mutant was killed)")
+ case testErr == nil && m.equivalent != "":
+ expectedSurvivors++
+ fmt.Printf("SURVIVED %s (expected: %s)\n", m.name, m.equivalent)
+ default:
+ fmt.Printf("SURVIVED %s <-- TEST GAP\n", m.name)
+ unexpected = append(unexpected, m.name)
+ }
+ }
+ fmt.Printf("\n%d mutants: %d killed, %d expected-equivalent survivors, %d UNEXPECTED\n",
+ len(mutants), killed, expectedSurvivors, len(unexpected))
+ if len(unexpected) > 0 {
+ for _, n := range unexpected {
+ fmt.Printf(" unexpected: %s\n", n)
+ }
+ os.Exit(1)
+ }
+}
diff --git a/oracle_hardening_test.go b/oracle_hardening_test.go
new file mode 100644
index 0000000..b36572c
--- /dev/null
+++ b/oracle_hardening_test.go
@@ -0,0 +1,147 @@
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
+// This file is part of go-algorand
+//
+// go-algorand is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Affero General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// go-algorand is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with go-algorand. If not, see .
+
+package sortition
+
+import (
+ "encoding/json"
+ "math/big"
+ "os"
+ "testing"
+)
+
+// TestSelectHighPrecisionCertifiedVectors makes the ordinary big.Float
+// tolerance oracle earn its precision assumption: 256, 512, and 1024-bit
+// walks must converge to every independently Arb-certified large-money
+// quantile. This catches precision-sensitive or shared harness defects even
+// when the production result itself still matches the checked-in want value.
+func TestSelectHighPrecisionCertifiedVectors(t *testing.T) {
+ raw, err := os.ReadFile("testdata/f128_arb_certificates.json")
+ if err != nil {
+ t.Fatal(err)
+ }
+ var certificates arbCertificateFile
+ if err := json.Unmarshal(raw, &certificates); err != nil {
+ t.Fatal(err)
+ }
+ if len(certificates.Vectors) < 10 {
+ t.Fatalf("too few certified vectors: %d", len(certificates.Vectors))
+ }
+ for _, vector := range certificates.Vectors {
+ t.Run(vector.Label, func(t *testing.T) {
+ var d Digest
+ decoded, ok := new(big.Int).SetString(vector.Digest, 16)
+ if !ok || decoded.Sign() < 0 || decoded.BitLen() > 256 {
+ t.Fatalf("invalid digest %q", vector.Digest)
+ }
+ decoded.FillBytes(d[:])
+ for _, prec := range []uint{256, 512, 1024} {
+ if got := selectAtPrecision(vector.Money, vector.TotalMoney, vector.ExpectedSize, d, prec); got != vector.Want {
+ t.Fatalf("%d-bit oracle selected %d, Arb certificate wants %d", prec, got, vector.Want)
+ }
+ }
+ })
+ }
+}
+
+func exactDigestBoundary(point, digestDen, cdfDen *big.Int, cdfNums []*big.Int) bool {
+ lhs := new(big.Int).Mul(point, cdfDen)
+ for _, cdfNum := range cdfNums {
+ if lhs.Cmp(new(big.Int).Mul(cdfNum, digestDen)) == 0 {
+ return true
+ }
+ }
+ return false
+}
+
+// TestSelectHighPrecisionExhaustiveSmallDomain independently checks the
+// high-precision harness over a complete finite parameter/boundary grid. At
+// non-equality points the 512/1024-bit walks must equal exact rational
+// mathematics; on exact boundaries their rounded recurrence may land one
+// inclusive boundary later, but the two high precisions must still converge.
+func TestSelectHighPrecisionExhaustiveSmallDomain(t *testing.T) {
+ digestDen := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+ const maxMoney = uint64(10)
+ const maxTotal = uint64(20)
+ var parameters, candidates, asserted, exactBoundaries, coarseDifferences uint64
+
+ for money := uint64(1); money <= maxMoney; money++ {
+ for total := uint64(2); total <= maxTotal; total++ {
+ cdfDen := new(big.Int).Exp(new(big.Int).SetUint64(total), new(big.Int).SetUint64(money), nil)
+ for expected := uint64(1); expected < total; expected++ {
+ parameters++
+ cdfNums := make([]*big.Int, money)
+ points := map[string]*big.Int{
+ "0": new(big.Int),
+ digestDen.Text(16): new(big.Int).Set(digestDen),
+ }
+ for j := uint64(0); j < money; j++ {
+ cdfNum := exactCDFNum(money, total, expected, j)
+ cdfNums[j] = cdfNum
+ scaled := new(big.Int).Mul(cdfNum, digestDen)
+ at := new(big.Int).Quo(scaled, cdfDen)
+ if at.Sign() > 0 {
+ below := new(big.Int).Sub(new(big.Int).Set(at), big.NewInt(1))
+ points[below.Text(16)] = below
+ }
+ points[at.Text(16)] = new(big.Int).Set(at)
+ if at.Cmp(digestDen) < 0 {
+ above := new(big.Int).Add(new(big.Int).Set(at), big.NewInt(1))
+ points[above.Text(16)] = above
+ }
+ }
+
+ for _, point := range points {
+ candidates++
+ d := digestFromBigInt(t, point, digestDen)
+ if inFrozenSliver(money, digestRatioBig(d, 512)) {
+ continue
+ }
+ asserted++
+ want := exactSelectFromCDFNumerators(money, cdfNums, cdfDen, digestDen, point)
+ got256 := selectAtPrecision(money, total, expected, d, 256)
+ got512 := selectAtPrecision(money, total, expected, d, 512)
+ got1024 := selectAtPrecision(money, total, expected, d, 1024)
+ if got512 != got1024 {
+ t.Fatalf("oracle failed to converge: 512=%d 1024=%d (money=%d total=%d expected=%d digest=%x)",
+ got512, got1024, money, total, expected, d)
+ }
+ if got256 != got512 {
+ coarseDifferences++
+ if uint64Distance(got256, got512) > 1 {
+ t.Fatalf("256-bit oracle differs by more than one boundary: 256=%d 512=%d", got256, got512)
+ }
+ }
+ if exactDigestBoundary(point, digestDen, cdfDen, cdfNums) {
+ exactBoundaries++
+ if uint64Distance(got1024, want) > 1 {
+ t.Fatalf("1024-bit oracle=%d vs exact boundary=%d", got1024, want)
+ }
+ } else if got1024 != want {
+ t.Fatalf("1024-bit oracle=%d vs exact=%d (money=%d total=%d expected=%d digest=%x)",
+ got1024, want, money, total, expected, d)
+ }
+ }
+ }
+ }
+ }
+ if parameters < 1_000 || asserted < 20_000 || exactBoundaries == 0 || coarseDifferences < 1_000 {
+ t.Fatalf("anti-vacuity: parameters=%d candidates=%d asserted=%d exactBoundaries=%d coarseDifferences=%d",
+ parameters, candidates, asserted, exactBoundaries, coarseDifferences)
+ }
+ t.Logf("oracle grid: parameters=%d candidates=%d asserted=%d exactBoundaries=%d 256vs512=%d",
+ parameters, candidates, asserted, exactBoundaries, coarseDifferences)
+}
diff --git a/sortition.go b/sortition.go
index 8ea34b9..77dd5cd 100644
--- a/sortition.go
+++ b/sortition.go
@@ -1,4 +1,4 @@
-// Copyright (C) 2019-2023 Algorand, Inc.
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
// This file is part of go-algorand
//
// go-algorand is free software: you can redistribute it and/or modify
@@ -58,6 +58,96 @@ func Select(money uint64, totalMoney uint64, expectedSize float64, vrfOutput Dig
return uint64(C.sortition_binomial_cdf_walk(C.double(binomialN), C.double(binomialP), C.double(cratio), C.uint64_t(money)))
}
+// SelectF128MaxMoney bounds SelectF128's money argument. The f128 stored
+// exponent of a value v is about log2(v) - 127 (the mantissa is normalized to
+// [2^127, 2^128)), and the smallest representable 1-p exceeds 2^-64, so
+// pmf(0) = (1-p)^money carries a stored exponent no lower than
+// -64*money - 128, and the worst intermediate -- the a.exp+b.exp sum inside a
+// multiply, whose value parts total at most money -- stays above
+// -64*money - 256. With money < 2^56 every such quantity is bounded by
+// ~2^62+2^8 in magnitude, comfortably inside int64. (A 2^57 bound is NOT
+// safe: money = 2^57-1 with 1-p = 1/(2^64-1) needs a stored exponent near
+// -2^63-63 and wraps.)
+//
+// The bound is ~7x (about 2.8 bits) above Algorand's 10^16 microalgo supply.
+// Behavior above it is undefined, and SelectF128 does not check it at
+// runtime. Consumers should assert their supply invariants against this
+// constant in a test, so a future economics change fails loudly there instead
+// of silently misrounding consensus.
+const SelectF128MaxMoney = uint64(1) << 56
+
+// SelectF128 is a deterministic sortition function. It evaluates both the VRF
+// ratio and binomial CDF at f128 precision using software integer arithmetic, so
+// its result is bit-reproducible across platforms. money must be below
+// SelectF128MaxMoney.
+//
+// Unlike Select, SelectF128 takes the committee size as the exact uint64 it is
+// in the protocol rather than a float64. The distribution constants are then
+// single correctly-rounded f128 divides of exact integers -- float64(totalMoney)
+// is inexact above 2^53, and a float64 p = expectedSize/totalMoney would stack
+// further roundings -- and invalid probabilities (NaN, Inf, negative,
+// fractional) are unrepresentable. p >= 1, i.e. expectedSize >= totalMoney, is
+// an exact integer comparison, handled as all probability mass at money.
+//
+// CONSENSUS / MIGRATION NOTE. SelectF128 is not bit-identical to the deployed
+// Boost-double Select. They agree on the overwhelming majority of inputs but can
+// differ at knife-edge VRF outputs near a CDF boundary. SelectF128 also preserves
+// the digest ratio at f128 precision instead of first rounding it to float64;
+// this avoids the multi-step tail divergence when a near-maximum digest rounds
+// to 1.0 in float64. The success probability is likewise formed at f128
+// precision from the integer expectedSize/totalMoney instead of a float64
+// quotient. Generic CDF-boundary differences can still change a selection by
+// one, so replacing Select with SelectF128 remains a protocol-gated,
+// network-coordinated consensus change.
+//
+// One tail edge: the top ~2^-129 of digest space rounds to an f128 ratio of
+// exactly 1.0 (only the all-0xff digest IS exactly 1.0; the rest of the
+// interval rounds up to it). With the threshold fixed at 1.0, the exact CDF is
+// < 1 for every j < money and the walk never evaluates cdf(money) == 1, so the
+// exact-CDF count is money -- and the walk returns money unless the accumulated
+// f128 CDF happens to round up to exactly 1.0 at an earlier j (a rounding
+// artifact), in which case it returns that j. Both outcomes occur, decided
+// per-distribution at ulp granularity: the money=1954 case in
+// TestSelectF128RatioExactlyOne stops at j=3, while the same distribution with
+// total=2_000_000_000_000_000 falls through to money. Both match the 128-bit
+// big.Float oracle. Boost's double CDF -- evaluated independently per j via
+// ibetac rather than accumulated -- can also saturate to 1.0 on its far coarser
+// grid, potentially at a different (typically earlier) j. At these near-maximum
+// digests the two implementations can therefore return wildly different counts:
+// one may stop within a few steps of the binomial tail while the other returns
+// the full trial count money. A uniform VRF output lands in this interval with
+// probability about 2^-129 per credential. The event is possible, but the
+// consensus threat model treats it as negligible and assumes the registered
+// key and unpredictable seed prevent an adversary from targeting it. This is
+// therefore a documented statistical edge rather than a case special-cased by
+// the implementation.
+//
+// The same holds in a wider sliver just below 1.0. pmf(0) = (1-p)^money
+// amplifies the 2^-129 rounding of 1-p by up to the trial count, and pmf(0)
+// scales every PMF term, so the accumulated CDF settles at a plateau that can
+// sit as much as ~money*2^-129 below 1 (~2^-78 at 2e15 microalgos of stake,
+// and ~2^-76 at the 10^16-microalgo mainnet supply ceiling). A digest ratio
+// between that plateau and 1.0 sits above every boundary without rounding to
+// 1.0; the walk detects the frozen CDF and immediately returns money, the same
+// result the plain walk would reach after up to money no-op iterations.
+//
+// These outputs are possible under current go-algorand committee and balance
+// bounds; TestSelectF128CurrentConsensusFrozenTail pins examples from the base
+// account minimum through the mainnet supply ceiling. For an account with
+// stake m, the affected interval is approximately m*2^-129, and summing that
+// first-order bound over all online accounts gives approximately
+// totalMoney*2^-129 per committee selection, independent of how stake is
+// split. The consensus rationale for accepting the edge is probabilistic, not
+// impossibility: registered keys and an unpredictable seed are assumed to
+// prevent targeting the interval. Within it the count is DEFINED as money
+// rather than the exact binomial-tail crossing. Computing pmf(0) with guard
+// bits could narrow the interval, at the cost of additional consensus-critical
+// arithmetic and audit surface.
+func SelectF128(money uint64, totalMoney uint64, expectedSize uint64, vrfOutput Digest) uint64 {
+ ratio := f128FromDigestRatio(vrfOutput)
+ return binomialCDFWalkF128(expectedSize, totalMoney, ratio, money)
+}
+
func init() {
var b int
var err error
diff --git a/sortition_test.go b/sortition_test.go
index 8b6fc07..9550312 100644
--- a/sortition_test.go
+++ b/sortition_test.go
@@ -1,4 +1,4 @@
-// Copyright (C) 2019-2023 Algorand, Inc.
+// Copyright (C) 2019-2026 Algorand Foundation Ltd.
// This file is part of go-algorand
//
// go-algorand is free software: you can redistribute it and/or modify
@@ -34,6 +34,10 @@ func BenchmarkSortition(b *testing.B) {
}
func TestSortitionBasic(t *testing.T) {
+ // The 2% tolerance below is only ~2 sigma for N=1000, so a randomly
+ // seeded rng flakes about 4% of the time (21/500 in one measurement);
+ // seed deterministically to keep the sanity check stable.
+ rng := rand.New(rand.NewSource(42))
hitcount := uint64(0)
const N = 1000
const expectedSize = 20
@@ -41,7 +45,7 @@ func TestSortitionBasic(t *testing.T) {
const totalMoney = 200
for i := 0; i < N; i++ {
var vrfOutput Digest
- rand.Read(vrfOutput[:])
+ rng.Read(vrfOutput[:])
selected := Select(myMoney, totalMoney, expectedSize, vrfOutput)
hitcount += selected
}
@@ -58,3 +62,16 @@ func TestSortitionBasic(t *testing.T) {
t.Errorf("wanted %d selections but got %d, d=%d, maxd=%d", expected, hitcount, d, maxd)
}
}
+
+// TestSelectF128MaxMoneyHeadroom is a tripwire on the exported domain bound:
+// it must keep at least two bits of headroom over the 10^16 microalgo supply
+// (the bound itself sits ~7x above it). Consumers are expected to run the
+// mirror-image check (their maximum stake against SelectF128MaxMoney) in
+// their own test suites.
+func TestSelectF128MaxMoneyHeadroom(t *testing.T) {
+ const mainnetSupply = uint64(10_000_000_000_000_000)
+ if SelectF128MaxMoney < 4*mainnetSupply {
+ t.Fatalf("SelectF128MaxMoney=%d leaves less than 2 bits of headroom over the %d supply",
+ SelectF128MaxMoney, mainnetSupply)
+ }
+}
diff --git a/testdata/f128_arb_certificates.json b/testdata/f128_arb_certificates.json
new file mode 100644
index 0000000..dc5b170
--- /dev/null
+++ b/testdata/f128_arb_certificates.json
@@ -0,0 +1,139 @@
+{
+ "format": 1,
+ "generator": "python-flint 0.6.0 / Arb regularized incomplete beta",
+ "semantics": "exact binomial CDF; frozen-tail definition is tested separately",
+ "vectors": [
+ {
+ "cdf_previous_upper": {"exponent": -128, "mantissa": "156599776420698413224136027028177551535"},
+ "cdf_selected_lower": {"exponent": -128, "mantissa": "183682590500240050239238580403590659885"},
+ "digest": "8000000000000000000000000000000000000000000000000000000000000000",
+ "expected_size": 100,
+ "label": "toy symmetric median",
+ "money": 100,
+ "precision_bits": 128,
+ "total_money": 200,
+ "want": 50
+ },
+ {
+ "cdf_previous_upper": {"exponent": 0, "mantissa": "0"},
+ "cdf_selected_lower": {"exponent": -127, "mantissa": "132489513636925215833268058479265445677"},
+ "digest": "5555555555555555555555555555555555555555555555555555555555555555",
+ "expected_size": 1,
+ "label": "toy low probability",
+ "money": 250,
+ "precision_bits": 128,
+ "total_money": 1000,
+ "want": 0
+ },
+ {
+ "cdf_previous_upper": {"exponent": -129, "mantissa": "336977676150904036158535031540153095825"},
+ "cdf_selected_lower": {"exponent": -128, "mantissa": "173445286874191883749825300413014048653"},
+ "digest": "8000000000000000000000000000000000000000000000000000000000000000",
+ "expected_size": 1500,
+ "label": "realistic half stake median",
+ "money": 1000000000000000,
+ "precision_bits": 128,
+ "total_money": 2000000000000000,
+ "want": 750
+ },
+ {
+ "cdf_previous_upper": {"exponent": 0, "mantissa": "0"},
+ "cdf_selected_lower": {"exponent": -128, "mantissa": "340282349906820542765149466846391461685"},
+ "digest": "8000000000000000000000000000000000000000000000000000000000000000",
+ "expected_size": 5000,
+ "label": "base minimum ordinary ratio",
+ "money": 100000,
+ "precision_bits": 128,
+ "total_money": 10000000000000000,
+ "want": 0
+ },
+ {
+ "cdf_previous_upper": {"exponent": -120, "mantissa": "1329227995776608823028674855379832097"},
+ "cdf_selected_lower": {"exponent": -128, "mantissa": "340282366920933148852584800746384209155"},
+ "digest": "fffffffffeffffffffffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 5000,
+ "label": "payout minimum upper tail",
+ "money": 30000000000,
+ "precision_bits": 128,
+ "total_money": 10000000000000000,
+ "want": 5
+ },
+ {
+ "cdf_previous_upper": {"exponent": -255, "mantissa": "57896044618658097653266946477002718428901149743768787149754063757417995575359"},
+ "cdf_selected_lower": {"exponent": -256, "mantissa": "115792089237316195329040965944463479586722667708360581904416723911341536385337"},
+ "digest": "ffffffffffffffefffffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 1500,
+ "label": "online stake 1500 tail",
+ "money": 2000000000000000,
+ "precision_bits": 256,
+ "total_money": 2000000000000000,
+ "want": 1852
+ },
+ {
+ "cdf_previous_upper": {"exponent": -254, "mantissa": "28948022309329048449299151382105534347536068724501085061549015824078971801291"},
+ "cdf_selected_lower": {"exponent": -253, "mantissa": "14474011154664524245255387440208678493560896567140995207334418326168362551413"},
+ "digest": "fffffffffffffeffffffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 6000,
+ "label": "online stake 6000 tail",
+ "money": 2000000000000000,
+ "precision_bits": 256,
+ "total_money": 2000000000000000,
+ "want": 6667
+ },
+ {
+ "cdf_previous_upper": {"exponent": -1023, "mantissa": "89884656743115795385242970878554523733564596984016408477625334557586105188621702535322718752296048392454972761686449997199852199991795019341736990637396753510447422984243425760068850563107190941993659593943613760903625081088191641240624303691268980002798728649789643649421203935109649491495961034724279304257"},
+ "cdf_selected_lower": {"exponent": -1023, "mantissa": "89884656743115795385388413133669867415371296882976122087520179112241806050593573000179506316104827266705484847876967722941298388590715560014325342305736617517667199330567254403677229988571370026180049092144739446189306956680122827434420356103071786941032054663434032334314616392702364816670890868039570748871"},
+ "digest": "ffffffffffffffffbfffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 5000,
+ "label": "supply ceiling 5000 tail",
+ "money": 10000000000000000,
+ "precision_bits": 1024,
+ "total_money": 10000000000000000,
+ "want": 5667
+ },
+ {
+ "cdf_previous_upper": {"exponent": -128, "mantissa": "340282366920938462250232204093755176997"},
+ "cdf_selected_lower": {"exponent": -126, "mantissa": "85070591730234615630680912901280710237"},
+ "digest": "ffffffffffffffbfffffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 2990,
+ "label": "one third supply 2990 tail",
+ "money": 3333333333333333,
+ "precision_bits": 128,
+ "total_money": 10000000000000000,
+ "want": 1281
+ },
+ {
+ "cdf_previous_upper": {"exponent": -512, "mantissa": "13407807929942495810394297316189667165625824202373296128999426718380328133165363386009866257047385263178159879922554263125227642700345201700995329471586391"},
+ "cdf_selected_lower": {"exponent": -512, "mantissa": "13407807929942597033636683745655549224439377075387229830050959142212430928372199307480576084781478841759129035867135163847950930971121237959356775927134835"},
+ "digest": "ffffffffffffbfffffffffffffffffffffffffffffffffffffffffffffffffff",
+ "expected_size": 1,
+ "label": "domain edge tiny probability",
+ "money": 72057594037927935,
+ "precision_bits": 512,
+ "total_money": 18446744073709551615,
+ "want": 5
+ },
+ {
+ "cdf_previous_upper": {"exponent": -129, "mantissa": "267764159092786813970757076638904795559"},
+ "cdf_selected_lower": {"exponent": 0, "mantissa": "1"},
+ "digest": "8000000000000000000000000000000000000000000000000000000000000000",
+ "expected_size": 6000,
+ "label": "probability near one",
+ "money": 3000,
+ "precision_bits": 128,
+ "total_money": 6001,
+ "want": 3000
+ },
+ {
+ "cdf_previous_upper": {"exponent": -129, "mantissa": "335325780441179917990356248262715260671"},
+ "cdf_selected_lower": {"exponent": -128, "mantissa": "172619476700348505274287757720384776963"},
+ "digest": "8000000000000000000000000000000000000000000000000000000000000000",
+ "expected_size": 9223372036854775807,
+ "label": "large balanced probability",
+ "money": 3000,
+ "precision_bits": 128,
+ "total_money": 18446744073709551615,
+ "want": 1500
+ }
+ ]
+}
diff --git a/tools/generate_arb_oracle.py b/tools/generate_arb_oracle.py
new file mode 100644
index 0000000..44f6384
--- /dev/null
+++ b/tools/generate_arb_oracle.py
@@ -0,0 +1,225 @@
+#!/usr/bin/env python3
+"""Generate rigorous large-money binomial-CDF certificates with Arb.
+
+This is deliberately an offline oracle, not a dependency of `go test`.
+python-flint exposes Arb ball arithmetic; its regularized incomplete beta
+implementation is structurally independent of SelectF128's q^money-seeded
+forward PMF recurrence. Each emitted vector carries dyadic Arb endpoints that
+prove, using exact integer comparison,
+
+ CDF(j-1) < digest/(2^256-1) <= CDF(j).
+
+Setup and regenerate from the repository root:
+
+ python3 -m venv .venv-oracle
+ .venv-oracle/bin/pip install -r tools/requirements-oracle.txt
+ .venv-oracle/bin/python tools/generate_arb_oracle.py \
+ > testdata/f128_arb_certificates.json
+
+Verify the checked-in corpus without rewriting it:
+
+ .venv-oracle/bin/python tools/generate_arb_oracle.py \
+ --check testdata/f128_arb_certificates.json
+"""
+
+import argparse
+import difflib
+import json
+import sys
+import time
+from dataclasses import dataclass
+
+import flint
+from flint import arb, ctx, fmpq
+
+
+DIGEST_DEN = (1 << 256) - 1
+MAX_MONEY = 1 << 56
+PRECISIONS = (128, 256, 512, 1024, 2048, 4096)
+
+
+@dataclass(frozen=True)
+class Case:
+ label: str
+ money: int
+ total: int
+ expected: int
+ digest: int
+
+
+def max_minus_power(bit: int) -> int:
+ return DIGEST_DEN - (1 << bit)
+
+
+CASES = (
+ Case("toy symmetric median", 100, 200, 100, 1 << 255),
+ Case("toy low probability", 250, 1000, 1, DIGEST_DEN // 3),
+ Case("realistic half stake median", 1_000_000_000_000_000, 2_000_000_000_000_000, 1500, 1 << 255),
+ Case("base minimum ordinary ratio", 100_000, 10_000_000_000_000_000, 5000, 1 << 255),
+ Case("payout minimum upper tail", 30_000_000_000, 10_000_000_000_000_000, 5000, max_minus_power(216)),
+ Case("online stake 1500 tail", 2_000_000_000_000_000, 2_000_000_000_000_000, 1500, max_minus_power(196)),
+ Case("online stake 6000 tail", 2_000_000_000_000_000, 2_000_000_000_000_000, 6000, max_minus_power(200)),
+ Case("supply ceiling 5000 tail", 10_000_000_000_000_000, 10_000_000_000_000_000, 5000, max_minus_power(190)),
+ Case("one third supply 2990 tail", 10_000_000_000_000_000 // 3, 10_000_000_000_000_000, 2990, max_minus_power(198)),
+ Case("domain edge tiny probability", MAX_MONEY - 1, DIGEST_DEN >> 192, 1, max_minus_power(206)),
+ Case("probability near one", 3000, 6001, 6000, 1 << 255),
+ Case("large balanced probability", 3000, DIGEST_DEN >> 192, (DIGEST_DEN >> 192) // 2, 1 << 255),
+)
+
+
+def dyadic(x: arb) -> tuple[int, int]:
+ """Return the exact dyadic point x as (integer mantissa, exponent)."""
+ mantissa, exponent = x.man_exp()
+ return int(mantissa), int(exponent)
+
+
+def compare_rat_dyadic(num: int, den: int, point: tuple[int, int]) -> int:
+ """Compare num/den with mantissa*2^exponent using integer arithmetic."""
+ mantissa, exponent = point
+ if exponent >= 0:
+ left, right = num, (mantissa << exponent) * den
+ else:
+ left, right = num << (-exponent), mantissa * den
+ return (left > right) - (left < right)
+
+
+def cdf_ball(case: Case, j: int, precision: int) -> arb:
+ ctx.prec = precision
+ if j < 0:
+ return arb(0)
+ if j >= case.money:
+ return arb(1)
+ q = arb(fmpq(case.total - case.expected, case.total))
+ return q.beta_lower(case.money - j, j + 1, regularized=True)
+
+
+def classify(case: Case, j: int) -> bool:
+ """Return whether exact ratio <= CDF(j), raising precision as needed."""
+ for precision in PRECISIONS:
+ cdf = cdf_ball(case, j, precision)
+ lower, upper = dyadic(cdf.lower()), dyadic(cdf.upper())
+ if compare_rat_dyadic(case.digest, DIGEST_DEN, lower) <= 0:
+ return True
+ if compare_rat_dyadic(case.digest, DIGEST_DEN, upper) > 0:
+ return False
+ raise RuntimeError(f"could not separate ratio from CDF({j}) for {case.label}")
+
+
+def select(case: Case) -> int:
+ if not 0 < case.expected < case.total:
+ raise ValueError(f"certificate case must have 0 < expected < total: {case.label}")
+ if not 0 < case.digest <= DIGEST_DEN:
+ raise ValueError(f"certificate case must have 0 < digest <= max: {case.label}")
+ # Do not binary-search the entire [0,n] interval immediately. For the
+ # protocol-shaped tiny probabilities, asking incomplete beta for j near
+ # n/2 is vastly harder than evaluating it around the committee-scale
+ # quantile. Bracket exponentially from the nearer tail first.
+ if 2 * case.expected <= case.total:
+ if classify(case, 0):
+ return 0
+ lo, hi = 0, 1
+ while hi < case.money and not classify(case, hi):
+ lo, hi = hi, min(case.money, 2 * hi)
+ else:
+ if not classify(case, case.money - 1):
+ return case.money
+ hi, delta = case.money - 1, 2
+ lo = max(-1, case.money - delta)
+ while lo >= 0 and classify(case, lo):
+ hi, delta = lo, 2 * delta
+ lo = max(-1, case.money - delta)
+
+ while hi - lo > 1:
+ mid = lo + (hi - lo) // 2
+ if classify(case, mid):
+ hi = mid
+ else:
+ lo = mid
+ return hi
+
+
+def certify(case: Case, selected: int) -> dict:
+ for precision in PRECISIONS:
+ previous = cdf_ball(case, selected - 1, precision).upper()
+ current = cdf_ball(case, selected, precision).lower()
+ previous_dyadic, current_dyadic = dyadic(previous), dyadic(current)
+ if (
+ compare_rat_dyadic(case.digest, DIGEST_DEN, previous_dyadic) > 0
+ and compare_rat_dyadic(case.digest, DIGEST_DEN, current_dyadic) <= 0
+ ):
+ return {
+ "label": case.label,
+ "money": case.money,
+ "total_money": case.total,
+ "expected_size": case.expected,
+ "digest": f"{case.digest:064x}",
+ "want": selected,
+ "precision_bits": precision,
+ "cdf_previous_upper": {
+ "mantissa": str(previous_dyadic[0]),
+ "exponent": previous_dyadic[1],
+ },
+ "cdf_selected_lower": {
+ "mantissa": str(current_dyadic[0]),
+ "exponent": current_dyadic[1],
+ },
+ }
+ raise RuntimeError(f"could not certify selected boundary for {case.label}")
+
+
+def main() -> None:
+ parser = argparse.ArgumentParser()
+ parser.add_argument("--case", help="generate only the case whose label contains this text")
+ parser.add_argument(
+ "--check",
+ metavar="PATH",
+ help="regenerate the full corpus and fail if its parsed JSON differs from PATH",
+ )
+ args = parser.parse_args()
+ if args.case is not None and args.check is not None:
+ parser.error("--case and --check cannot be combined")
+ cases = tuple(case for case in CASES if args.case is None or args.case in case.label)
+ if not cases:
+ parser.error("--case did not match any vector")
+
+ vectors = []
+ for case in cases:
+ started = time.monotonic()
+ print(f"certifying {case.label}...", file=sys.stderr, flush=True)
+ vectors.append(certify(case, select(case)))
+ print(f"certified {case.label} in {time.monotonic() - started:.3f}s", file=sys.stderr, flush=True)
+ output = {
+ "format": 1,
+ "generator": f"python-flint {flint.__version__} / Arb regularized incomplete beta",
+ "semantics": "exact binomial CDF; frozen-tail definition is tested separately",
+ "vectors": vectors,
+ }
+ rendered = json.dumps(output, indent=2, sort_keys=True)
+ if args.check is None:
+ print(rendered)
+ return
+
+ try:
+ with open(args.check, encoding="utf-8") as certificate_file:
+ committed = json.load(certificate_file)
+ except (OSError, json.JSONDecodeError) as error:
+ raise SystemExit(f"cannot read certificate corpus {args.check}: {error}") from error
+ if committed != output:
+ committed_rendered = json.dumps(committed, indent=2, sort_keys=True)
+ print(
+ "".join(
+ difflib.unified_diff(
+ committed_rendered.splitlines(keepends=True),
+ rendered.splitlines(keepends=True),
+ fromfile=args.check,
+ tofile="fresh Arb generation",
+ )
+ ),
+ file=sys.stderr,
+ )
+ raise SystemExit("Arb certificate corpus is stale or non-reproducible")
+ print(f"verified {len(vectors)} Arb certificates against {args.check}", file=sys.stderr)
+
+
+if __name__ == "__main__":
+ main()
diff --git a/tools/requirements-oracle.txt b/tools/requirements-oracle.txt
new file mode 100644
index 0000000..638065c
--- /dev/null
+++ b/tools/requirements-oracle.txt
@@ -0,0 +1 @@
+python-flint==0.6.0