1/**
2 * Standard normal CDF and inverse for the truncated-Gaussian sampler,
3 * standing in for scipy.special.ndtr/ndtri in the Python package.
4 *
5 * Φ comes from W. J. Cody's rational erfc (SPECFUN calerf), full double
6 * precision with correct relative accuracy deep in the lower tail — which is
7 * what the mirrored truncated sampler relies on. Φ⁻¹ is Acklam's
8 * approximation (~1.15e-9 relative) polished by one Halley step against this
9 * Φ, which reaches ~1e-13 everywhere the step is taken: down to x ≈ -37.5
10 * (below, exp(x²/2) overflows) and up to p = 1 - 1e-9 (beyond, the spacing
11 * of doubles near 1 limits any method to about what Acklam already gives).
12 */
14const THRESH = 0.46875
15const SQRPI = 5.6418958354775628695e-1 // 1/√π
17const ERF_A = [3.1611237438705656, 1.13864154151050156e2, 3.77485237685302021e2, 3.20937758913846947e3, 1.85777706184603153e-1]
18const ERF_B = [2.36012909523441209e1, 2.44024637934444173e2, 1.28261652607737228e3, 2.84423683343917062e3]
19const ERF_C = [5.64188496988670089e-1, 8.88314979438837594, 6.61191906371416295e1, 2.98635138197400131e2, 8.8195222124176909e2, 1.71204761263407058e3, 2.05107837782607147e3, 1.23033935479799725e3, 2.15311535474403846e-8]
20const ERF_D = [1.57449261107098347e1, 1.17693950891312499e2, 5.37181101862009858e2, 1.62138957456669019e3, 3.29079923573345963e3, 4.36261909014324716e3, 3.43936767414372164e3, 1.23033935480374942e3]
21const ERF_P = [3.05326634961232344e-1, 3.60344899949804439e-1, 1.25781726111229246e-1, 1.60837851487422766e-2, 6.58749161529837803e-4, 1.63153871373020978e-2]
22const ERF_Q = [2.56852019228982242, 1.87295284992346047, 5.27905102951428412e-1, 6.05183413124413191e-2, 2.33520497626869185e-3]
24/** erf(x) for |x| ≤ 0.46875. */
25function erfSmall(x: number): number {
26 const z = x * x
27 let num = ERF_A[4] * z
28 let den = z
29 for (let i = 0; i < 3; i++) {
30 num = (num + ERF_A[i]) * z
31 den = (den + ERF_B[i]) * z
32 }
33 return (x * (num + ERF_A[3])) / (den + ERF_B[3])
34}
36/** erfc(y) for y > 0.46875 (with the split exp(-y²) trick for accuracy). */
37function erfcLarge(y: number): number {
38 let result: number
39 if (y <= 4) {
40 let num = ERF_C[8] * y
41 let den = y
42 for (let i = 0; i < 7; i++) {
43 num = (num + ERF_C[i]) * y
44 den = (den + ERF_D[i]) * y
45 }
46 result = (num + ERF_C[7]) / (den + ERF_D[7])
47 } else {
48 const z = 1 / (y * y)
49 let num = ERF_P[5] * z
50 let den = z
51 for (let i = 0; i < 4; i++) {
52 num = (num + ERF_P[i]) * z
53 den = (den + ERF_Q[i]) * z
54 }
55 result = (z * (num + ERF_P[4])) / (den + ERF_Q[4])
56 result = (SQRPI - result) / y
57 }
58 const ysq = Math.trunc(y * 16) / 16
59 const del = (y - ysq) * (y + ysq)
60 return Math.exp(-ysq * ysq) * Math.exp(-del) * result
61}
63/** Complementary error function, double precision over the whole line. */
64export function erfc(x: number): number {
65 const y = Math.abs(x)
66 if (y <= THRESH) return 1 - erfSmall(x)
67 const tail = y > 26.6 ? 0 : erfcLarge(y)
68 return x < 0 ? 2 - tail : tail
69}
71/** Standard normal CDF Φ(t), accurate relative to its size in both tails. */
72export function ndtr(t: number): number {
73 return 0.5 * erfc(-t / Math.SQRT2)
74}
76const ACK_A = [-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2, 1.38357751867269e2, -3.066479806614716e1, 2.506628277459239]
77const ACK_B = [-5.447609879822406e1, 1.615858368580409e2, -1.556989798598866e2, 6.680131188771972e1, -1.328068155288572e1]
78const ACK_C = [-7.784894002430293e-3, -3.223964580411365e-1, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783]
79const ACK_D = [7.784695709041462e-3, 3.224671290700398e-1, 2.445134137142996, 3.754408661907416]
81const SQRT_2PI = Math.sqrt(2 * Math.PI)
83/** Standard normal quantile Φ⁻¹(p) for p in (0, 1). */
84export function ndtri(p: number): number {
85 const pLow = 0.02425
86 let x: number
87 if (p >= pLow && p <= 1 - pLow) {
88 const q = p - 0.5
89 const r = q * q
90 x =
91 ((((((ACK_A[0] * r + ACK_A[1]) * r + ACK_A[2]) * r + ACK_A[3]) * r + ACK_A[4]) * r + ACK_A[5]) * q) /
92 (((((ACK_B[0] * r + ACK_B[1]) * r + ACK_B[2]) * r + ACK_B[3]) * r + ACK_B[4]) * r + 1)
93 } else {
94 const lower = p < pLow
95 const q = Math.sqrt(-2 * Math.log(lower ? p : 1 - p))
96 x =
97 (((((ACK_C[0] * q + ACK_C[1]) * q + ACK_C[2]) * q + ACK_C[3]) * q + ACK_C[4]) * q + ACK_C[5]) /
98 ((((ACK_D[0] * q + ACK_D[1]) * q + ACK_D[2]) * q + ACK_D[3]) * q + 1)
99 if (!lower) x = -x
100 }
101 // One Halley step against the accurate Φ (see the header for the guards).
102 if (p <= 1 - 1e-9 && x > -37.5) {
103 const e = ndtr(x) - p
104 const u = e * SQRT_2PI * Math.exp((x * x) / 2)
105 x -= u / (1 + (x * u) / 2)
106 }
107 return x
108}