math.erfinv(1.0f) returned -infinity and erfinv(-1.0f) returned
+infinity: at value = +-1, w = -log((1-value)(1+value)) = -log(0) =
+infinity, the tail branch's leading coefficient is negative, and
`return p * value` delivered the inverted infinity. erf saturates to
exactly 1.0f at 5.543 sigma, so any interval edge beyond that reaches
this bug in ordinary use.
Guard both edges explicitly before the polynomial runs: erfinv(1) is
now +infinity, erfinv(-1) is -infinity, and anything strictly outside
[-1, 1] (including +-infinity) is NaN, which is the correct answer and
was already the fallthrough behavior. Documented the domain and the
edge behavior in a comment, matching the house style `sign` already
uses at math.cs:84. The approximation itself is unchanged.
Also recorded the true measured accuracy beside erf and erfinv: worst
absolute error 6.621e-7 and 5.066e-7 respectively, both inside the
2e-6 bar this package's tests pin.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
fc9ce970b8