fix: restore full precision to erf's interval constants

[agene0001]

erf_impl reconstructs erfc(z) for z >= 0.5 as exp(-z²)/z · (b + r(z)), where r is a minimax rational fit and b a per-interval constant. All 13 b constants carry only 10 significant decimal digits:

0.3440242112,   // interval [0.5, 0.75)
0.419990927,    // interval [0.75, 1.25)

Solving for the b each interval's fit implies (at 60-digit precision) gives values constant to ~19 digits across each interval — so the rational fits are good to ~1e-19 and the truncated constants were the only error source. The recovered values are exactly f32-representable: Boost declared them as float literals and the port that statrs inherits transcribed the printed decimal. This PR restores the full values, e.g. 0.3440242112 → 0.3440242111682891845703125.

It also compensates the squaring: z * z is rounded before exp sees it, costing ~z² ulps (~460 by z = 27). A Dekker product recovers the rounding error of the square exactly (Dekker rather than f64::mul_add, which is a slow software FMA on targets without the instruction, e.g. baseline x86-64).

Measured against mpmath at 45 dps over a 700-point sweep plus every interval boundary:

	before	after
erf	~1e-10 rel	0.24 median / 1.1 max ulp
erfc	~1e-10 rel	0.52 median / 2.8 max ulp
Normal::cdf	~5e-11 rel	0.32 median

This propagates to everything built on erfc — Normal cdf/sf, LogNormal, Levy. The suite masked it with epsilon = 1e-9/1e-11 tolerances, now tightened to 4 ulp of each expected value.

Six reference literals in the tests were themselves fitted to the old output and are replaced with mpmath values: three LogNormal::sf expectations (off by up to 5.1e-12) and three erfc_inv expectations — one of which, erfc_inv(1e-10), was off by 8.8e-9 and had needed epsilon = 1e-7 to pass. erf_inv/erfc_inv themselves were already ≤2 ulp; only the expectations were wrong.

🤖 Generated with Claude Code

[day01]

isnt better keep erf fixes in https://github.com/rust-lang/compiler-builtins/tree/main/libm/src/math ?

[agene0001]

Fixing it there wouldn't change anything here, because statrs doesn't use libm's erf.

src/function/erf.rs is a self-contained implementation — 28 of its own rational-approximation coefficient tables, inherited from the Math.NET Numerics port the crate began as. libm is in the dependency tree, via num-traits' libm feature, but only for basic float ops under no_std; erf.rs never calls into it. Rust's std doesn't expose erf at all, so a stats crate has to carry its own either way.

The constants in this PR are interval boundaries in that implementation — the cutoffs selecting which rational approximation to use, e.g. 0.3440242112 → 0.3440242111682891845703125. libm's erf is FDLIBM-derived: a different algorithm with different constants, so there's nothing there these correspond to.

Switching statrs to libm's erf outright would be a real option, but it's a much bigger change — different accuracy characteristics, and a behavioural change for every distribution built on it — so probably its own discussion rather than part of this fix.