Add `erfcx` (scaled complementary error function) as a Stan Math primitive

[avehtari]

I think it would be useful to add erfcx(x) = exp(x²)·erfc(x) to prim/rev/fwd and the OpenCL kernel generator, as a first-class function. I have PR ready, but would like ask first whether others agree that this would be worth doing.

This would be the generalisation of #2803, which already identified erfcx as the fix for exp_mod_normal overflow. The same primitive would fixe or improve at least five other distribution functions (see below).

Boost has erf and erfc, but doesn't have erfcx. SciPy, Julia, libcerf, and Fortran 2008 has this so it's known to be useful. Good thing is that Stan Math already ships the algorithm and needed coefficients in normal_lcdf, but the bad thing is that currently nothing else can reach them.

#2803 reports exp_mod_normal_lpdf(1 | 0, 1, 40) returning -Inf where the correct value is -1.394277. The cause is that the density is formed as exp(a²/2 − az) · erfc((a−z)/√2) with a = λσ: above a ≈ 37.6 the first factor overflows while the second underflows, so the product is inf · 0 = NaN, or -Inf after the log. The scaled form never separates the factors.

The same cancellation appears wherever a normal tail is formed. Two identities cover most cases, both free of cancellation and overflow:

• log Φ(x) = LOG_HALF + log(erfcx(−x·INV_SQRT_TWO)) − x²/2
• φ(x)/Φ(x) = SQRT_TWO_OVER_SQRT_PI / erfcx(−x·INV_SQRT_TWO) (inverse Mills ratio)

skew_normal_lpdf returns -Inf with infinite gradients in the tail. It computes log(erfc(−α·z/√2)) and √(2/π)·exp(−(αz)²/2)/erfc(−αz/√2), which are exactly log(2Φ(αz)) and the inverse Mills ratio. With α = 1:

skew_normal_lpdf(-40, 0, 1, 1)   current: -Inf, gradients Inf
                                 correct: -1604.8342333664

std::erfc(28.284) is exactly 0, so the log is -Inf and exp(-sq - log_erfc) is Inf. Both lines become one call to the identities above.

std_normal_lcdf / normal_lcdf gradients are ~1e-5 accurate. The gradient is the inverse Mills ratio, currently computed by roughly 120 lines of A&S 7.1.26 plus five Taylor expansions plus four residual-fit corrections. normal_lcdf.hpp's own documentation tabulates the worst in-range relative error of those branches as 7.61e-06 to 6.09e-05. The erfcx identity gives the same quantity to machine precision and would delete the table of magic constants. (normal_lcdf.hpp already cites SciPy's log_ndtr as needing "no rational approximation or Taylor patches at all below x = -1, because erfcx never forms exp(+t²)".)

Others in the same shape: normal_lccdf, std_normal_lccdf, lognormal_{cdf,lcdf,lccdf} (normal on the log scale), ordered_probit_lpmf (differences of normal CDFs), and the erfc half of skew_normal_{cdf,lcdf,lccdf}.

Prior art

software	approach	licence
Faddeeva (S. G. Johnson)	continued fraction for x > 50; Chebyshev under y = 4/(4+x), 100 subintervals with lookup table	MIT
SciPy	Faddeeva (still cited in v1.18)	BSD-3 + vendored MIT
Julia	Faddeeva, vendored in JuliaMath/openspecfun	MIT
libcerf	piecewise Chebyshev; rewritten in 2.5, guarantees 2.34e-16 relative for x < 12	—
Fortran 2008	erfc_scaled intrinsic, typically Cody (1969)	—
Boost	no erfcx or erfc_scaled (checked 1.87, the vendored version)	—

Notably, no serious implementation forms exp(x²)·erfc(x) as its primary method.

Implementation

#2803 includes a candidate implementation reproduced from a StackOverflow answer. It is likely excellent, but StackOverflow content is CC BY-SA, which is generally treated as incompatible with Stan Math's BSD-3, so vendoring it would need an explicit grant from its author. Worth pursuing separately.

Faddeeva is MIT and therefore usable, but vendoring it means carrying a source file and a 100-entry Chebyshev table.

The proposal here needs neither, because Cody (1969) is already in the repo. normal_lcdf.hpp uses those coefficients for its tail value, and rearranging its own expression yields erfcx directly, with no exp and no erfc call at all:

erfcx(w) = (INV_SQRT_PI + (p/q)/w²) / w

So the implementation is two branches:

• x ≥ 4 — Cody rational, Horner form in u = 1/x². The crossover matches the one normal_lcdf already documents and cross-checks against R's pnorm.
• x < 4 — exp(x*x)*erfc(x), with x*x split so the argument handed to exp is exact. Without the split, the rounding of x*x is amplified by exp into a relative error of about x²·eps.

Measured against a long double reference, 40 000 points per range:

range	naive exp(x*x)*erfc(x)	proposed
[−26.6, −10]	512.4 ulp	2.73
[−4, 0]	17.5	3.81
[0, 4]	11.3	4.99
[4, 9]	65.9	2.00
[9, 20]	257.2	2.06
[20, 50]	NaN (inf·0)	2.49
[50, 200]	—	5.42

Worst case 5.4 ulp overall, 2.0–2.5 ulp in the upper tail where it matters. Speed at x ∈ [5,30]: 5.4 ns/call versus 21.2 for the naive product.

Derivative, which reuses the value so it is exact wherever the value is:

d/dx erfcx(x) = 2·x·erfcx(x) − 2/√π

Edge cases: erfcx(0) = 1, erfcx(+inf) = 0, erfcx(-inf) = +inf, NaN preserved.

Follow-ups

• One PR per distribution with bigger improvement in accuracy: exp_mod_normal, skew_normal, lognormal, ordered_probit
• One PR per disttribution with smaller improvement in accuracy + dropping second copy of the same coefficients: normal, std_normal

[avehtari]

ping @andrjohns, @WardBrian, @spinkney

[WardBrian]

Seems reasonable. I’d be happy to re-use the existing code we have or vendor the MIT one if it’s just one file and has some advantage in speed or accuracy

[avehtari]

It seems the Stan's version inside normal_lcdf is accurate and fast, and already battle tested in Stan