Use the asymptotic Mills ratio in the far lower tail (DLMF 7.12.1)

The scaled_diff < -29 branch took the Abramowitz & Stegun negative-tail
approximation and added a cubic residual fit whose leading term is
0.0015065154280332 * x2. That term grows quadratically in scaled_diff.

DLMF 7.12.1 gives, for y = -scaled_diff > 0,

  erfc(y) ~ (exp(-y^2) / (sqrt(pi) y)) * sum_m (-1)^m (1/2)_m / y^(2m)

and since dncdf_log = (2/sqrt(pi)) exp(-y^2) / erfc(y), this is

  dncdf_log ~ 2y / (1 - 1/(2y^2) + 3/(4y^4) - 15/(8y^6) + ...)

so dncdf_log grows LINEARLY, as -2*scaled_diff. A quadratic residual fit
cannot track a linear asymptote; it must eventually dominate. Measured
relative error of the old branch against a 60-digit mpmath reference:

  scaled_diff   -29       -40       -75      -100     -1000
  old         9.5e-06   1.9e-04   1.0e-02   2.3e-02   6.8e-01
  new         1.3e-11   1.0e-12   6.5e-15   7.0e-16   7.8e-18

At scaled_diff = -1000 the old branch is 68% wrong. In y terms (mu=0,
sigma=1) the gradient passes 1e-4 relative error at y = -54.46 and 1e-2 at
y = -106.06. This is a plain first-order var gradient, so it corrupts HMC
transitions, not just Hessians -- and PR #3363's reflection newly exposes it
on the lccdf side at large positive y.

Replacing the branch with the truncated asymptotic is also exponential-free:
no erf, no exp, nothing that can overflow. DLMF states the remainder is
bounded by the first neglected term for real argument, which at
scaled_diff = -29 is 6.5625/29^8 = 1.3e-11, matching the measured error.

The seam at scaled_diff = -29 moves from 3.3e-06 to 6.2e-06 relative, still
well inside the existing inter-branch jumps (~2.7e-5 at scaled_diff = 2.9).

Also corrects the sqrt(2) scaling in the new rev sweep assertions: dncdf_log
is d/d(scaled_diff), while the reported partial is dncdf_log/(sigma*sqrt(2)).
The discrete mpmath-referenced cases were unaffected. Verified the corrected
sweep is still red without this fix and green with it.

All P1, P2 and P3 tests now pass across prim, rev, fwd and mix.

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

Author: SteveBronder

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

@@ -121,6 +121,14 @@ static constexpr const char* std_normal_lcdf_device_function
               dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
                        + 0.04700850676 * t2 * t - 0.03478651979 * t4
                        - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
+            } else if (scaled_y < -29.0) {
+              // asymptotic Mills ratio, DLMF 7.12.1: grows linearly as
+              // -2*scaled_y, so no quadratic residual fit can track it
+              const double inv_x2 = 1.0 / x2;
+              dnlcdf
+                  = -2.0 * scaled_y
+                    / (1.0
+                       + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875)));
             } else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) {
               t = 1.0 / (1.0 - 0.3275911 * scaled_y);
               t2 = t * t;
@@ -129,10 +137,7 @@ static constexpr const char* std_normal_lcdf_device_function
                   = M_2_SQRTPI
                     / (0.254829592 * t - 0.284496736 * t2 + 1.421413741 * t2 * t
                        - 1.453152027 * t4 + 1.061405429 * t4 * t);
-              if (scaled_y < -29.0) {
-                dnlcdf += 0.0015065154280332 * x2
-                          - 0.3993154819705530 * scaled_y - 4.2919418242931700;
-              } else if (scaled_y < -17.0) {
+              if (scaled_y < -17.0) {
                 dnlcdf += 0.0001263257217272 * x2 * scaled_y
                           + 0.0123586859488623 * x2
                           - 0.0860505264736028 * scaled_y - 1.252783383752970;

stan/math/opencl/prim/normal_lcdf.hpp

@@ -60,6 +60,7 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY(
 const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
     double normal_ldncdf = 0.0; double t = 0.0; double t2 = 0.0;
     double t4 = 0.0; double normal_lcdf_exp_m_x2 = 0.0;
+    double normal_lcdf_inv_x2 = 0.0;
 
     // calculate using piecewise function
     // (due to instability / inaccuracy in the various approximations)
@@ -120,6 +121,17 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
       normal_ldncdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
                       + 0.04700850676 * t2 * t - 0.03478651979 * t4
                       - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
+    } else if (normal_lcdf_deriv_scaled_diff < -29.0) {
+      // asymptotic Mills ratio, DLMF 7.12.1: grows linearly as -2*scaled_diff,
+      // so no quadratic residual fit can track it
+      normal_lcdf_inv_x2 = 1.0 / x2;
+      normal_ldncdf
+          = -2.0 * normal_lcdf_deriv_scaled_diff
+            / (1.0
+               + normal_lcdf_inv_x2
+                     * (-0.5
+                        + normal_lcdf_inv_x2
+                              * (0.75 + normal_lcdf_inv_x2 * -1.875)));
     } else if (10.0 * log(fabs(normal_lcdf_deriv_scaled_diff)) < log(DBL_MAX)) {
       // approximation derived from Abramowitz and Stegun (1964) 7.1.26
       // use fact that erf(x)=-erf(-x)
@@ -134,11 +146,7 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
                - 1.453152027 * t4 + 1.061405429 * t4 * t);
       // check if we need to add a correction term
       // (from cubic fit of residuals)
-      if (normal_lcdf_deriv_scaled_diff < -29.0) {
-        normal_ldncdf += 0.0015065154280332 * x2
-                         - 0.3993154819705530 * normal_lcdf_deriv_scaled_diff
-                         - 4.2919418242931700;
-      } else if (normal_lcdf_deriv_scaled_diff < -17.0) {
+      if (normal_lcdf_deriv_scaled_diff < -17.0) {
         normal_ldncdf += 0.0001263257217272 * x2 * normal_lcdf_deriv_scaled_diff
                          + 0.0123586859488623 * x2
                          - 0.0860505264736028 * normal_lcdf_deriv_scaled_diff

stan/math/prim/prob/normal_lcdf.hpp

@@ -191,6 +191,13 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lcdf(const T_y& y,
         dncdf_log = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
                     + 0.04700850676 * t2 * t - 0.03478651979 * t4
                     - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
+      } else if (scaled_diff < -29.0) {
+        // asymptotic Mills ratio, DLMF 7.12.1: dncdf_log grows linearly as
+        // -2*scaled_diff, so no quadratic residual fit can track it
+        const T_partials_return inv_x2 = 1.0 / x2;
+        dncdf_log
+            = -2.0 * scaled_diff
+              / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875)));
       } else if (10.0 * log(fabs(scaled_diff))
                  < log(std::numeric_limits<T_partials_return>::max())) {
         // approximation derived from Abramowitz and Stegun (1964) 7.1.26
@@ -207,10 +214,7 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lcdf(const T_y& y,
                     - 1.453152027 * t4 + 1.061405429 * t4 * t));
         // check if we need to add a correction term
         // (from cubic fit of residuals)
-        if (scaled_diff < -29.0) {
-          dncdf_log += 0.0015065154280332 * x2
-                       - 0.3993154819705530 * scaled_diff - 4.2919418242931700;
-        } else if (scaled_diff < -17.0) {
+        if (scaled_diff < -17.0) {
           dncdf_log += 0.0001263257217272 * x2 * scaled_diff
                        + 0.0123586859488623 * x2
                        - 0.0860505264736028 * scaled_diff - 1.252783383752970;

stan/math/prim/prob/std_normal_lcdf.hpp

@@ -169,6 +169,12 @@ inline return_type_t<T_y> std_normal_lcdf(const T_y& y) {
         dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
                  + 0.04700850676 * t2 * t - 0.03478651979 * t4
                  - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
+      } else if (scaled_y < -29.0) {
+        // asymptotic Mills ratio, DLMF 7.12.1: dnlcdf grows linearly as
+        // -2*scaled_y, so no quadratic residual fit can track it
+        const T_partials_return inv_x2 = 1.0 / x2;
+        dnlcdf = -2.0 * scaled_y
+                 / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875)));
       } else if (10.0 * log(fabs(scaled_y))
                  < log(std::numeric_limits<T_partials_return>::max())) {
         // approximation derived from Abramowitz and Stegun (1964) 7.1.26
@@ -183,10 +189,7 @@ inline return_type_t<T_y> std_normal_lcdf(const T_y& y) {
                     - 1.453152027 * t4 + 1.061405429 * t4 * t);
         // check if we need to add a correction term
         // (from cubic fit of residuals)
-        if (scaled_y < -29.0) {
-          dnlcdf += 0.0015065154280332 * x2 - 0.3993154819705530 * scaled_y
-                    - 4.2919418242931700;
-        } else if (scaled_y < -17.0) {
+        if (scaled_y < -17.0) {
           dnlcdf += 0.0001263257217272 * x2 * scaled_y + 0.0123586859488623 * x2
                     - 0.0860505264736028 * scaled_y - 1.252783383752970;
         } else if (scaled_y < -7.0) {

test/unit/math/rev/prob/normal_ccdf_log_test.cpp

@@ -49,7 +49,8 @@ TEST(ProbDistributionsNormalLccdf, var_gradient_matches_mills_ratio) {
     const double mills
         = -2.0 * s
           / (1.0 + inv * (-0.5 + inv * (0.75 + inv * (-1.875 + inv * 6.5625))));
-    EXPECT_LT(std::fabs(yv.adj() / -mills - 1.0), 1e-9)
+    // dncdf_log is d/ds; the partial is dncdf_log / (sigma*sqrt(2))
+    EXPECT_LT(std::fabs(yv.adj() / (-mills / std::sqrt(2.0)) - 1.0), 1e-9)
         << "gradient drifts from the asymptotic Mills ratio at y = " << y;
     stan::math::set_zero_all_adjoints();
   }

test/unit/math/rev/prob/normal_cdf_log_test.cpp

@@ -60,7 +60,8 @@ TEST(ProbDistributionsNormalLcdf, var_gradient_matches_mills_ratio) {
     const double mills
         = -2.0 * s
           / (1.0 + inv * (-0.5 + inv * (0.75 + inv * (-1.875 + inv * 6.5625))));
-    EXPECT_LT(std::fabs(yv.adj() / mills - 1.0), 1e-9)
+    // dncdf_log is d/ds; the partial is dncdf_log / (sigma*sqrt(2))
+    EXPECT_LT(std::fabs(yv.adj() / (mills / std::sqrt(2.0)) - 1.0), 1e-9)
         << "gradient drifts from the asymptotic Mills ratio at y = " << y;
     stan::math::set_zero_all_adjoints();
   }

test/unit/math/rev/prob/std_normal_ccdf_log_test.cpp

@@ -48,7 +48,8 @@ TEST(ProbDistributionsStdNormalLccdf, var_gradient_matches_mills_ratio) {
     const double mills
         = -2.0 * s
           / (1.0 + inv * (-0.5 + inv * (0.75 + inv * (-1.875 + inv * 6.5625))));
-    EXPECT_LT(std::fabs(yv.adj() / -mills - 1.0), 1e-9)
+    // dncdf_log is d/ds; the partial is dncdf_log / (sigma*sqrt(2))
+    EXPECT_LT(std::fabs(yv.adj() / (-mills / std::sqrt(2.0)) - 1.0), 1e-9)
         << "gradient drifts from the asymptotic Mills ratio at y = " << y;
     stan::math::set_zero_all_adjoints();
   }

test/unit/math/rev/prob/std_normal_cdf_log_test.cpp

@@ -50,7 +50,8 @@ TEST(ProbDistributionsStdNormalLcdf, var_gradient_matches_mills_ratio) {
     const double mills
         = -2.0 * s
           / (1.0 + inv * (-0.5 + inv * (0.75 + inv * (-1.875 + inv * 6.5625))));
-    EXPECT_LT(std::fabs(yv.adj() / mills - 1.0), 1e-9)
+    // dncdf_log is d/ds; the partial is dncdf_log / (sigma*sqrt(2))
+    EXPECT_LT(std::fabs(yv.adj() / (mills / std::sqrt(2.0)) - 1.0), 1e-9)
         << "gradient drifts from the asymptotic Mills ratio at y = " << y;
     stan::math::set_zero_all_adjoints();
   }