fix: Stabilize beta CDF and extreme quantiles

[day01]

Summary

• use exact symmetry and Temme's uniform expansion for large-parameter beta CDFs
• replace the inverse-beta floor and unchecked value-domain Newton step with safeguarded log-domain and shape-two solvers
• certify adjacent and subnormal rounding while preserving monotonicity
• split the beta implementation into focused production modules of at most 250 lines and add regression benchmarks

Root cause

The forward path formed ill-conditioned gamma and power differences in f64; cancellation eventually drove a valid CDF outside [0, 1]. The inverse path clamped its initial estimate away from zero and iterated on an underflowing value-domain CDF, causing the fixed ~1e-17 floor, non-monotonic results, non-termination, and an unchecked error panic.

The replacement uses DLMF 8.17, DLMF 8.18.9–12, TOMS Algorithm 708, and Temme, Special Functions (1996), section 11.3.3.2. Boost-derived kernels are attributed and distributed under BSL-1.0 in THIRD_PARTY_NOTICES.md.

TDD and independent numerical references

The issue reproducers fail on upstream and pass with this change. Inputs are passed to every reference as their exact binary64 values. The columns below are signed ULP error against 500-digit MPFR/Boost.Multiprecision references; the endpoint uses the exact I_x(2,b) identity at 1100 digits. 0 is correctly rounded.

Cases: center I_0.5(1e8, 1e8); asymmetric I_x(1e8, 2e8) at x=0.33330611678068106; inverse cases (a,b,p)=(200,2,1e-60), (200,2,1e-165), (0.1,500,1e-30), and the endpoint (2,1e308,0.5).

Implementation	Center	Asymmetric	p=1e-60	p=1e-165	Tiny shape	Endpoint	Outcome
upstream statrs	outside [0,1]	+2,922,987,503	timeout	panic	+4,288,298,043,698,456,868	timeout	bugs reproduced
statrs (this PR)	0	-1	0	0	0	0	correct, 0–1 ULP
Boost.Math 1.90	+15	+15,949,984	+12	0	+3	returns 0	less accurate on five cases
SciPy 1.18.0	+15	+15,949,984	+12	0	+3	returns 0	same Boost-backed results
R 4.2.1	0	+891	0	0	-234	timeout	mixed
mpmath 1.4.1, 500 digits	timeout	timeout	0	0	0	generic inverse failed	independent cross-check
MPFR/Boost.Multiprecision, 500/1100 digits	reference	reference	reference	reference	reference	reference	rounding oracle

SciPy's beta kernels are Boost-backed, so SciPy and Boost are not independent references. mpmath is independent; MPFR/Boost.Multiprecision is used only as the high-precision arithmetic oracle, not as the Boost.Math implementation under test. A separate 1763-case endpoint audit around rounding thresholds found no error above 1 ULP; a theoretically unresolved adjacent-float overlap selects the monotone ties-to-even result without panicking.

Performance

Case	Upstream	This change	Difference
#434 center	1.909 us, outside [0,1]	3.675 ns, 0 ULP	99.8% faster
#434 large asymmetric CDF	1.899 us, +2,922,987,503 ULP	485.581 ns, -1 ULP	74.4% faster
#435 non-termination input	timeout	363.448 ns, 0 ULP	terminates
#435 panic input	panic	372.712 ns, 0 ULP	finite result
#435 tiny-shape input	5.811 us, wrong/non-monotone	628.094 ns, 0 ULP	89.2% faster
#435 large-b endpoint	timeout	121.324 ns, 0 ULP	terminates
#435 exact midpoint fallback	panic	14.727 us, 0 ULP	finite, correctly rounded
typical finite CDF	153.811 ns, -46 ULP	223.018 ns, +1 ULP	45.0% slower
typical finite inverse	507.253 ns, +10 ULP	602.560 ns, 0 ULP	18.8% slower

Implementation	Execution mode	#434 asymmetric CDF	#435 p=1e-60	#435 tiny shape	#435 endpoint
statrs (this PR)	Rust release, target-cpu=native	485.581 ns, -1 ULP	363.448 ns, 0 ULP	628.094 ns, 0 ULP	121.324 ns, 0 ULP
Boost.Math 1.90	C++17, -O3 -march=native	1.966 us, +15,949,984 ULP	321.625 ns, +12 ULP	205.119 ns, +3 ULP	411.673 ns, returns 0
R 4.2.1	scalar R call	920.000 ns, +891 ULP	1.573 us, 0 ULP	1.120 us, -234 ULP	timeout
mpmath 1.4.1, 500 digits	scalar Python reference	timeout	4.095 ms, 0 ULP	21.022 ms, 0 ULP	generic inverse failed
SciPy 1.18.0	scalar Python call	2.692 us, +15,949,984 ULP	1.099 us, +12 ULP	985.318 ns, +3 ULP	1.216 us, returns 0

Timings are scalar-call latency; a faster time does not imply a correct result—the accuracy table applies. mpmath is an accuracy oracle rather than a production-speed competitor. Batched SciPy/R throughput is measured separately and is not presented as scalar latency.

Related issues

Fixes #434 — Beta::cdf loses accuracy above shapes of ~1e4 and returns values outside [0, 1] by 1e8.

Fixes #435 — inv_beta_reg panics, does not return, and is non-monotone for small x.

[codecov]

Codecov Report

❌ Patch coverage is 95.22498% with 156 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.13%. Comparing base (92819b6) to head (240c15a).

Files with missing lines	Patch %	Lines
src/function/beta/inverse/shape_two/endpoint.rs	84.48%	18 Missing ⚠️
src/function/beta/inverse/shape_two/adjacent.rs	88.46%	15 Missing ⚠️
.../function/beta/inverse/shape_two/endpoint_fixed.rs	93.46%	13 Missing ⚠️
src/function/beta/quantile/bounds.rs	92.65%	13 Missing ⚠️
...ction/beta/inverse/shape_two/endpoint_certified.rs	92.59%	12 Missing ⚠️
src/function/beta/inverse/shape_two/value.rs	95.28%	10 Missing ⚠️
src/function/beta/log_forward.rs	94.36%	8 Missing ⚠️
src/function/beta/asymptotic.rs	94.11%	7 Missing ⚠️
src/function/beta/fraction.rs	93.63%	7 Missing ⚠️
...nction/beta/inverse/shape_two/endpoint_beta_two.rs	94.16%	7 Missing ⚠️
... and 11 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #447      +/-   ##
==========================================
+ Coverage   95.07%   95.13%   +0.06%     
==========================================
  Files          62       90      +28     
  Lines       14191    17128    +2937     
==========================================
+ Hits        13492    16295    +2803     
- Misses        699      833     +134     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[day01]

@YeungOnion i take a boost impl and few articles impl.
Comparison is correct, please take a look especially on licensing things.