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