Accuracy characterization vs mpmath (#46) - #101
Merged
Conversation
Adds the oracle/report half of #46 (SIMD accuracy characterization vs mpmath): tools/accuracy_vs_mpmath.py and docs/ACCURACY_CHARACTERIZATION.md. The sweep-generation half (tools/accuracy_sweep.cpp, tools/CMakeLists.txt) is owned by a sibling change on this branch. tools/accuracy_vs_mpmath.py: - mpmath (mp.dps=50) reference pdf/logpdf/cdf/quantile for all 19 distributions, parameter order read from each include/libstats/distributions/<dist>.h Doxygen constructor comment. - Gaussian/LogNormal CDF reuse the erfc closed form from scripts/gen_gaussian_cdf_vectors.py / gen_lognormal_cdf_vectors.py. Von Mises CDF reuses scripts/gen_vonmises_cdf_vectors.py's direct quadrature and its bit-exact double wrap_to_pi (including the PI_D anchoring rationale, copied rather than rediscovered). - Discrete distributions (poisson, binomial, negative_binomial, geometric, discrete) use mpmath's regularized incomplete gamma/beta identities for pmf/cdf; quantile by monotone integer search. - Deep-tail cdf rows (F < 1e-3) report law_frac against the achievable-accuracy law rel(F) ~ |ln F| * 2^-52 from tests/test_lognormal_cdf_accuracy.cpp's law_budget comment, kept separate from max_rel so a law-limited tail doesn't read as a defect. - NaN/Inf rows are contract-checked (NaN in -> NaN out, etc.) separately from the error statistics. - 38 oracle self-checks (2+ per distribution: median/symmetry identities and closed-form spot values) run before any row is compared; the process exits non-zero if any fails. All 38 pass. Verified against a hand-made 34-row fixture (not checked in; built with mpmath-rounded scalar/batch values so error is near the double-rounding floor, ~1e-17 to ~1e-16, across every distribution/method combination, plus a deep-tail lognormal cdf row, a quantile row with batch_bits="-", and a NaN-input row): - Every legitimate row measured max_rel in the 1e-18 to 1e-16 range, as expected for a fixture built from the oracle's own rounded output. - One row (exponential cdf, x=1.5) was deliberately sabotaged (scalar_bits inflated ~25%) and the report correctly isolated it: max_rel=0.2526 against every neighboring row's ~1e-17, with batch_vs_scalar flagging the same row (abs=0.24, rel=0.2017) since only scalar_bits was corrupted. - A separate negative test corrupted a NaN row's scalar_bits to a finite value; the contract-violation check caught it ("NaN input did not produce NaN scalar output") and excluded it from the error stats, as designed. docs/ACCURACY_CHARACTERIZATION.md: title, PROVISIONAL banner (one machine -- Zen 4/AVX-512, MSVC Release; scalar libm = MSVC UCRT; Kaby Lake/AVX2 and M1/NEON legs pending), regeneration recipe, the law-of-F caveat for deep-tail cdf columns, cross-references to the four existing pinned accuracy gates (test_trig_ulp_gates, test_vonmises_cdf_ accuracy, test_lognormal_cdf_accuracy, test_gaussian_cdf_accuracy), and an empty BEGIN/END GENERATED block for the orchestrator's verification pass to populate against the real sweep. Judgment calls (see final report for the full list): Von Mises quantile oracle bisects the mu=0-centered CDF then shifts/wraps, which mirrors the library's approach but not its 2049-point grid interpolation, so scalar-vs-oracle disagreement there is expected to include interpolation error, not just double rounding; Student's t CDF uses the standard regularized-incomplete-beta closed form (self-checked via symmetry and cdf(0)=0.5 rather than derived from first principles in this message). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds tools/accuracy_sweep.cpp: a standalone binary that evaluates every
distribution's scalar AND batch (span/SIMD, FORCE_VECTORIZED) pdf/logpdf/
cdf/quantile paths over fixed characterization grids and dumps bit-exact
results (uint64 hex of the IEEE-754 bit pattern) to a CSV, for the sibling
tools/accuracy_vs_mpmath.py oracle to compare against mpmath. This
replaces the issue's original pylibstats approach: pylibstats pins to a
released libstats and would characterize the wrong code, not the commit
under test.
Grid construction is support-aware for free: continuous x-grids are built
by calling each distribution's own getQuantile() at a fixed, symmetric,
log-spaced probability grid (down to p=1e-300 in both tails), so no
per-family tail formulas are hand-derived. Beta/uniform (the only
both-sides-bounded continuous supports) additionally get explicit
boundary +/- {1e-15, 1e-9, 1e-4} points. Discrete grids are landmark
points (0, 1, 2, mean, mean +/- t*sd for a deterministic multiplier
ladder up to 16000 sd, and support edges) driven by each distribution's
own getMean()/getVariance(), not sampling. NaN/+inf/-inf rows are added
for every continuous dist's pdf/logpdf/cdf. Batch quantile is always "-"
-- no distribution exposes a span quantile overload. 3 parameter
instances per distribution (typical / small-boundary / stressed, e.g.
von Mises kappa=100, lognormal sigma=2) within each header's documented
valid range.
Registers accuracy_sweep in tools/CMakeLists.txt via add_libstats_tool,
next to the other correctness/stress tools.
Determinism: ran twice into separate files and diffed byte-for-byte
identical (5928 rows, 0 quantile skips, MSVC Release/AVX-512). Spot-
checked by hand: gaussian pdf(0; mean=0,sigma=1) = 0x3fd9884533d43650
(0.3989422804014327 = 1/sqrt(2pi), matches to the last bit against an
independent computation); exponential cdf(x=9.992...e-16; lambda=1)
matches 1-exp(-x) to 6 ULP (expected -- the kernel uses an expm1-style
form for small-x accuracy, which is more accurate than the naive
reference, not less); gaussian quantile(p=0.5; mean=0,sigma=1) = 0x0
(median at the mean). Format cross-checked against tools/
accuracy_vs_mpmath.py's parse_csv()/DIST_PARAMS (already landed by the
sibling change): column order, "0x0" p2 sentinel for 1-parameter dists,
"-" batch sentinel, dist/method vocabulary, and the literal header line
all match its input contract.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…results Verification pass over the two agent-built halves of #46. The sweep tool survived unchanged (two-process byte-identical determinism confirmed); the mpmath oracle needed substantial hardening before it could survive the real 5928-row sweep -- every failure below was hit in sequence by actual sweep rows, none by the fixture: - Quantile inversion: mp.findroot secant stalls permanently on deep-tail plateaus (gamma at p ~ 1-1e-16: |f| ~ 1e-15 vs tol 2.6e-54). Pure log-space bisection fixed that but costs ~200 CDF evaluations/row -- hours when one evaluation is a large-parameter incomplete beta. Replaced with safeguarded false position in (log x, log F) coordinates (bracket invariant kept, anti-stagnation forcing, exact in one step on power-law tails), doubling bracket expansion with an exponent cap (gamma alpha=0.01 at p=1e-300 has its quantile near 1e-30000), and asymptotic-normal seeds where erfinv degenerates (2p-1 rounds to +-1 at dps 50 below p ~ 1e-50). - mp.betainc is unusable at large parameters: huge-a far-from-root escalates hypercomb precision for minutes then raises ValueError; large-BOTH-parameters (binomial n=1e6 at the mean) hangs even in the central region. Added _betainc_reg: far-tail lead guard (exp(lead) below every positive double; also the exact log-coordinate the solver interpolates on) + clean-room Lentz continued-fraction incomplete beta for min(a,b) >= 5000, written from the textbook recurrence, converging in O(sqrt(ab/(a+b))) mpf iterations. Anchored by new self-checks (bit-agreement vs mp.betainc at moderate parameters, I_{1/2}(a,a)=1/2 at a=1e4) and cross-validated against direct quadrature at 1e-44. - Lower gammainc raises NoConvergence (mpmath's own class, NOT a ValueError) right of the mean at large alpha (chi-squared k=1e5): route y > alpha through the upper-gamma complement; guard the lower tail with the same lead form. - Von Mises quantile findroot call passed the bisect bracket twice (TypeError on mpmath 1.4); boundary-limit fixes: beta pdf/logpdf edge tri-state (integrable singularity for shape < 1 -- the library's +inf at x=1 was correct, the oracle's -inf was not), weibull logpdf(0) tri-state, geometric logpdf k=0 at p=1 (0*log(0) NaN). - Report: +-inf inputs are scored against their universal limits (pdf->0, logpdf->-inf, cdf->0/1) instead of NaN-ing in the reference formulas; von Mises +-inf rows check scalar/batch consistency only (periodic: no limit exists; the library's saturation convention is recorded as such); non-finite library values against finite references are contract rows, never folded into error stats; batch values are now checked on non-finite-reference rows too; the full violation list is emitted as a durable appendix in the generated doc block. Results (docs/ACCURACY_CHARACTERIZATION.md, generated + Findings prose): the four pinned gates reconcile -- gaussian law_frac 1.319 vs pinned 0.287 of the 1e-15 budget (= 1.29 bare-law), lognormal 1.784 <= 2.2 equivalent, von Mises max_abs 1.7e-16/5.1e-16 vs 2.2e-16/8.9e-16 budgets. 86 contract findings: batch NaN propagation inconsistent across 8 distributions (uniform batch pdf(NaN) returns the in-support density), NaN returns at +-inf where limits exist, batch logpdf(+inf) clamp -4605.0, quantile NaN/saturation at extreme p, and large-parameter CDF accuracy limits (binomial n=1e6 off 1.3e-2 at the mean; the parked corvus adoption items #47/#52 are the remedy class for these). test_uniform_enhanced's timing assertion fails pre-existing on clean main (SIMD speedup 1.44x vs adaptive 1.8x threshold); this change set touches only tools/ and docs/. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This was referenced Aug 21, 2026
GCC 13's -Werror=strict-overflow fires a pointer-wraparound false
positive ("assuming pointer wraparound does not occur when comparing
P +- C1 with P +- C2") when it constexpr-expands std::sort's introsort
inside xConstructionPGrid()'s immediately-invoked initializer, failing
the strict-warnings runner. The grid is strictly ascending by
construction (both exponent lists descend, and 10^-0.52 ~ 0.302 < 0.35),
so the sort+unique was purely defensive; an assert loop keeps the check
without giving GCC a sort to expand. Regenerated CSV verified
byte-identical to the pre-change sweep.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #46.
What this is
The v2.3.0 accuracy-characterization deliverable: a deterministic bit-exact sweep binary (
tools/accuracy_sweep.cpp), an mp.dps-50 mpmath oracle/report generator (tools/accuracy_vs_mpmath.py, 42 unconditional self-checks), and the generateddocs/ACCURACY_CHARACTERIZATION.mdwith per-distribution error tables and a durable contract-findings appendix. Replaces the issue's pylibstats route, which pins to released v2.2.0 and would characterize the wrong code.Highlights
pdf/logpdf/cdf/quantile, support-aware grids top = 1e-300; byte-identical across independent processes.mp.betainchangs or raises (min(a,b) >= 5000— binomial n=1e6 at the mean hangs mpmath even centrally), upper-gamma complement right of the mean (chi-squared k=1e5 raisesNoConvergence), far-tail asymptotic guards, boundary-limit tri-states. Cross-validated against direct quadrature at 1e-44.logpdf(+inf)clamp −4605.0; quantile NaN/saturation at extreme p; large-parameter CDF limits (binomial n=1e6 off 1.3e-2 at the mean — the parked corvus adoption items bessel.h Tier 2 fallback (A&S polynomial) limits VonMises accuracy to ~10⁻⁷ on macOS/AppleClang #47/BinomialDistribution CDF slower than scipy; PDF near-parity — PMF summation and scalar lgamma are limiting factors #52 are the remedy class).Follow-up candidates (not filed yet)
Notes
test_uniform_enhanced's timing assertion fails pre-existing on clean main (SIMD speedup 1.44x vs adaptive 1.8x); this branch touches onlytools/anddocs/.🤖 Generated with Claude Code