Scaled complementary error function erfcx(x) = exp(x²) · erfc(x),
specialized to x ≥ 0. Returns 1.0 for x ≤ 0 and 0.0 for
x = +∞. For 0 < x < 26 uses the direct exp(x²)·erfc(x) form;
beyond that the (otherwise overflowing) exp(x²) is replaced by a
4-term asymptotic expansion (1/(x√π))·(1 − 1/(2x²) + 3/(4x⁴) − …),
keeping relative accuracy near machine epsilon. The non-negative
restriction lets the caller skip the reflection identity.
Numerically stable ln Φ(x) for the standard normal CDF. For x ≥ 0
computes ln(Φ(x)) directly with a small floor against underflow; for
x < 0 rewrites
ln Φ(x) = −u² + ln(½·erfcx(u)), u = −x/√2,
which preserves digits all the way into the deep left tail (no
ln(0)). Returns ±∞ and NaN at the corresponding inputs.
Numerically stable ln(1 − Φ(x)) = ln Φ(−x) for the standard normal
survival function. Delegates to normal_logcdf(-x) so the deep-right
tail benefits from the same erfcx-based representation.
Numerically stable signed log-sum-exp. Given pairs
(log|aⱼ|, sign(aⱼ)) (with signs[j] ∈ {−1, 0, +1}), returns
(log|S|, sign(S)) for S = Σⱼ signs[j]·exp(log_mags[j]). Positive
and negative magnitudes are reduced separately with the standard
log-sum-exp trick (subtract the max, sum, log, add back); the two
partial sums are then combined via log(|p − n|) = max(log p, log n) + log1mexp(|log p − log n|), preserving accuracy
even when p ≈ n (catastrophic cancellation regime). When all
signs are zero or all magnitudes are −∞, returns
(NEG_INFINITY, 0.0).
Joint evaluation of ln Φ(x) and the Mills-ratio analogue
φ(x) / Φ(x), signed for the symmetric branch. Used by the latent
probit families where the inverse-link gradient needs the ratio and
the likelihood needs the log-CDF on the same x; computing both in
one call shares the erfcx evaluation that dominates the cost in the
deep tail.