Module erfcx 

erfcx(x) = e^{x^2} erfc(x), the scaled complementary error function.

Motivation

erfc underflows to zero at x ~ 27 in binary64 and x ~ 9 in binary32, where the true value is e^{-x^2}/(x sqrt(pi)), nonzero and merely unrepresentable. Every Gaussian tail, importance weight and log-likelihood past that point silently becomes zero. erfcx removes the exponential and is O(1/x), so it stays representable for every finite argument and carries full relative accuracy the whole way.

Algorithm

The Faddeeva function restricted to the imaginary axis: w(ix) = erfcx(x) exactly. Weideman’s approximation (see [crate::tables::weideman]) is

Z = (L + iz)/(L - iz),   w(z) = 1/(sqrt(pi)(L - iz)) + 2 P(Z)/(L - iz)^2

with P real. Substituting z = ix for real x makes L - iz = L + x and Z = (L - x)/(L + x), both real: every complex operation in the method disappears and what is left is one reciprocal and one real Horner. There are no transcendentals at all on the non-negative side, which makes this cheaper than the erfc it complements.

The domain is well conditioned throughout. L + x >= L > 0 for every finite x >= 0, so the reciprocal needs no guard, and Z runs monotonically over (-1, 1], so the Horner stays inside the unit disc the coefficients were fitted on. Only the infinities fall outside that, Z there being (L - inf) * 0, and they are named explicitly under check_overflow.

Measured against mpmath at 50 digits with the N = 40 table, the worst relative error over x from 0 to 1e15 is 1.22 ulp.

Negative arguments

erfcx(-x) = 2 e^{x^2} - erfcx(x), which genuinely overflows for x below about -26.6 (binary64). erfcx grows like e^{x^2} to the left, so the infinity is the correct answer rather than a failure. This is the only branch, and it is the only place an exp appears.

Functions

erfcx_internal

erfcx_with, with N and the table chosen by the precision policy.

erfcx_with

erfcx by the N-term Weideman approximation on the imaginary axis.