Fix erf/erfc catastrophic cancellation producing impossible values (erf > 1, erfc < 0, negative Normal survival)

[gaoflow]

Problem

Special::erf / Special::erfc return mathematically impossible values in the tail, and this propagates into Normal::cdf:

Special::erf(3.99);              // 1.0000004145   erf must be < 1
Special::erfc(3.9);              // -6.78e-9       erfc must be >= 0
Special::erfc(3.99);             // -4.14e-7       erfc must be >= 0

$n = new Normal(0, 1);
1 - $n->cdf(5.6);                // -8.82e-8       a survival probability, must be >= 0
1 - $n->cdf(5.643);             // -2.08e-7       negative probability

No reference implementation is needed to see these are wrong: erf is bounded in (-1, 1), erfc in [0, 2], and a probability cannot be negative.

Root cause

All three come from one place, errorFunction():

• For 0.01 < |x| <= 4, erf is summed from the alternating Maclaurin series Σ (-1)^n x^(2n+1) / (n!(2n+1)). Its terms grow to about e^(x^2) before decaying, so near x = 4 individual terms reach ~1e7 while the sum is ~1 — losing about 11 digits to cancellation. erfc(x) = 1 - erf(x) then dumps that error straight into a value of order 1e-8, flipping its sign.
• The |x| > 4 branch used erfcAsymptoticSeries() truncated at 4 terms (it drops the +105/16 x^8 term). That divergent series is worst at the x = 4 threshold (~2.5%) and still ~1.5e-5 at x = 5.
• Normal::cdf() calls erf, so the tail CDF / survival and every Normal-derived path inherit the defect.

The existing tests missed it because none exercise 3 < |x| < 4, and the erfc tests near x = 5 use an absolute tolerance far larger than the value.

Fix

Route erf/erfc through the file's own numerically stable continued-fraction gamma, using the identity erfc(x) = Γ(1/2, x²)/√π for x >= 0 (reusing the already-tested upperIncompleteGamma) instead of an unstable/truncated series. The 8-term Maclaurin branch is kept for |x| <= 0.01, where 1 - erfc would instead cancel. erfc(-x) = 2 - erfc(x) and erf(-x) = -erf(x) stay exact. Net effect is less code, not more.

Verification

Checked against mpmath (mp.dps = 40, cross-checked at dps 35) across every dispatch branch:

• relative error < 1e-12 for both erf and erfc at every grid point, including the small-x region that was already correct (no regression);
• erf in [-1, 1] and erfc in [0, 2] on a dense sweep 0 <= x <= 8 (both signs);
• the standard-normal survival 1 - cdf(z) is non-negative for z in the previously-broken 4.5-7 range.

Added tests cover the mpmath reference values, the invariants, and the Normal tail. Full suite green locally (AllButLinearAlgebra + LinearAlgebra), plus parallel-lint, phpcs, and phpstan.

[markrogoyski]

Hi @gaoflow,

Thanks for your interest in MathPHP. And thanks for the PR to improve the code.

I'm currently in the middle of a version upgrade of the library, so please be patient as it may take me some time to get to this and decide whether to apply any fixes forward or to the current version as well.

Again, thanks for the PR. And apologies if it takes me a bit to look through this.
Mark

[gaoflow]

The red jobs aren't from this change: all 53 failures are ReflectionException from ReflectionMethod invocations on private Piecewise/Polynomial/Finance/chebyshevEval methods in files this PR doesn't touch, and #492 (a disjoint diff) fails with the identical list. The same tests pass locally with the repo's pinned PHPUnit 8.5.53, so it looks like a CI runner/PHP-build issue.