Module ndtr 

The standard normal CDF ndtr, its logarithm log_ndtr, and logerfc = ln erfc.

Motivation

ndtr(x) = erfc(-x/sqrt 2)/2 is one line, and is here so the forward CDF exists beside its inverse probit. The two log forms are not re-expressions: ndtr underflows to zero near x = -38.6 (binary64) and -14.4 (binary32), and erfc near 27 / 9.3, so ln(ndtr(x)) and ln(erfc(x)) return -inf exactly where a log-likelihood, a censored-data model or a Bayesian-optimization acquisition function needs them most. Both logs are perfectly ordinary numbers there (log_ndtr(-100) = -5004.6). These kernels carry them.

Algorithm

Three arms per function, chosen by where each spelling is accurate. For log_ndtr, with u = |x|/sqrt 2:

x >  0:       ln_1p(-erfc(u)/2)           ndtr -> 1, the complement is the small side
-5.7 < x <= 0: ln(erfc(u)/2)               bit-identical to ln(ndtr(x))
x <= -5.7:    ln(erfcx(u)/2) - u^2        the tail; erfcx has no underflow

The moderate region runs on erfc deliberately. erfcx is a Weideman rational whose term count follows the policy, and at the default tier it is 4.2e-10 relative. The log turns relative error into absolute, so ln(erfcx(u)/2) at x = -1 would be 1e6 ulp off the ln(ndtr(x)) a caller could write by hand. erfc is a rational-times- exp fit that holds a few ulp at every tier, and its one weakness (the x^2 under the exp amplifies the argument’s rounding by x^2) is exactly what the log absorbs: an error of c x^2 epsilon relative to erfc is c x^2 epsilon absolute in the log, against a result of -x^2/2.

In the tail the same absorption is what makes erfcx affordable: its error goes into the log as an absolute delta, against a result dominated by -u^2, so the tail arm runs erfcx two rungs above the caller’s tier (LogTailPolicy, N = 40 from the default tier) and starts at u = 4, where 8.7e-16 / 16 is a quarter ulp. The threshold is where erfc still has a hundred orders of magnitude of headroom in either format. The arms share one ln: the argument is lane-selected between erfc/2 and erfcx/2, and only the tail subtracts u^2.

u^2 is formed as (|x|/2)|x| rather than x*x/2, so the intermediate does not overflow before the result does: log_ndtr(-1.3e154) is a representable -8.45e307.

logerfc is the same shape with the tail on the right (ln(erfcx(x)) - x^2 above x = 4, ln(erfc(x)) between 1/2 and 4) and a bounded left side. There erfc(-|x|) = 1 + erf(|x|), so the arm is ln_1p(erf(|x|)), which is 2|x|/sqrt(pi) near zero and must not be formed from erfc. The same cancellation sits in erfc(x) = 1 - erf(x) for small positive x (it rounds 1.128e-8 to 1.1e-8 at x = 1e-8), so below x = 1/2 the right side is ln_1p(-erf(x)) as well. A python model of the seam puts both spellings at 1-2 ulp on either side of it.

Each transcendental is evaluated only when some lane’s arm needs it, or unconditionally under avoid_branching. A packet of one sign in the moderate region, the common case, pays one erfc and one log.

Functions

inv_log_ndtr_impl

The inverse of log_ndtr_impl: the x with ln ndtr(x) = y, for y <= 0.

log_ndtr_impl

ln(ndtr(x)), finite for every finite x.

log_ndtr_with_deriv_impl

(ln ndtr(x), phi(x)/ndtr(x)): the value and its derivative, the inverse Mills ratio.

logerfc_impl

ln(erfc(x)), finite for every finite x.

ndtr_impl

ndtr(x) = erfc(-x/sqrt 2)/2, the standard normal CDF.

residual_tolerance

The residual tolerance the Newton inverses stop at: a tier-dependent number of ulps of scale, which the caller sets to the size of the quantity the residual is measured in.

Type Aliases

LogTailPolicy

The policy the tail arms evaluate erfcx under: two precision rungs above the caller’s, so the default tier takes the full N = 40 Weideman table.