| 371 | // erf(y) − ax = (1 − erfc(y)) − ax = (1 − ax) − erfc(y), |
| 372 | // computing 1 − ax once (exact) and erfc(y) directly (its own continued |
| 373 | // fraction). (NU-P1-9) |
| 374 | const q = 1 - ax; // complementary input, exact |
| 375 | const useComplement = ax > 0.9; |
| 376 | for (let i = 0; i < 5; i++) { |
| 377 | const residual = useComplement ? q - erfc(y) : erf(y) - ax; |
| 378 | y -= residual * c * Math.exp(y * y); |
| 379 | } |
| 380 | |
| 381 | return sign * y; |
| 382 | } |
| 383 | |
| 384 | /** |
| 385 | * Complementary error function, erfc(x) = 1 - erf(x), accurate to full |
| 386 | * machine (double) precision. |
| 387 | * |
| 388 | * For |x| < 2 the value is computed as `1 - erf(x)` (no significant |
| 389 | * cancellation). For larger |x|, `1 - erf(x)` would lose all precision |
| 390 | * (erf(x) ≈ 1), so erfc is computed directly from a continued fraction |
| 391 | * (DLMF 7.9.3) evaluated with the modified Lentz algorithm. |
| 392 | * |
| 393 | * References: |
| 394 | * - NIST DLMF: https://dlmf.nist.gov/7.9 |
| 395 | */ |
| 396 | export function erfc(x: number): number { |
| 397 | if (Number.isNaN(x)) return NaN; |
| 398 | if (!Number.isFinite(x)) return x > 0 ? 0 : 2; |
| 399 | if (x < 0) return 2 - erfc(-x); |
| 400 | if (x < 2) return 1 - erf(x); |
| 401 | |