Widen the Cody branch past the erfc^2 underflow (Cody 1969)

The lcdf value branch used log(erfc(-s)) for scaled_diff in (-20, 0]. In
forward-over-reverse mode (fvar<fvar<var>>) that produces a NaN third
derivative: fvar's log rule forms x.d_ / x.val_ with x.val_ = erfc(|s|), and
once erfc(|s|)^2 leaves the normal range the quotient chain evaluates
0 * finite / 0.

Measured, not predicted:
  erfc(|s|)^2 leaves the normal range at |s| = 18.7268381191
  NaN window bisects to s in (-20.0000000000, -19.2103473480)
    i.e. y in (-28.2842712475, -27.1675337575) at mu=0, sigma=1

Reachable ONLY in fvar<fvar<var>>, i.e. expect_ad's grad_hessian(); fvar<var>
and fvar<fvar<double>> are both clean, and nothing shows up as an inf in the
forward value slot. The existing y = -20*sqrt(2) test point sits exactly on
the lower edge and takes the Cody branch, which is why this was never caught.

Move the boundary to -12, handing s <= -12 to the Cody (1969) rational
approximation that already sits below it. That clears the underflow window
with 6.7 units of margin. Cody's relative error over the newly handed-over
range s in [-20, -12] measures at worst 9.593e-17 against a 60-digit mpmath
reference -- machine precision, and indistinguishable from the erfc branch it
replaces. The dense R-generated value tables in
prim/prob/{normal,std_normal}_cdf_log_test.cpp cover y in [-37.5, 10], which
spans the re-branched region, and no expectation moves.

P1 and P2 tests go green; P3 remains red.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: SteveBronder

stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp

@@ -29,7 +29,7 @@ static constexpr const char* std_normal_lcdf_device_function
               if (isnan(lcdf_n)) {
                 lcdf_n = 0;
               }
-            } else if (scaled_y > -20.0) {
+            } else if (scaled_y > -12.0) {
               // CDF(x) = 1/2 - 1/2 erf(-x) = 1/2 erfc(-x)
               lcdf_n = log(erfc(-scaled_y)) - M_LN2;
             } else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) {

stan/math/opencl/prim/normal_lcdf.hpp

@@ -30,7 +30,7 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY(
       if (isnan(normal_lcdf_n)) {
         normal_lcdf_n = 0;
       }
-    } else if (normal_lcdf_scaled_diff > -20.0) {
+    } else if (normal_lcdf_scaled_diff > -12.0) {
       // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x)
       normal_lcdf_n = log(erfc(-normal_lcdf_scaled_diff)) - M_LN2;
     } else if (10.0 * log(fabs(normal_lcdf_scaled_diff)) < log(DBL_MAX)) {

stan/math/prim/prob/normal_lcdf.hpp

@@ -91,8 +91,11 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lcdf(const T_y& y,
       if (!is_not_nan(cdf_log)) {
         cdf_log = 0;
       }
-    } else if (scaled_diff > -20.0) {
+    } else if (scaled_diff > -12.0) {
       // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x)
+      // bounded away from -18.727, where erfc(-x)^2 leaves the normal
+      // range and the autodiff quotient rule forms 0/0; Cody (1969)
+      // below is accurate to 9.6e-17 over the range handed to it
       cdf_log += log(erfc(-scaled_diff)) + LOG_HALF;
     } else if (10.0 * log(fabs(scaled_diff))
                < log(std::numeric_limits<T_partials_return>::max())) {

stan/math/prim/prob/std_normal_lcdf.hpp

@@ -71,7 +71,10 @@ inline return_type_t<T_y> std_normal_lcdf(const T_y& y) {
       if (!is_not_nan(lcdf)) {
         lcdf = 0;
       }
-    } else if (scaled_y > -20.0) {
+    } else if (scaled_y > -12.0) {
+      // bounded away from -18.727, where erfc(-x)^2 leaves the normal
+      // range and the autodiff quotient rule forms 0/0; Cody (1969)
+      // below is accurate to 9.6e-17 over the range handed to it
       // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x)
       lcdf += log(erfc(-scaled_y)) + LOG_HALF;
     } else if (10.0 * log(fabs(scaled_y))