[avehtari]

Closes #352.

Summary

The functions reference says in six places that Phi_approx, log(Phi_approx(...)) and log1m(Phi_approx(...)) are more robust in the tails than normal_cdf, normal_lcdf, normal_lccdf and their std_normal_ versions. It also says that the log functions underflow or overflow below -37.5 and above 8.25. Both statements are wrong. This PR replaces them and adds a tail warning to Phi_approx.

Evidence

Relative error against mpmath (60 digits) on Stan Math develop of 2026-10-22:

function	z range	relative error
normal_lcdf	-1e5 to 37	at most 1e-13
normal_lccdf	-37.5 to 40	at most 1e-13
Phi, normal_cdf	-37.5 to 8	at most 4e-11
log(Phi_approx(z))	-2 to 2	at most 3e-3
log(Phi_approx(z))	-21.5 to -5	0.1 to 2
log(Phi_approx(z))	5 to 7.1	0.8 to 1
log(Phi_approx(z))	below -21.6	-inf (the true value is finite)
log(Phi_approx(z))	above 7.1	1 (returns 0)

With stan-dev/math#3406 the first three rows are at most 4e-16. Phi_approx is inv_logit(0.07056 x^3 + 1.5976 x), and its log behaves like 0.07056 x^3 in the lower tail instead of -x^2/2. For example, Phi_approx(-10) is 2.6e-38 and Phi(-10) is 7.6e-24.

The statements on the log functions have been out of date since stan-dev/math#1411 (2019) for normal_lcdf and since stan-dev/math#3363 (2026) for normal_lccdf.

[bob-carpenter]

Thanks, Aki---this is super helpful