From a335aa2eddf6065fc09b9c1bb05f6712018d6841 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 11:25:51 +0200
Subject: [PATCH 01/10] Tests: Phi_approx derivative in the upper tail

---
 test/unit/math/mix/fun/Phi_approx_test.cpp    | 74 +++++++++++++++++++
 test/unit/math/opencl/rev/Phi_approx_test.cpp | 29 ++++++++
 2 files changed, 103 insertions(+)

diff --git a/test/unit/math/mix/fun/Phi_approx_test.cpp b/test/unit/math/mix/fun/Phi_approx_test.cpp
index b89df1935b7..1bb709bc9da 100644
--- a/test/unit/math/mix/fun/Phi_approx_test.cpp
+++ b/test/unit/math/mix/fun/Phi_approx_test.cpp
@@ -1,4 +1,7 @@
 #include <test/unit/math/test_ad.hpp>
+#include <cmath>
+#include <limits>
+#include <vector>
 
 TEST(mathMixMatFun, PhiApprox) {
   auto f = [](const auto& x1) { return stan::math::Phi_approx(x1); };
@@ -23,3 +26,74 @@ TEST(mathMixMatFun, PhiApprox_varmat) {
   }
   expect_ad_vector_matvar(f, A);
 }
+
+namespace Phi_approx_test {
+// d/dx inv_logit(0.07056 x^3 + 1.5976 x) by mpmath at 60 digits; at x = 40
+// and -40 the true value, about 4e-1987, underflows to 0
+inline std::vector<double> args() {
+  return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
+}
+inline std::vector<double> derivs() {
+  return {0.0,
+          5.9589784505245725e-37,
+          8.7097694481792379e-21,
+          5.1340418352583332e-15,
+          0.24152728992165639,
+          0.3994,
+          0.24152728992165639,
+          5.1340418352583332e-15,
+          8.7097694481792379e-21,
+          5.9589784505245725e-37,
+          0.0};
+}
+// Relative check, absolute below the smallest normal double.
+inline void expect_rel(double expected, double actual, double x,
+                       double tol = 1e-13) {
+  EXPECT_NEAR(
+      expected, actual,
+      tol * std::fmax(std::fabs(expected), std::numeric_limits<double>::min()))
+      << "x = " << x;
+}
+}  // namespace Phi_approx_test
+
+TEST(mathMixMatFun, PhiApproxScalarDerivativeTails) {
+  using stan::math::var;
+  std::vector<double> xs = Phi_approx_test::args();
+  std::vector<double> ds = Phi_approx_test::derivs();
+  for (size_t i = 0; i < xs.size(); ++i) {
+    var a = xs[i];
+    stan::math::Phi_approx(a).grad();
+    Phi_approx_test::expect_rel(ds[i], a.adj(), xs[i]);
+    stan::math::recover_memory();
+  }
+}
+
+TEST(mathMixMatFun, PhiApproxDerivativeTails) {
+  using stan::math::Phi_approx;
+  using stan::math::sum;
+  using stan::math::var;
+  std::vector<double> xs = Phi_approx_test::args();
+  std::vector<double> ds = Phi_approx_test::derivs();
+  Eigen::VectorXd x = stan::math::to_vector(xs);
+
+  stan::math::var_value<Eigen::VectorXd> xv(x);
+  sum(Phi_approx(xv)).grad();
+  for (int i = 0; i < x.size(); ++i) {
+    Phi_approx_test::expect_rel(ds[i], xv.adj()(i), xs[i]);
+  }
+  stan::math::recover_memory();
+
+  Eigen::Matrix<var, Eigen::Dynamic, 1> xm = x;
+  sum(Phi_approx(xm)).grad();
+  for (int i = 0; i < x.size(); ++i) {
+    Phi_approx_test::expect_rel(ds[i], xm(i).adj(), xs[i]);
+  }
+  stan::math::recover_memory();
+
+  std::vector<var> xs_v(xs.begin(), xs.end());
+  sum(Phi_approx(xs_v)).grad();
+  for (size_t i = 0; i < xs.size(); ++i) {
+    Phi_approx_test::expect_rel(ds[i], xs_v[i].adj(), xs[i]);
+  }
+  stan::math::recover_memory();
+}
diff --git a/test/unit/math/opencl/rev/Phi_approx_test.cpp b/test/unit/math/opencl/rev/Phi_approx_test.cpp
index 51c2bc14a49..d3e06c16fd5 100644
--- a/test/unit/math/opencl/rev/Phi_approx_test.cpp
+++ b/test/unit/math/opencl/rev/Phi_approx_test.cpp
@@ -13,6 +13,35 @@ TEST(OpenCL_Phi_approx, prim_rev_values_small) {
   stan::math::test::compare_cpu_opencl_prim_rev(Phi_approx_functor, a);
 }
 
+/**
+ * Device derivative in the tails against fixed references, which
+ * `compare_cpu_opencl_prim_rev` cannot check: its relative tolerance floors
+ * the denominator at 1. References by mpmath at 60 digits.
+ */
+TEST(OpenCL_Phi_approx, rev_tail_derivative_against_references) {
+  const int N = 8;
+  Eigen::VectorXd x(N);
+  x << -40, -10, -7, 1, 7, 8, 10, 40;
+  // at x = -40 and 40 the true value, about 4e-1987, underflows to 0
+  Eigen::VectorXd expected(N);
+  expected << 0.0, 5.9589784505245725e-37, 5.1340418352583332e-15,
+      0.24152728992165639, 5.1340418352583332e-15, 8.7097694481792379e-21,
+      5.9589784505245725e-37, 0.0;
+
+  stan::math::var_value<stan::math::matrix_cl<double>> x_cl(
+      stan::math::to_matrix_cl(x));
+  stan::math::sum(stan::math::Phi_approx(x_cl)).grad();
+
+  const Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj());
+  for (int i = 0; i < N; ++i) {
+    EXPECT_NEAR(
+        expected[i], adj[i],
+        1e-13 * std::fmax(expected[i], std::numeric_limits<double>::min()))
+        << "x = " << x[i];
+  }
+  stan::math::recover_memory();
+}
+
 TEST(OpenCL_Phi_approx, prim_rev_size_0) {
   int N = 0;
 

From 198d939115c4dbe6393a2589225d6667f050fec6 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 11:25:51 +0200
Subject: [PATCH 02/10] Fix Phi_approx reverse-mode derivative cancellation in
 the upper tail

---
 stan/math/opencl/rev/Phi_approx.hpp |  7 +++++--
 stan/math/rev/fun/Phi_approx.hpp    | 26 +++++++++++++-------------
 2 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/stan/math/opencl/rev/Phi_approx.hpp b/stan/math/opencl/rev/Phi_approx.hpp
index 4653b7e542f..708d435b9ed 100644
--- a/stan/math/opencl/rev/Phi_approx.hpp
+++ b/stan/math/opencl/rev/Phi_approx.hpp
@@ -20,9 +20,12 @@ template <typename T,
 inline var_value<matrix_cl<double>> Phi_approx(const var_value<T>& A) {
   return make_callback_var(
       Phi_approx(A.val()), [A](vari_value<matrix_cl<double>>& res) mutable {
+        // e / (1 + e)^2, not res * (1 - res), which cancels for u >> 0
+        auto x_sq = square(A.val());
+        auto e = exp(-fabs(elt_multiply(A.val(), 0.07056 * x_sq + 1.5976)));
         A.adj() += elt_multiply(
-            elt_multiply(elt_multiply(res.adj(), res.val()), 1 - res.val()),
-            3.0 * 0.07056 * square(A.val()) + 1.5976);
+            elt_multiply(res.adj(), elt_divide(e, square(1.0 + e))),
+            3.0 * 0.07056 * x_sq + 1.5976);
       });
 }
 
diff --git a/stan/math/rev/fun/Phi_approx.hpp b/stan/math/rev/fun/Phi_approx.hpp
index 0219c92853d..1d3fe29af20 100644
--- a/stan/math/rev/fun/Phi_approx.hpp
+++ b/stan/math/rev/fun/Phi_approx.hpp
@@ -3,8 +3,8 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
-#include <stan/math/rev/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/Phi_approx.hpp>
+#include <cmath>
 
 namespace stan {
 namespace math {
@@ -46,26 +46,26 @@ namespace math {
  */
 inline var Phi_approx(const var& a) {
   double av_squared = a.val() * a.val();
-  double f = inv_logit(0.07056 * a.val() * av_squared + 1.5976 * a.val());
-  double da = f * (1 - f) * (3.0 * 0.07056 * av_squared + 1.5976);
+  double u = 0.07056 * a.val() * av_squared + 1.5976 * a.val();
+  // e / (1 + e)^2, not f * (1 - f), which cancels for u >> 0
+  double e = std::exp(-std::fabs(u));
+  double f = (u >= 0 ? 1.0 : e) / (1.0 + e);
+  double da
+      = e / ((1.0 + e) * (1.0 + e)) * (3.0 * 0.07056 * av_squared + 1.5976);
   return make_callback_var(
       f, [a, da](auto& vi) mutable { a.adj() += vi.adj() * da; });
 }
 
 template <typename T, require_var_matrix_t<T>* = nullptr>
 inline auto Phi_approx(const T& a) {
+  const auto& x = a.val().array();
   arena_t<value_type_t<T>> f(a.rows(), a.cols());
   arena_t<value_type_t<T>> da(a.rows(), a.cols());
-  for (Eigen::Index j = 0; j < a.cols(); ++j) {
-    for (Eigen::Index i = 0; i < a.rows(); ++i) {
-      const auto a_val = a.val().coeff(i, j);
-      const auto av_squared = a_val * a_val;
-      f.coeffRef(i, j) = inv_logit(0.07056 * a_val * av_squared
-                                   + 1.5976 * a.val().coeff(i, j));
-      da.coeffRef(i, j) = f.coeff(i, j) * (1 - f.coeff(i, j))
-                          * (3.0 * 0.07056 * av_squared + 1.5976);
-    }
-  }
+  // e / (1 + e)^2, not f * (1 - f), which cancels for u >> 0
+  da.array() = (-(0.07056 * x * x.square() + 1.5976 * x).abs()).exp();
+  f.array() = (x >= 0).select(1.0, da.array()) / (1.0 + da.array());
+  da.array()
+      *= (3.0 * 0.07056 * x.square() + 1.5976) / (1.0 + da.array()).square();
   return make_callback_var(f, [a, da](auto& vi) mutable {
     a.adj().array() += vi.adj().array() * da.array();
   });

From c305f3ba733c26d25e878b69488d96ee78a0f564 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 11:25:53 +0200
Subject: [PATCH 03/10] Tests: inv_logit and log_inv_logit_diff derivatives at
 large arguments

---
 .../mix/constraint/lub_constrain_1_test.cpp   | 10 +++
 .../unit/math/mix/fun/Phi_approx_fwd_test.cpp | 66 ++++++++++++++
 .../unit/math/mix/fun/inv_logit_tail_refs.hpp | 86 +++++++++++++++++++
 test/unit/math/mix/fun/inv_logit_test.cpp     | 85 ++++++++++++++++++
 .../math/mix/fun/log1m_inv_logit_test.cpp     |  9 ++
 test/unit/math/mix/fun/log1p_exp_test.cpp     | 46 ++++++++++
 .../math/mix/fun/log_inv_logit_diff_test.cpp  | 65 ++++++++++++++
 test/unit/math/mix/fun/log_inv_logit_test.cpp |  9 ++
 test/unit/math/opencl/rev/inv_logit_test.cpp  | 34 ++++++++
 .../opencl/rev/log_inv_logit_diff_test.cpp    | 30 +++++++
 10 files changed, 440 insertions(+)
 create mode 100644 test/unit/math/mix/fun/Phi_approx_fwd_test.cpp
 create mode 100644 test/unit/math/mix/fun/inv_logit_tail_refs.hpp

diff --git a/test/unit/math/mix/constraint/lub_constrain_1_test.cpp b/test/unit/math/mix/constraint/lub_constrain_1_test.cpp
index d76175a441c..09ffa4b2e9c 100644
--- a/test/unit/math/mix/constraint/lub_constrain_1_test.cpp
+++ b/test/unit/math/mix/constraint/lub_constrain_1_test.cpp
@@ -1,5 +1,6 @@
 #include <test/unit/math/test_ad.hpp>
 #include <test/unit/math/mix/constraint/lub_constrain_helpers.hpp>
+#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
 
 TEST(mathMixMatFun, lub_constrain_scalars) {
   double x1 = 0.7;
@@ -24,3 +25,12 @@ TEST(mathMixMatFun, lub_constrain_scalars) {
   lub_constrain_tests::expect(x1, lb_inf, ub_inf);
   lub_constrain_tests::expect(x2, lb_inf, ub_inf);
 }
+
+TEST(mathMixMatFun, lub_constrain_fwd_derivative_tails) {
+  using inv_logit_tail_refs::ref;
+  // derivatives of -2 + 5 inv_logit(x)
+  inv_logit_tail_refs::expect_fwd_derivatives(
+      [](const auto& x) { return stan::math::lub_constrain(x, -2.0, 3.0); },
+      [](const ref& r) { return 5 * r.d1; },
+      [](const ref& r) { return 5 * r.d2; });
+}
diff --git a/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp b/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp
new file mode 100644
index 00000000000..9c175a1c417
--- /dev/null
+++ b/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp
@@ -0,0 +1,66 @@
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
+#include <vector>
+
+// Forward-mode derivatives of Phi_approx(x) = inv_logit(0.07056 x^3 + 1.5976 x)
+// come from the fvar overload of inv_logit; they were 0 wherever its own
+// derivative was.
+namespace Phi_approx_fwd_test {
+inline std::vector<double> args() {
+  return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
+}
+// first derivatives by mpmath at 60 digits; at x = 40 and -40 the true
+// value, about 4e-1987, underflows to 0
+inline std::vector<double> derivs() {
+  return {0.0,
+          5.9589784505245725e-37,
+          8.7097694481792379e-21,
+          5.1340418352583332e-15,
+          0.24152728992165639,
+          0.3994,
+          0.24152728992165639,
+          5.1340418352583332e-15,
+          8.7097694481792379e-21,
+          5.9589784505245725e-37,
+          0.0};
+}
+// second derivatives, computed the same way
+inline std::vector<double> second_derivs() {
+  return {0.0,
+          1.3455155964317485e-35,
+          1.2996275105706113e-19,
+          6.018298087621098e-14,
+          0.24180038682517271,
+          0.0,
+          -0.24180038682517271,
+          -6.018298087621098e-14,
+          -1.2996275105706113e-19,
+          -1.3455155964317485e-35,
+          0.0};
+}
+}  // namespace Phi_approx_fwd_test
+
+TEST(mathMixMatFun, PhiApproxFwdDerivativeTails) {
+  using inv_logit_tail_refs::expect_rel;
+  using stan::math::fvar;
+  using stan::math::Phi_approx;
+  using stan::math::var;
+  std::vector<double> xs = Phi_approx_fwd_test::args();
+  std::vector<double> ds = Phi_approx_fwd_test::derivs();
+  std::vector<double> d2s = Phi_approx_fwd_test::second_derivs();
+  for (size_t i = 0; i < xs.size(); ++i) {
+    expect_rel(ds[i], Phi_approx(fvar<double>(xs[i], 1)).d_, xs[i]);
+
+    fvar<fvar<double>> ff(fvar<double>(xs[i], 1), fvar<double>(1, 0));
+    fvar<fvar<double>> y_ff = Phi_approx(ff);
+    expect_rel(ds[i], y_ff.d_.val_, xs[i]);
+    expect_rel(d2s[i], y_ff.d_.d_, xs[i]);
+
+    fvar<var> fv(xs[i], 1);
+    fvar<var> y_fv = Phi_approx(fv);
+    expect_rel(ds[i], y_fv.d_.val(), xs[i]);
+    y_fv.d_.grad();
+    expect_rel(d2s[i], fv.val_.adj(), xs[i]);
+    stan::math::recover_memory();
+  }
+}
diff --git a/test/unit/math/mix/fun/inv_logit_tail_refs.hpp b/test/unit/math/mix/fun/inv_logit_tail_refs.hpp
new file mode 100644
index 00000000000..25163d6b495
--- /dev/null
+++ b/test/unit/math/mix/fun/inv_logit_tail_refs.hpp
@@ -0,0 +1,86 @@
+#ifndef TEST_UNIT_MATH_MIX_FUN_INV_LOGIT_TAIL_REFS_HPP
+#define TEST_UNIT_MATH_MIX_FUN_INV_LOGIT_TAIL_REFS_HPP
+
+#include <stan/math/mix.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+namespace inv_logit_tail_refs {
+
+// inv_logit(x), inv_logit(-x), inv_logit'(x) and inv_logit''(x), computed
+// with mpmath at 60 digits and rounded to double (0 where they underflow).
+struct ref {
+  double x;
+  double s;
+  double s_neg;
+  double d1;
+  double d2;
+};
+
+inline std::vector<ref> refs() {
+  return {{-800, 0, 1.0, 0, 0},
+          {-40, 4.248354255291589e-18, 1.0, 4.248354255291589e-18,
+           4.248354255291589e-18},
+          {-37, 8.533047625744065e-17, 0.9999999999999999,
+           8.533047625744065e-17, 8.533047625744063e-17},
+          {-30, 9.357622968839299e-14, 0.9999999999999064,
+           9.357622968838423e-14, 9.357622968836672e-14},
+          {-20, 2.0611536181902037e-09, 0.9999999979388464,
+           2.061153613941849e-09, 2.061153605445141e-09},
+          {-1, 0.2689414213699951, 0.7310585786300049, 0.19661193324148185,
+           0.09085774767294841},
+          {0, 0.5, 0.5, 0.25, 0},
+          {1, 0.7310585786300049, 0.2689414213699951, 0.19661193324148185,
+           -0.09085774767294841},
+          {20, 0.9999999979388464, 2.0611536181902037e-09,
+           2.061153613941849e-09, -2.061153605445141e-09},
+          {30, 0.9999999999999064, 9.357622968839299e-14, 9.357622968838423e-14,
+           -9.357622968836672e-14},
+          {37, 0.9999999999999999, 8.533047625744065e-17, 8.533047625744065e-17,
+           -8.533047625744063e-17},
+          {40, 1.0, 4.248354255291589e-18, 4.248354255291589e-18,
+           -4.248354255291589e-18},
+          {100, 1.0, 3.720075976020836e-44, 3.720075976020836e-44,
+           -3.720075976020836e-44},
+          {700, 1.0, 9.85967654375977e-305, 9.85967654375977e-305,
+           -9.85967654375977e-305},
+          {800, 1.0, 0, 0, 0}};
+}
+
+// Relative check, absolute below the smallest normal double.
+inline void expect_rel(double expected, double actual, double x,
+                       double tol = 1e-13) {
+  EXPECT_NEAR(
+      expected, actual,
+      tol * std::fmax(std::fabs(expected), std::numeric_limits<double>::min()))
+      << "x = " << x;
+}
+
+// First and second derivatives of f at every reference point through
+// fvar<double>, fvar<fvar<double>> and fvar<var>, against d1(r) and d2(r).
+template <typename F, typename D1, typename D2>
+void expect_fwd_derivatives(const F& f, const D1& d1, const D2& d2) {
+  using stan::math::fvar;
+  using stan::math::var;
+  for (const auto& r : refs()) {
+    expect_rel(d1(r), f(fvar<double>(r.x, 1)).d_, r.x);
+
+    fvar<fvar<double>> ff(fvar<double>(r.x, 1), fvar<double>(1, 0));
+    fvar<fvar<double>> y_ff = f(ff);
+    expect_rel(d1(r), y_ff.d_.val_, r.x);
+    expect_rel(d2(r), y_ff.d_.d_, r.x);
+
+    fvar<var> fv(r.x, 1);
+    fvar<var> y_fv = f(fv);
+    expect_rel(d1(r), y_fv.d_.val(), r.x);
+    y_fv.d_.grad();
+    expect_rel(d2(r), fv.val_.adj(), r.x);
+    stan::math::recover_memory();
+  }
+}
+
+}  // namespace inv_logit_tail_refs
+
+#endif
diff --git a/test/unit/math/mix/fun/inv_logit_test.cpp b/test/unit/math/mix/fun/inv_logit_test.cpp
index 8fee66f7eeb..8a6a058461c 100644
--- a/test/unit/math/mix/fun/inv_logit_test.cpp
+++ b/test/unit/math/mix/fun/inv_logit_test.cpp
@@ -1,4 +1,5 @@
 #include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
 
 TEST(mathMixMatFun, invLogit) {
   auto f = [](const auto& x1) { return stan::math::inv_logit(x1); };
@@ -12,3 +13,87 @@ TEST(mathMixMatFun, invLogit) {
   stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args));
   stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args));
 }
+
+TEST(mathMixMatFun, invLogitDerivativesLargeArgs) {
+  using inv_logit_tail_refs::ref;
+  for (const auto& r : inv_logit_tail_refs::refs()) {
+    stan::math::var a = r.x;
+    stan::math::inv_logit(a).grad();
+    inv_logit_tail_refs::expect_rel(r.d1, a.adj(), r.x);
+    stan::math::recover_memory();
+  }
+  inv_logit_tail_refs::expect_fwd_derivatives(
+      [](const auto& x) { return stan::math::inv_logit(x); },
+      [](const ref& r) { return r.d1; }, [](const ref& r) { return r.d2; });
+}
+
+TEST(mathMixMatFun, invLogitDerivativesLargeArgsContainers) {
+  using inv_logit_tail_refs::expect_rel;
+  using stan::math::fvar;
+  using stan::math::inv_logit;
+  using stan::math::sum;
+  using stan::math::var;
+  std::vector<inv_logit_tail_refs::ref> ps = inv_logit_tail_refs::refs();
+  int n = ps.size();
+  Eigen::VectorXd x(n);
+  Eigen::VectorXd d1(n);
+  Eigen::VectorXd d2(n);
+  for (int i = 0; i < n; ++i) {
+    x(i) = ps[i].x;
+    d1(i) = ps[i].d1;
+    d2(i) = ps[i].d2;
+  }
+
+  stan::math::var_value<Eigen::VectorXd> xv(x);
+  sum(inv_logit(xv)).grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), xv.adj()(i), x(i));
+  }
+  stan::math::recover_memory();
+
+  stan::math::var_value<Eigen::RowVectorXd> xr(x.transpose());
+  sum(inv_logit(xr)).grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), xr.adj()(i), x(i));
+  }
+  stan::math::recover_memory();
+
+  // inv_logit' is even
+  Eigen::MatrixXd x2(n, 2);
+  x2 << x, -x;
+  stan::math::var_value<Eigen::MatrixXd> xmv(x2);
+  sum(inv_logit(xmv)).grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), xmv.adj()(i, 0), x(i));
+    expect_rel(d1(i), xmv.adj()(i, 1), -x(i));
+  }
+  stan::math::recover_memory();
+
+  Eigen::Matrix<var, Eigen::Dynamic, 1> xm = x;
+  sum(inv_logit(xm)).grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), xm(i).adj(), x(i));
+  }
+  stan::math::recover_memory();
+
+  std::vector<var> xs(x.data(), x.data() + n);
+  sum(inv_logit(xs)).grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), xs[i].adj(), x(i));
+  }
+  stan::math::recover_memory();
+
+  Eigen::Matrix<fvar<var>, Eigen::Dynamic, 1> xf(n);
+  for (int i = 0; i < n; ++i) {
+    xf(i) = fvar<var>(x(i), 1);
+  }
+  Eigen::Matrix<fvar<var>, Eigen::Dynamic, 1> yf = inv_logit(xf);
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d1(i), yf(i).d_.val(), x(i));
+  }
+  sum(yf).d_.grad();
+  for (int i = 0; i < n; ++i) {
+    expect_rel(d2(i), xf(i).val_.adj(), x(i));
+  }
+  stan::math::recover_memory();
+}
diff --git a/test/unit/math/mix/fun/log1m_inv_logit_test.cpp b/test/unit/math/mix/fun/log1m_inv_logit_test.cpp
index 175ead491ac..6e730bfb0ae 100644
--- a/test/unit/math/mix/fun/log1m_inv_logit_test.cpp
+++ b/test/unit/math/mix/fun/log1m_inv_logit_test.cpp
@@ -1,4 +1,5 @@
 #include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
 
 TEST(mathMixMatFun, log1mInvLogit) {
   auto f = [](const auto& x1) { return stan::math::log1m_inv_logit(x1); };
@@ -25,3 +26,11 @@ TEST(mathMixMatFun, log1minvlogit_varmat) {
   }
   expect_ad_vector_matvar(f, A);
 }
+
+TEST(mathMixMatFun, log1mInvLogitFwdDerivativeTails) {
+  using inv_logit_tail_refs::ref;
+  // d/dx = -inv_logit(x), d^2/dx^2 = -inv_logit'(x)
+  inv_logit_tail_refs::expect_fwd_derivatives(
+      [](const auto& x) { return stan::math::log1m_inv_logit(x); },
+      [](const ref& r) { return -r.s; }, [](const ref& r) { return -r.d1; });
+}
diff --git a/test/unit/math/mix/fun/log1p_exp_test.cpp b/test/unit/math/mix/fun/log1p_exp_test.cpp
index 9641c036c0e..13167f8a464 100644
--- a/test/unit/math/mix/fun/log1p_exp_test.cpp
+++ b/test/unit/math/mix/fun/log1p_exp_test.cpp
@@ -1,4 +1,7 @@
 #include <test/unit/math/test_ad.hpp>
+#include <cmath>
+#include <limits>
+#include <vector>
 
 TEST(mathMixMatFun, log1pExp) {
   auto f = [](const auto& x1) { return stan::math::log1p_exp(x1); };
@@ -13,3 +16,46 @@ TEST(mathMixMatFun, log1pExp) {
   stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args));
   stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args));
 }
+
+TEST(mathMixMatFun, log1pExpHigherDerivativesLargeArgs) {
+  using stan::math::fvar;
+  using stan::math::log1p_exp;
+  using stan::math::var;
+  // {x, log1p_exp''(x), log1p_exp'''(x)}: mpmath, 50 digits; 0 where the
+  // value underflows (e^-800)
+  std::vector<std::vector<double>> points{
+      {-800, 0, 0},
+      {-40, 4.2483542552915889592e-18, 4.2483542552915889231e-18},
+      {-20, 2.0611536139418493437e-9, 2.0611536054451408856e-9},
+      {-1, 1.9661193324148185254e-1, 9.0857747672948409442e-2},
+      {0, 0.25, 0},
+      {1, 1.9661193324148185254e-1, -9.0857747672948409442e-2},
+      {20, 2.0611536139418493437e-9, -2.0611536054451408856e-9},
+      {30, 9.3576229688384233028e-14, -9.3576229688366720006e-14},
+      {37, 8.533047625744064338e-17, -8.5330476257440628818e-17},
+      {40, 4.2483542552915889592e-18, -4.2483542552915889231e-18},
+      {100, 3.720075976020835963e-44, -3.720075976020835963e-44},
+      {700, 9.8596765437597708567e-305, -9.8596765437597708567e-305},
+      {800, 0, 0}};
+  for (const auto& p : points) {
+    double x = p[0];
+    double d2 = p[1];
+    double d3 = p[2];
+    double tol2 = 1e-13 * std::fmax(d2, std::numeric_limits<double>::min());
+    double tol3
+        = 1e-13 * std::fmax(std::fabs(d3), std::numeric_limits<double>::min());
+
+    fvar<fvar<double>> ff(fvar<double>(x, 1), fvar<double>(1, 0));
+    EXPECT_NEAR(log1p_exp(ff).d_.d_, d2, tol2) << "x = " << x;
+
+    fvar<var> fv(x, 1);
+    log1p_exp(fv).d_.grad();
+    EXPECT_NEAR(fv.val_.adj(), d2, tol2) << "x = " << x;
+    stan::math::recover_memory();
+
+    fvar<fvar<var>> ffv(fvar<var>(x, 1), fvar<var>(1, 0));
+    log1p_exp(ffv).d_.d_.grad();
+    EXPECT_NEAR(ffv.val_.val_.adj(), d3, tol3) << "x = " << x;
+    stan::math::recover_memory();
+  }
+}
diff --git a/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp b/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp
index c5c4f31f474..ca7bde6a05b 100644
--- a/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp
+++ b/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp
@@ -1,5 +1,6 @@
 #include <test/unit/math/test_ad.hpp>
 #include <limits>
+#include <vector>
 
 TEST(mathMixScalFun, logInvLogitDiff) {
   auto f = [](const auto& x1, const auto& x2) {
@@ -28,3 +29,67 @@ TEST(mathMixScalFun, logInvLogitDiff_vec) {
   in2 << 0.5, -3.4;
   stan::test::expect_ad_vectorized_binary(f, in1, in2);
 }
+
+TEST(mathMixScalFun, logInvLogitDiffGradientLargeArgs) {
+  using stan::math::fvar;
+  using stan::math::log_inv_logit_diff;
+  using stan::math::var;
+  // {x, y, d/dx, d/dy}: mpmath, 50 digits
+  std::vector<std::vector<double>> points{
+      {40, 0, 8.4967085105831779907e-18, -5.0000000000000000425e-1},
+      {100, 1, 1.3832290902125321262e-43, -7.3105857863000487925e-1},
+      {30, -30, 9.3576229688401746049e-14, -9.3576229688401746049e-14},
+      {37, 36, 5.8197670686932650972e-1, -1.5819767068693261924},
+      {-30, -40, 1.0000454019909161115, -4.5401991009692016683e-5},
+      {1, 0, 8.5091812823932154513e-1, -1.0819767068693264244},
+      {700, 20, 4.7835719068902113962e-296, -9.9999999793884638181e-1},
+      {0.5, -1.0, 6.647575855870136797e-1, -5.5615833815886336509e-1}};
+  for (const auto& p : points) {
+    double x = p[0];
+    double y = p[1];
+    double dx = p[2];
+    double dy = p[3];
+
+    var xv = x;
+    var yv = y;
+    log_inv_logit_diff(xv, yv).grad();
+    EXPECT_NEAR(xv.adj(), dx, 1e-13 * dx) << x << ", " << y;
+    EXPECT_NEAR(yv.adj(), dy, -1e-13 * dy) << x << ", " << y;
+    stan::math::recover_memory();
+
+    var xv2 = x;
+    log_inv_logit_diff(xv2, y).grad();
+    EXPECT_NEAR(xv2.adj(), dx, 1e-13 * dx) << x << ", " << y;
+    stan::math::recover_memory();
+
+    EXPECT_NEAR(log_inv_logit_diff(fvar<double>(x, 1), y).d_, dx, 1e-13 * dx)
+        << x << ", " << y;
+    EXPECT_NEAR(log_inv_logit_diff(fvar<double>(x, 1), fvar<double>(y, 0)).d_,
+                dx, 1e-13 * dx)
+        << x << ", " << y;
+    EXPECT_NEAR(log_inv_logit_diff(fvar<double>(x, 0), fvar<double>(y, 1)).d_,
+                dy, -1e-13 * dy)
+        << x << ", " << y;
+  }
+}
+
+TEST(mathMixScalFun, logInvLogitDiffGradientEqualArgs) {
+  using stan::math::fvar;
+  using stan::math::log_inv_logit_diff;
+  using stan::math::var;
+  // behaviour change: d/dx at x == y is +inf, the limit from x > y (was -inf)
+  double inf = std::numeric_limits<double>::infinity();
+  var xv = 2.0;
+  var yv = 2.0;
+  log_inv_logit_diff(xv, yv).grad();
+  EXPECT_EQ(xv.adj(), inf);
+  EXPECT_EQ(yv.adj(), -inf);
+  stan::math::recover_memory();
+
+  var xv2 = 2.0;
+  log_inv_logit_diff(xv2, 2.0).grad();
+  EXPECT_EQ(xv2.adj(), inf);
+  stan::math::recover_memory();
+
+  EXPECT_EQ(log_inv_logit_diff(fvar<double>(2, 1), 2.0).d_, inf);
+}
diff --git a/test/unit/math/mix/fun/log_inv_logit_test.cpp b/test/unit/math/mix/fun/log_inv_logit_test.cpp
index f94478e21c9..3bad8322319 100644
--- a/test/unit/math/mix/fun/log_inv_logit_test.cpp
+++ b/test/unit/math/mix/fun/log_inv_logit_test.cpp
@@ -1,4 +1,5 @@
 #include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
 
 TEST(mathMixMatFun, logInvLogit) {
   auto f = [](const auto& x1) { return stan::math::log_inv_logit(x1); };
@@ -25,3 +26,11 @@ TEST(mathMixMatFun, logInvLogitVarMat) {
   }
   expect_ad_vector_matvar(f, A);
 }
+
+TEST(mathMixMatFun, logInvLogitFwdDerivativeTails) {
+  using inv_logit_tail_refs::ref;
+  // d/dx = inv_logit(-x), d^2/dx^2 = -inv_logit'(x)
+  inv_logit_tail_refs::expect_fwd_derivatives(
+      [](const auto& x) { return stan::math::log_inv_logit(x); },
+      [](const ref& r) { return r.s_neg; }, [](const ref& r) { return -r.d1; });
+}
diff --git a/test/unit/math/opencl/rev/inv_logit_test.cpp b/test/unit/math/opencl/rev/inv_logit_test.cpp
index b8af60f0778..8549b411aef 100644
--- a/test/unit/math/opencl/rev/inv_logit_test.cpp
+++ b/test/unit/math/opencl/rev/inv_logit_test.cpp
@@ -27,4 +27,38 @@ TEST(OpenCLinv_logit, prim_rev_values_large) {
   stan::math::test::compare_cpu_opencl_prim_rev(inv_logit_functor, a);
 }
 
+TEST(OpenCLinv_logit, rev_derivative_large_args) {
+  using stan::math::matrix_cl;
+  using stan::math::var_value;
+  Eigen::VectorXd x(8);
+  x << -40, -20, 0, 20, 30, 37, 40, 100;
+  // inv_logit'(x): mpmath, 50 digits
+  Eigen::VectorXd d(8);
+  d << 4.2483542552915889592e-18, 2.0611536139418493437e-9, 0.25,
+      2.0611536139418493437e-9, 9.3576229688384233028e-14,
+      8.533047625744064338e-17, 4.2483542552915889592e-18,
+      3.720075976020835963e-44;
+  const var_value<matrix_cl<double>> x_cl = stan::math::to_matrix_cl(x);
+  stan::math::sum(stan::math::inv_logit(x_cl)).grad();
+  Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj());
+  for (int i = 0; i < x.size(); ++i) {
+    EXPECT_NEAR(adj(i), d(i), 1e-13 * d(i)) << "x = " << x(i);
+  }
+  stan::math::recover_memory();
+}
+
+TEST(OpenCLinv_logit, rev_non_const_lvalue) {
+  using stan::math::matrix_cl;
+  using stan::math::var_value;
+  Eigen::VectorXd x(2);
+  x << 0, 40;
+  var_value<matrix_cl<double>> x_cl = stan::math::to_matrix_cl(x);
+  var_value<matrix_cl<double>> y_cl = stan::math::inv_logit(x_cl);
+  stan::math::sum(y_cl).grad();
+  Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj());
+  EXPECT_NEAR(adj(0), 0.25, 1e-15);
+  EXPECT_NEAR(adj(1), 4.2483542552915889592e-18, 1e-30);
+  stan::math::recover_memory();
+}
+
 #endif
diff --git a/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp b/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp
index b8df60796da..3827d77ad0f 100644
--- a/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp
+++ b/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp
@@ -57,4 +57,34 @@ TEST(OpenCLMatrix_log_inv_logit_diff, prim_rev_scalar_values_large) {
                                                 a);
 }
 
+TEST(OpenCLMatrix_log_inv_logit_diff, rev_gradient_large_args) {
+  using stan::math::matrix_cl;
+  using stan::math::var_value;
+  Eigen::VectorXd x(8);
+  Eigen::VectorXd y(8);
+  x << 40, 100, 30, 37, -30, 1, 700, 0.5;
+  y << 0, 1, -30, 36, -40, 0, 20, -1;
+  // d/dx and d/dy: mpmath, 50 digits
+  Eigen::VectorXd dx(8);
+  Eigen::VectorXd dy(8);
+  dx << 8.4967085105831779907e-18, 1.3832290902125321262e-43,
+      9.3576229688401746049e-14, 5.8197670686932650972e-1,
+      1.0000454019909161115, 8.5091812823932154513e-1,
+      4.7835719068902113962e-296, 6.647575855870136797e-1;
+  dy << -5.0000000000000000425e-1, -7.3105857863000487925e-1,
+      -9.3576229688401746049e-14, -1.5819767068693261924,
+      -4.5401991009692016683e-5, -1.0819767068693264244,
+      -9.9999999793884638181e-1, -5.5615833815886336509e-1;
+  var_value<matrix_cl<double>> x_cl = stan::math::to_matrix_cl(x);
+  var_value<matrix_cl<double>> y_cl = stan::math::to_matrix_cl(y);
+  stan::math::sum(stan::math::log_inv_logit_diff(x_cl, y_cl)).grad();
+  Eigen::VectorXd x_adj = stan::math::from_matrix_cl(x_cl.adj());
+  Eigen::VectorXd y_adj = stan::math::from_matrix_cl(y_cl.adj());
+  for (int i = 0; i < x.size(); ++i) {
+    EXPECT_NEAR(x_adj(i), dx(i), 1e-13 * dx(i)) << x(i) << ", " << y(i);
+    EXPECT_NEAR(y_adj(i), dy(i), -1e-13 * dy(i)) << x(i) << ", " << y(i);
+  }
+  stan::math::recover_memory();
+}
+
 #endif

From 2eeebdc17972a4a296618eb10ec9c023d14f28d2 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 11:25:53 +0200
Subject: [PATCH 04/10] Accurate inv_logit and log_inv_logit_diff derivatives
 at large arguments

---
 stan/math/fwd/fun/inv_logit.hpp             |  5 ++-
 stan/math/fwd/fun/log_inv_logit_diff.hpp    | 11 ++---
 stan/math/opencl/rev/inv_logit.hpp          |  2 +-
 stan/math/opencl/rev/log_inv_logit_diff.hpp | 10 ++---
 stan/math/prim/fun/log_inv_logit_diff.hpp   |  4 +-
 stan/math/rev/fun/inv_logit.hpp             | 46 +++++++++++++--------
 stan/math/rev/fun/log_inv_logit_diff.hpp    | 16 +++----
 7 files changed, 53 insertions(+), 41 deletions(-)

diff --git a/stan/math/fwd/fun/inv_logit.hpp b/stan/math/fwd/fun/inv_logit.hpp
index 325a74a133a..2af1eb48962 100644
--- a/stan/math/fwd/fun/inv_logit.hpp
+++ b/stan/math/fwd/fun/inv_logit.hpp
@@ -18,8 +18,9 @@ namespace math {
  */
 template <typename T, require_fvar_t<T>* = nullptr>
 inline auto inv_logit(T&& x) {
-  return std::decay_t<T>(inv_logit(x.val_),
-                         x.d_ * inv_logit(x.val_) * (1 - inv_logit(x.val_)));
+  const auto val = inv_logit(x.val_);
+  // inv_logit(-x) rather than 1 - inv_logit(x), which cancels for x >> 0
+  return std::decay_t<T>(val, x.d_ * val * inv_logit(-x.val_));
 }
 
 }  // namespace math
diff --git a/stan/math/fwd/fun/log_inv_logit_diff.hpp b/stan/math/fwd/fun/log_inv_logit_diff.hpp
index e31f166d040..e7e962c5a00 100644
--- a/stan/math/fwd/fun/log_inv_logit_diff.hpp
+++ b/stan/math/fwd/fun/log_inv_logit_diff.hpp
@@ -21,11 +21,11 @@ namespace math {
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1}
+    \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1}
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1}
+    \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1}
    \f]
  *
  * @tparam T inner type of the fvar
@@ -36,16 +36,17 @@ namespace math {
  */
 template <typename T>
 inline fvar<T> log_inv_logit_diff(const fvar<T>& x, const fvar<T>& y) {
+  // both partials as sums of positive terms, no cancellation for x >> 0
+  const auto c = inv(expm1(x.val_ - y.val_));
   return fvar<T>(
       log_inv_logit_diff(x.val_, y.val_),
-      -x.d_ * (inv(expm1(y.val_ - x.val_)) + inv_logit(x.val_))
-          - y.d_ * (inv(expm1(x.val_ - y.val_)) + inv_logit(y.val_)));
+      x.d_ * (c + inv_logit(-x.val_)) - y.d_ * (c + inv_logit(y.val_)));
 }
 
 template <typename T>
 inline fvar<T> log_inv_logit_diff(const fvar<T>& x, double y) {
   return fvar<T>(log_inv_logit_diff(x.val_, y),
-                 -x.d_ * (inv(expm1(y - x.val_)) + inv_logit(x.val_)));
+                 x.d_ * (inv(expm1(x.val_ - y)) + inv_logit(-x.val_)));
 }
 
 template <typename T>
diff --git a/stan/math/opencl/rev/inv_logit.hpp b/stan/math/opencl/rev/inv_logit.hpp
index aef65ee3c4a..7602590bef6 100644
--- a/stan/math/opencl/rev/inv_logit.hpp
+++ b/stan/math/opencl/rev/inv_logit.hpp
@@ -21,7 +21,7 @@ inline var_value<matrix_cl<double>> inv_logit(const var_value<T>& A) {
   return make_callback_var(
       inv_logit(A.val()), [A](vari_value<matrix_cl<double>>& res) mutable {
         A.adj() += elt_multiply(res.adj(),
-                                elt_multiply(res.val(), 1.0 - res.val()));
+                                elt_multiply(res.val(), inv_logit(-A.val())));
       });
 }
 
diff --git a/stan/math/opencl/rev/log_inv_logit_diff.hpp b/stan/math/opencl/rev/log_inv_logit_diff.hpp
index e8eb5381f7a..82e87a534b0 100644
--- a/stan/math/opencl/rev/log_inv_logit_diff.hpp
+++ b/stan/math/opencl/rev/log_inv_logit_diff.hpp
@@ -38,13 +38,11 @@ inline var_value<matrix_cl<double>> log_inv_logit_diff(T_x&& x, T_y&& y) {
   return make_callback_var(
       log_inv_logit_diff(value_of(x_arena), value_of(y_arena)),
       [x_arena, y_arena](const vari_value<matrix_cl<double>>& res) mutable {
+        // both partials as sums of positive terms, no cancellation for x >> 0
+        auto c = inv(expm1(value_of(x_arena) - value_of(y_arena)));
         adjoint_results(x_arena, y_arena) += expressions(
-            -elt_multiply(res.adj(),
-                          inv(expm1(value_of(y_arena) - value_of(x_arena)))
-                              + inv_logit(value_of(x_arena))),
-            -elt_multiply(res.adj(),
-                          inv(expm1(value_of(x_arena) - value_of(y_arena)))
-                              + inv_logit(value_of(y_arena))));
+            elt_multiply(res.adj(), c + inv_logit(-value_of(x_arena))),
+            -elt_multiply(res.adj(), c + inv_logit(value_of(y_arena))));
       });
 }
 
diff --git a/stan/math/prim/fun/log_inv_logit_diff.hpp b/stan/math/prim/fun/log_inv_logit_diff.hpp
index 96bce37601b..ff12a5bb34f 100644
--- a/stan/math/prim/fun/log_inv_logit_diff.hpp
+++ b/stan/math/prim/fun/log_inv_logit_diff.hpp
@@ -20,11 +20,11 @@ namespace math {
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1}
+    \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1}
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1}
+    \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1}
    \f]
  *
  * @tparam T1 type of x argument
diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp
index 13b0e051b65..e739679132e 100644
--- a/stan/math/rev/fun/inv_logit.hpp
+++ b/stan/math/rev/fun/inv_logit.hpp
@@ -3,7 +3,11 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/prim/fun/abs.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
+#include <stan/math/prim/functor/apply_scalar_binary.hpp>
 
 namespace stan {
 namespace math {
@@ -15,20 +19,27 @@ namespace math {
  *
  * The derivative of inverse logit is
  *
- * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \mbox{logit}^{-1}(x) (1 -
- * \mbox{logit}^{-1}(x))\f$.
+ * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \frac{e}{(1 + e)^2}\f$
+ * with \f$e = \exp(-|x|)\f$.
  *
- * @tparam T Arithmetic or a type inheriting from `EigenBase`.
+ * @tparam T a `var` or a `var_value` of a matrix type
  * @param a Argument variable.
- * @return Inverse logit of argument.
+ * @return Inverse logit of argument, elementwise for a matrix.
  */
-template <typename T, require_var_t<T>* = nullptr>
+template <
+    typename T, require_var_t<T>* = nullptr,
+    require_all_not_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
 inline auto inv_logit(T&& a) {
-  return make_callback_var(inv_logit(a.val()), [a](auto& vi) mutable {
-    as_array_or_scalar(a).adj() += as_array_or_scalar(vi.adj())
-                                   * as_array_or_scalar(vi.val())
-                                   * (1.0 - as_array_or_scalar(vi.val()));
-  });
+  auto e = to_arena(exp(-abs(a.val())));
+  return make_callback_var(
+      apply_scalar_binary(
+          [](double x, double e) { return (x >= 0 ? 1.0 : e) / (1.0 + e); },
+          a.val(), e),
+      [a, e](auto& vi) mutable {
+        const auto& e_a = as_array_or_scalar(e);
+        as_array_or_scalar(a.adj())
+            += as_array_or_scalar(vi.adj()) * e_a / ((1.0 + e_a) * (1.0 + e_a));
+      });
 }
 
 /**
@@ -38,20 +49,21 @@ inline auto inv_logit(T&& a) {
  *
  * The derivative of inverse logit is
  *
- * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \mbox{logit}^{-1}(x) (1 -
- * \mbox{logit}^{-1}(x))\f$.
+ * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \frac{e}{(1 + e)^2}\f$
+ * with \f$e = \exp(-|x|)\f$.
  *
  * @tparam T type of Eigen expression
  * @param x Eigen expression
- * @return Inverse logit of argument.
+ * @return Inverse logit of each value in x.
  */
 template <typename T, require_eigen_vt<is_var, T>* = nullptr>
 inline auto inv_logit(T&& x) {
   auto x_arena = to_arena(std::forward<T>(x));
-  arena_t<T> ret = inv_logit(x_arena.val());
-  reverse_pass_callback([x_arena, ret]() mutable {
-    x_arena.adj().array()
-        += ret.adj().array() * ret.val().array() * (1.0 - ret.val().array());
+  const auto x_val = x_arena.val().array().eval();
+  arena_t<plain_type_t<decltype(x_val)>> e = (-x_val.abs()).exp();
+  arena_t<T> ret = ((x_val >= 0).select(1.0, e) / (1.0 + e)).matrix();
+  reverse_pass_callback([x_arena, ret, e]() mutable {
+    x_arena.adj().array() += ret.adj().array() * e / (1.0 + e).square();
   });
   return ret;
 }
diff --git a/stan/math/rev/fun/log_inv_logit_diff.hpp b/stan/math/rev/fun/log_inv_logit_diff.hpp
index 132c5d2c85b..b502c3a0614 100644
--- a/stan/math/rev/fun/log_inv_logit_diff.hpp
+++ b/stan/math/rev/fun/log_inv_logit_diff.hpp
@@ -21,11 +21,11 @@ namespace math {
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1}
+    \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1}
    \f]
 
    \f[
-    \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1}
+    \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1}
    \f]
  *
  * @tparam T1 Type of x argument
@@ -41,11 +41,10 @@ class log_inv_logit_diff_vv_vari : public op_vv_vari {
   log_inv_logit_diff_vv_vari(vari* avi, vari* bvi)
       : op_vv_vari(log_inv_logit_diff(avi->val_, bvi->val_), avi, bvi) {}
   void chain() {
-    avi_->adj_
-        -= adj_ * (inv(expm1(bvi_->val_ - avi_->val_)) + inv_logit(avi_->val_));
-
-    bvi_->adj_
-        -= adj_ * (inv(expm1(avi_->val_ - bvi_->val_)) + inv_logit(bvi_->val_));
+    // both partials as sums of positive terms, no cancellation for x >> 0
+    double c = inv(expm1(avi_->val_ - bvi_->val_));
+    avi_->adj_ += adj_ * (c + inv_logit(-avi_->val_));
+    bvi_->adj_ -= adj_ * (c + inv_logit(bvi_->val_));
   }
 };
 
@@ -54,7 +53,8 @@ class log_inv_logit_diff_vd_vari : public op_vd_vari {
   log_inv_logit_diff_vd_vari(vari* avi, double b)
       : op_vd_vari(log_inv_logit_diff(avi->val_, b), avi, b) {}
   void chain() {
-    avi_->adj_ -= adj_ * (inv(expm1(bd_ - avi_->val_)) + inv_logit(avi_->val_));
+    avi_->adj_
+        += adj_ * (inv(expm1(avi_->val_ - bd_)) + inv_logit(-avi_->val_));
   }
 };
 

From 36694eb9338ea63acb3ebc8c2b840866d146c60b Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 21:09:54 +0200
Subject: [PATCH 05/10] Phi_approx: name the coefficients in constants.hpp and
 use the double-precision minimax optimum

---
 stan/math/fwd/fun/Phi_approx.hpp              |  4 ++-
 .../kernels/device_functions/Phi_approx.hpp   |  7 +++--
 stan/math/opencl/rev/Phi_approx.hpp           |  6 +++--
 stan/math/prim/fun/Phi_approx.hpp             |  7 ++++-
 stan/math/prim/fun/constants.hpp              | 14 ++++++++++
 stan/math/rev/fun/Phi_approx.hpp              | 27 +++++++++++++------
 test/unit/math/mix/fun/Phi_approx_test.cpp    | 21 ++++++++-------
 test/unit/math/opencl/rev/Phi_approx_test.cpp |  6 ++---
 8 files changed, 65 insertions(+), 27 deletions(-)

diff --git a/stan/math/fwd/fun/Phi_approx.hpp b/stan/math/fwd/fun/Phi_approx.hpp
index 2509a44348a..00bb2fb7b2c 100644
--- a/stan/math/fwd/fun/Phi_approx.hpp
+++ b/stan/math/fwd/fun/Phi_approx.hpp
@@ -6,6 +6,7 @@
 #include <stan/math/fwd/fun/inv_logit.hpp>
 #include <stan/math/fwd/fun/pow.hpp>
 #include <stan/math/prim/fun/Phi_approx.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <cmath>
 
 namespace stan {
@@ -23,7 +24,8 @@ namespace math {
  */
 template <typename T>
 inline fvar<T> Phi_approx(const fvar<T>& x) {
-  return inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
+  return inv_logit(PHI_APPROX_CUBIC_COEF * pow(x, 3.0)
+                   + PHI_APPROX_LINEAR_COEF * x);
 }
 
 }  // namespace math
diff --git a/stan/math/opencl/kernels/device_functions/Phi_approx.hpp b/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
index a86c3563951..cb635ebf482 100644
--- a/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
+++ b/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
@@ -17,12 +17,15 @@ static constexpr const char* phi_approx_device_function
           // \endcond
           /** \ingroup opencl_kernels
            * Return the Phi_approx function applied to the specified
-           * argument.
+           * argument. The coefficients are PHI_APPROX_CUBIC_COEF and
+           * PHI_APPROX_LINEAR_COEF of stan/math/prim/fun/constants.hpp,
+           * repeated here because a kernel source string cannot use them.
            *
            * @return Phi_approx(x)
            */
           inline double Phi_approx(double x) {
-            return inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
+            return inv_logit(0.07056409656566716 * pow(x, 3.0)
+                             + 1.597602874240784 * x);
           }
           // \cond
           ) "\n#endif\n";  // NOLINT
diff --git a/stan/math/opencl/rev/Phi_approx.hpp b/stan/math/opencl/rev/Phi_approx.hpp
index 708d435b9ed..69870016d0d 100644
--- a/stan/math/opencl/rev/Phi_approx.hpp
+++ b/stan/math/opencl/rev/Phi_approx.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/opencl/kernel_generator.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 
 namespace stan {
 namespace math {
@@ -22,10 +23,11 @@ inline var_value<matrix_cl<double>> Phi_approx(const var_value<T>& A) {
       Phi_approx(A.val()), [A](vari_value<matrix_cl<double>>& res) mutable {
         // e / (1 + e)^2, not res * (1 - res), which cancels for u >> 0
         auto x_sq = square(A.val());
-        auto e = exp(-fabs(elt_multiply(A.val(), 0.07056 * x_sq + 1.5976)));
+        auto e = exp(-fabs(elt_multiply(
+            A.val(), PHI_APPROX_CUBIC_COEF * x_sq + PHI_APPROX_LINEAR_COEF)));
         A.adj() += elt_multiply(
             elt_multiply(res.adj(), elt_divide(e, square(1.0 + e))),
-            3.0 * 0.07056 * x_sq + 1.5976);
+            3.0 * PHI_APPROX_CUBIC_COEF * x_sq + PHI_APPROX_LINEAR_COEF);
       });
 }
 
diff --git a/stan/math/prim/fun/Phi_approx.hpp b/stan/math/prim/fun/Phi_approx.hpp
index 8351cb13e49..b6fa81716d1 100644
--- a/stan/math/prim/fun/Phi_approx.hpp
+++ b/stan/math/prim/fun/Phi_approx.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_FUN_PHI_APPROX_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <cmath>
@@ -12,6 +13,9 @@ namespace math {
 /**
  * Return an approximation of the unit normal CDF.
  *
+ * Computes \f$ \mbox{logit}^{-1}(\gamma x^3 + \beta x) \f$ with the
+ * coefficients <code>PHI_APPROX_CUBIC_COEF</code> and
+ * <code>PHI_APPROX_LINEAR_COEF</code>, the form of Bowling et al. (2009),
  * http://www.jiem.org/index.php/jiem/article/download/60/27
  *
  * This function can be used to implement the inverse link function
@@ -22,7 +26,8 @@ namespace math {
  */
 inline double Phi_approx(double x) {
   using std::pow;
-  return inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
+  return inv_logit(PHI_APPROX_CUBIC_COEF * pow(x, 3.0)
+                   + PHI_APPROX_LINEAR_COEF * x);
 }
 
 /**
diff --git a/stan/math/prim/fun/constants.hpp b/stan/math/prim/fun/constants.hpp
index fa2b64a3cde..38e0684430f 100644
--- a/stan/math/prim/fun/constants.hpp
+++ b/stan/math/prim/fun/constants.hpp
@@ -191,6 +191,20 @@ const double NEG_LOG_SQRT_TWO_PI = -std::log(SQRT_TWO_PI);
  */
 const double POISSON_MAX_RATE = std::pow(2.0, 30);
 
+/**
+ * Coefficients of the logistic approximation to the standard normal CDF
+ * computed by <code>Phi_approx</code>,
+ * \f$ \Phi(x) \approx \mbox{logit}^{-1}(\gamma x^3 + \beta x) \f$,
+ * the form of Bowling et al. (2009). They are the pair that minimises the
+ * maximum absolute deviation from \f$ \Phi \f$ over the real line, rounded
+ * to double precision: the deviation is 1.40313e-4, attained at
+ * \f$ x = \pm 0.4816, \pm 1.4764, \pm 2.6876 \f$. The paper prints the pair
+ * rounded to 0.07056 and 1.5976 (maximum deviation 1.41448e-4). The OpenCL
+ * device function repeats the values as literals.
+ */
+static constexpr double PHI_APPROX_CUBIC_COEF = 0.07056409656566716;
+static constexpr double PHI_APPROX_LINEAR_COEF = 1.597602874240784;
+
 /**
  * Return positive infinity.
  *
diff --git a/stan/math/rev/fun/Phi_approx.hpp b/stan/math/rev/fun/Phi_approx.hpp
index 1d3fe29af20..2d40ae6b7a1 100644
--- a/stan/math/rev/fun/Phi_approx.hpp
+++ b/stan/math/rev/fun/Phi_approx.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/prim/fun/Phi_approx.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <cmath>
 
 namespace stan {
@@ -31,27 +32,33 @@ namespace math {
    \f]
 
    \f[
-   \Phi_{\mbox{\footnotesize approx}}(x) = \mbox{logit}^{-1}(0.07056 \,
-   x^3 + 1.5976 \, x)
+   \Phi_{\mbox{\footnotesize approx}}(x) = \mbox{logit}^{-1}(\gamma \,
+   x^3 + \beta \, x)
    \f]
 
    \f[
    \frac{\partial \, \Phi_{\mbox{\footnotesize approx}}(x)}{\partial x}
-   = -\Phi_{\mbox{\footnotesize approx}}^2(x)
-   e^{-0.07056x^3 - 1.5976x}(-0.21168x^2-1.5976)
+   = \Phi_{\mbox{\footnotesize approx}}(x)
+   \left(1 - \Phi_{\mbox{\footnotesize approx}}(x)\right)
+   (3 \gamma x^2 + \beta)
    \f]
+
+   with \f$ \gamma \f$ = <code>PHI_APPROX_CUBIC_COEF</code> and
+   \f$ \beta \f$ = <code>PHI_APPROX_LINEAR_COEF</code>.
  *
  * @param a Variable argument.
  * @return The corresponding unit normal cdf approximation.
  */
 inline var Phi_approx(const var& a) {
   double av_squared = a.val() * a.val();
-  double u = 0.07056 * a.val() * av_squared + 1.5976 * a.val();
+  double u = PHI_APPROX_CUBIC_COEF * a.val() * av_squared
+             + PHI_APPROX_LINEAR_COEF * a.val();
   // e / (1 + e)^2, not f * (1 - f), which cancels for u >> 0
   double e = std::exp(-std::fabs(u));
   double f = (u >= 0 ? 1.0 : e) / (1.0 + e);
   double da
-      = e / ((1.0 + e) * (1.0 + e)) * (3.0 * 0.07056 * av_squared + 1.5976);
+      = e / ((1.0 + e) * (1.0 + e))
+        * (3.0 * PHI_APPROX_CUBIC_COEF * av_squared + PHI_APPROX_LINEAR_COEF);
   return make_callback_var(
       f, [a, da](auto& vi) mutable { a.adj() += vi.adj() * da; });
 }
@@ -62,10 +69,14 @@ inline auto Phi_approx(const T& a) {
   arena_t<value_type_t<T>> f(a.rows(), a.cols());
   arena_t<value_type_t<T>> da(a.rows(), a.cols());
   // e / (1 + e)^2, not f * (1 - f), which cancels for u >> 0
-  da.array() = (-(0.07056 * x * x.square() + 1.5976 * x).abs()).exp();
+  da.array()
+      = (-(PHI_APPROX_CUBIC_COEF * x * x.square() + PHI_APPROX_LINEAR_COEF * x)
+              .abs())
+            .exp();
   f.array() = (x >= 0).select(1.0, da.array()) / (1.0 + da.array());
   da.array()
-      *= (3.0 * 0.07056 * x.square() + 1.5976) / (1.0 + da.array()).square();
+      *= (3.0 * PHI_APPROX_CUBIC_COEF * x.square() + PHI_APPROX_LINEAR_COEF)
+         / (1.0 + da.array()).square();
   return make_callback_var(f, [a, da](auto& vi) mutable {
     a.adj().array() += vi.adj().array() * da.array();
   });
diff --git a/test/unit/math/mix/fun/Phi_approx_test.cpp b/test/unit/math/mix/fun/Phi_approx_test.cpp
index 1bb709bc9da..447c86b63fd 100644
--- a/test/unit/math/mix/fun/Phi_approx_test.cpp
+++ b/test/unit/math/mix/fun/Phi_approx_test.cpp
@@ -28,22 +28,23 @@ TEST(mathMixMatFun, PhiApprox_varmat) {
 }
 
 namespace Phi_approx_test {
-// d/dx inv_logit(0.07056 x^3 + 1.5976 x) by mpmath at 60 digits; at x = 40
+// d/dx Phi_approx(x) with PHI_APPROX_CUBIC_COEF and PHI_APPROX_LINEAR_COEF, by
+// mpmath at 60 digits; at x = 40
 // and -40 the true value, about 4e-1987, underflows to 0
 inline std::vector<double> args() {
   return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
 }
 inline std::vector<double> derivs() {
   return {0.0,
-          5.9589784505245725e-37,
-          8.7097694481792379e-21,
-          5.1340418352583332e-15,
-          0.24152728992165639,
-          0.3994,
-          0.24152728992165639,
-          5.1340418352583332e-15,
-          8.7097694481792379e-21,
-          5.9589784505245725e-37,
+          5.934767576595731e-37,
+          8.691773531236095e-21,
+          5.126988948527854e-15,
+          0.2415281648469019,
+          0.399400718560196,
+          0.2415281648469019,
+          5.126988948527854e-15,
+          8.691773531236095e-21,
+          5.934767576595731e-37,
           0.0};
 }
 // Relative check, absolute below the smallest normal double.
diff --git a/test/unit/math/opencl/rev/Phi_approx_test.cpp b/test/unit/math/opencl/rev/Phi_approx_test.cpp
index d3e06c16fd5..22378787f73 100644
--- a/test/unit/math/opencl/rev/Phi_approx_test.cpp
+++ b/test/unit/math/opencl/rev/Phi_approx_test.cpp
@@ -24,9 +24,9 @@ TEST(OpenCL_Phi_approx, rev_tail_derivative_against_references) {
   x << -40, -10, -7, 1, 7, 8, 10, 40;
   // at x = -40 and 40 the true value, about 4e-1987, underflows to 0
   Eigen::VectorXd expected(N);
-  expected << 0.0, 5.9589784505245725e-37, 5.1340418352583332e-15,
-      0.24152728992165639, 5.1340418352583332e-15, 8.7097694481792379e-21,
-      5.9589784505245725e-37, 0.0;
+  expected << 0.0, 5.934767576595731e-37, 5.126988948527854e-15,
+      0.2415281648469019, 5.126988948527854e-15, 8.691773531236095e-21,
+      5.934767576595731e-37, 0.0;
 
   stan::math::var_value<stan::math::matrix_cl<double>> x_cl(
       stan::math::to_matrix_cl(x));

From 32e48bd934df5a5e7ffef27b6c9a2f3227c26ea1 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 21:10:06 +0200
Subject: [PATCH 06/10] Move the Phi_approx forward-mode test to the Phi_approx
 branch

---
 .../unit/math/mix/fun/Phi_approx_fwd_test.cpp | 66 -------------------
 1 file changed, 66 deletions(-)
 delete mode 100644 test/unit/math/mix/fun/Phi_approx_fwd_test.cpp

diff --git a/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp b/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp
deleted file mode 100644
index 9c175a1c417..00000000000
--- a/test/unit/math/mix/fun/Phi_approx_fwd_test.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#include <test/unit/math/test_ad.hpp>
-#include <test/unit/math/mix/fun/inv_logit_tail_refs.hpp>
-#include <vector>
-
-// Forward-mode derivatives of Phi_approx(x) = inv_logit(0.07056 x^3 + 1.5976 x)
-// come from the fvar overload of inv_logit; they were 0 wherever its own
-// derivative was.
-namespace Phi_approx_fwd_test {
-inline std::vector<double> args() {
-  return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
-}
-// first derivatives by mpmath at 60 digits; at x = 40 and -40 the true
-// value, about 4e-1987, underflows to 0
-inline std::vector<double> derivs() {
-  return {0.0,
-          5.9589784505245725e-37,
-          8.7097694481792379e-21,
-          5.1340418352583332e-15,
-          0.24152728992165639,
-          0.3994,
-          0.24152728992165639,
-          5.1340418352583332e-15,
-          8.7097694481792379e-21,
-          5.9589784505245725e-37,
-          0.0};
-}
-// second derivatives, computed the same way
-inline std::vector<double> second_derivs() {
-  return {0.0,
-          1.3455155964317485e-35,
-          1.2996275105706113e-19,
-          6.018298087621098e-14,
-          0.24180038682517271,
-          0.0,
-          -0.24180038682517271,
-          -6.018298087621098e-14,
-          -1.2996275105706113e-19,
-          -1.3455155964317485e-35,
-          0.0};
-}
-}  // namespace Phi_approx_fwd_test
-
-TEST(mathMixMatFun, PhiApproxFwdDerivativeTails) {
-  using inv_logit_tail_refs::expect_rel;
-  using stan::math::fvar;
-  using stan::math::Phi_approx;
-  using stan::math::var;
-  std::vector<double> xs = Phi_approx_fwd_test::args();
-  std::vector<double> ds = Phi_approx_fwd_test::derivs();
-  std::vector<double> d2s = Phi_approx_fwd_test::second_derivs();
-  for (size_t i = 0; i < xs.size(); ++i) {
-    expect_rel(ds[i], Phi_approx(fvar<double>(xs[i], 1)).d_, xs[i]);
-
-    fvar<fvar<double>> ff(fvar<double>(xs[i], 1), fvar<double>(1, 0));
-    fvar<fvar<double>> y_ff = Phi_approx(ff);
-    expect_rel(ds[i], y_ff.d_.val_, xs[i]);
-    expect_rel(d2s[i], y_ff.d_.d_, xs[i]);
-
-    fvar<var> fv(xs[i], 1);
-    fvar<var> y_fv = Phi_approx(fv);
-    expect_rel(ds[i], y_fv.d_.val(), xs[i]);
-    y_fv.d_.grad();
-    expect_rel(d2s[i], fv.val_.adj(), xs[i]);
-    stan::math::recover_memory();
-  }
-}

From 5b4e614f284acadd4cb87b0f99eb738eba93654b Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 21:24:44 +0200
Subject: [PATCH 07/10] Tests: Phi_approx tail references shared by the mix and
 OpenCL tests, forward-mode derivatives

---
 .../math/mix/fun/Phi_approx_tail_refs.hpp     | 58 ++++++++++++++
 test/unit/math/mix/fun/Phi_approx_test.cpp    | 77 +++++++++----------
 test/unit/math/opencl/rev/Phi_approx_test.cpp | 21 ++---
 3 files changed, 102 insertions(+), 54 deletions(-)
 create mode 100644 test/unit/math/mix/fun/Phi_approx_tail_refs.hpp

diff --git a/test/unit/math/mix/fun/Phi_approx_tail_refs.hpp b/test/unit/math/mix/fun/Phi_approx_tail_refs.hpp
new file mode 100644
index 00000000000..963fbe1f6f5
--- /dev/null
+++ b/test/unit/math/mix/fun/Phi_approx_tail_refs.hpp
@@ -0,0 +1,58 @@
+#ifndef TEST_UNIT_MATH_MIX_FUN_PHI_APPROX_TAIL_REFS_HPP
+#define TEST_UNIT_MATH_MIX_FUN_PHI_APPROX_TAIL_REFS_HPP
+
+#include <gtest/gtest.h>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+// First and second derivatives of Phi_approx(x) = inv_logit(c3 x^3 + c1 x),
+// computed with mpmath at 60 digits from the double values of
+// PHI_APPROX_CUBIC_COEF (c3) and PHI_APPROX_LINEAR_COEF (c1). At x = 40 and
+// -40 the true values, about 3e-1987 and 1e-1984, underflow to 0.
+namespace Phi_approx_tail_refs {
+
+inline std::vector<double> args() {
+  return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
+}
+
+inline std::vector<double> derivs() {
+  return {0.0,
+          5.934767576595741e-37,
+          8.691773531236104e-21,
+          5.126988948527858e-15,
+          0.2415281648469019,
+          0.399400718560196,
+          0.2415281648469019,
+          5.126988948527858e-15,
+          8.691773531236104e-21,
+          5.934767576595741e-37,
+          0.0};
+}
+
+inline std::vector<double> second_derivs() {
+  return {0.0,
+          1.3401219329217354e-35,
+          1.297010750153944e-19,
+          6.010339719762636e-14,
+          0.24180176878701165,
+          0.0,
+          -0.24180176878701165,
+          -6.010339719762636e-14,
+          -1.297010750153944e-19,
+          -1.3401219329217354e-35,
+          0.0};
+}
+
+// Relative check, absolute below the smallest normal double.
+inline void expect_rel(double expected, double actual, double x,
+                       double tol = 1e-13) {
+  EXPECT_NEAR(
+      expected, actual,
+      tol * std::fmax(std::fabs(expected), std::numeric_limits<double>::min()))
+      << "x = " << x;
+}
+
+}  // namespace Phi_approx_tail_refs
+
+#endif
diff --git a/test/unit/math/mix/fun/Phi_approx_test.cpp b/test/unit/math/mix/fun/Phi_approx_test.cpp
index 447c86b63fd..1d7117a54f6 100644
--- a/test/unit/math/mix/fun/Phi_approx_test.cpp
+++ b/test/unit/math/mix/fun/Phi_approx_test.cpp
@@ -1,6 +1,5 @@
 #include <test/unit/math/test_ad.hpp>
-#include <cmath>
-#include <limits>
+#include <test/unit/math/mix/fun/Phi_approx_tail_refs.hpp>
 #include <vector>
 
 TEST(mathMixMatFun, PhiApprox) {
@@ -27,74 +26,72 @@ TEST(mathMixMatFun, PhiApprox_varmat) {
   expect_ad_vector_matvar(f, A);
 }
 
-namespace Phi_approx_test {
-// d/dx Phi_approx(x) with PHI_APPROX_CUBIC_COEF and PHI_APPROX_LINEAR_COEF, by
-// mpmath at 60 digits; at x = 40
-// and -40 the true value, about 4e-1987, underflows to 0
-inline std::vector<double> args() {
-  return {-40, -10, -8, -7, -1, 0, 1, 7, 8, 10, 40};
-}
-inline std::vector<double> derivs() {
-  return {0.0,
-          5.934767576595731e-37,
-          8.691773531236095e-21,
-          5.126988948527854e-15,
-          0.2415281648469019,
-          0.399400718560196,
-          0.2415281648469019,
-          5.126988948527854e-15,
-          8.691773531236095e-21,
-          5.934767576595731e-37,
-          0.0};
-}
-// Relative check, absolute below the smallest normal double.
-inline void expect_rel(double expected, double actual, double x,
-                       double tol = 1e-13) {
-  EXPECT_NEAR(
-      expected, actual,
-      tol * std::fmax(std::fabs(expected), std::numeric_limits<double>::min()))
-      << "x = " << x;
-}
-}  // namespace Phi_approx_test
-
 TEST(mathMixMatFun, PhiApproxScalarDerivativeTails) {
   using stan::math::var;
-  std::vector<double> xs = Phi_approx_test::args();
-  std::vector<double> ds = Phi_approx_test::derivs();
+  std::vector<double> xs = Phi_approx_tail_refs::args();
+  std::vector<double> ds = Phi_approx_tail_refs::derivs();
   for (size_t i = 0; i < xs.size(); ++i) {
     var a = xs[i];
     stan::math::Phi_approx(a).grad();
-    Phi_approx_test::expect_rel(ds[i], a.adj(), xs[i]);
+    Phi_approx_tail_refs::expect_rel(ds[i], a.adj(), xs[i]);
     stan::math::recover_memory();
   }
 }
 
 TEST(mathMixMatFun, PhiApproxDerivativeTails) {
+  using Phi_approx_tail_refs::expect_rel;
   using stan::math::Phi_approx;
   using stan::math::sum;
   using stan::math::var;
-  std::vector<double> xs = Phi_approx_test::args();
-  std::vector<double> ds = Phi_approx_test::derivs();
+  std::vector<double> xs = Phi_approx_tail_refs::args();
+  std::vector<double> ds = Phi_approx_tail_refs::derivs();
   Eigen::VectorXd x = stan::math::to_vector(xs);
 
   stan::math::var_value<Eigen::VectorXd> xv(x);
   sum(Phi_approx(xv)).grad();
   for (int i = 0; i < x.size(); ++i) {
-    Phi_approx_test::expect_rel(ds[i], xv.adj()(i), xs[i]);
+    expect_rel(ds[i], xv.adj()(i), xs[i]);
   }
   stan::math::recover_memory();
 
   Eigen::Matrix<var, Eigen::Dynamic, 1> xm = x;
   sum(Phi_approx(xm)).grad();
   for (int i = 0; i < x.size(); ++i) {
-    Phi_approx_test::expect_rel(ds[i], xm(i).adj(), xs[i]);
+    expect_rel(ds[i], xm(i).adj(), xs[i]);
   }
   stan::math::recover_memory();
 
   std::vector<var> xs_v(xs.begin(), xs.end());
   sum(Phi_approx(xs_v)).grad();
   for (size_t i = 0; i < xs.size(); ++i) {
-    Phi_approx_test::expect_rel(ds[i], xs_v[i].adj(), xs[i]);
+    expect_rel(ds[i], xs_v[i].adj(), xs[i]);
   }
   stan::math::recover_memory();
 }
+
+// Forward mode differentiates through the fvar overload of inv_logit, whose
+// derivative was 0 wherever the reverse-mode one was.
+TEST(mathMixMatFun, PhiApproxFwdDerivativeTails) {
+  using Phi_approx_tail_refs::expect_rel;
+  using stan::math::fvar;
+  using stan::math::Phi_approx;
+  using stan::math::var;
+  std::vector<double> xs = Phi_approx_tail_refs::args();
+  std::vector<double> ds = Phi_approx_tail_refs::derivs();
+  std::vector<double> d2s = Phi_approx_tail_refs::second_derivs();
+  for (size_t i = 0; i < xs.size(); ++i) {
+    expect_rel(ds[i], Phi_approx(fvar<double>(xs[i], 1)).d_, xs[i]);
+
+    fvar<fvar<double>> ff(fvar<double>(xs[i], 1), fvar<double>(1, 0));
+    fvar<fvar<double>> y_ff = Phi_approx(ff);
+    expect_rel(ds[i], y_ff.d_.val_, xs[i]);
+    expect_rel(d2s[i], y_ff.d_.d_, xs[i]);
+
+    fvar<var> fv(xs[i], 1);
+    fvar<var> y_fv = Phi_approx(fv);
+    expect_rel(ds[i], y_fv.d_.val(), xs[i]);
+    y_fv.d_.grad();
+    expect_rel(d2s[i], fv.val_.adj(), xs[i]);
+    stan::math::recover_memory();
+  }
+}
diff --git a/test/unit/math/opencl/rev/Phi_approx_test.cpp b/test/unit/math/opencl/rev/Phi_approx_test.cpp
index 22378787f73..53cbd886d42 100644
--- a/test/unit/math/opencl/rev/Phi_approx_test.cpp
+++ b/test/unit/math/opencl/rev/Phi_approx_test.cpp
@@ -3,6 +3,8 @@
 #include <stan/math.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/opencl/util.hpp>
+#include <test/unit/math/mix/fun/Phi_approx_tail_refs.hpp>
+#include <vector>
 
 auto Phi_approx_functor
     = [](const auto& a) { return stan::math::Phi_approx(a); };
@@ -16,28 +18,19 @@ TEST(OpenCL_Phi_approx, prim_rev_values_small) {
 /**
  * Device derivative in the tails against fixed references, which
  * `compare_cpu_opencl_prim_rev` cannot check: its relative tolerance floors
- * the denominator at 1. References by mpmath at 60 digits.
+ * the denominator at 1. References in Phi_approx_tail_refs.hpp.
  */
 TEST(OpenCL_Phi_approx, rev_tail_derivative_against_references) {
-  const int N = 8;
-  Eigen::VectorXd x(N);
-  x << -40, -10, -7, 1, 7, 8, 10, 40;
-  // at x = -40 and 40 the true value, about 4e-1987, underflows to 0
-  Eigen::VectorXd expected(N);
-  expected << 0.0, 5.934767576595731e-37, 5.126988948527854e-15,
-      0.2415281648469019, 5.126988948527854e-15, 8.691773531236095e-21,
-      5.934767576595731e-37, 0.0;
+  const Eigen::VectorXd x = stan::math::to_vector(Phi_approx_tail_refs::args());
+  const std::vector<double> expected = Phi_approx_tail_refs::derivs();
 
   stan::math::var_value<stan::math::matrix_cl<double>> x_cl(
       stan::math::to_matrix_cl(x));
   stan::math::sum(stan::math::Phi_approx(x_cl)).grad();
 
   const Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj());
-  for (int i = 0; i < N; ++i) {
-    EXPECT_NEAR(
-        expected[i], adj[i],
-        1e-13 * std::fmax(expected[i], std::numeric_limits<double>::min()))
-        << "x = " << x[i];
+  for (int i = 0; i < x.size(); ++i) {
+    Phi_approx_tail_refs::expect_rel(expected[i], adj[i], x[i]);
   }
   stan::math::recover_memory();
 }

From d2742ffbaff268dad17585bc296801864d225d17 Mon Sep 17 00:00:00 2001
From: "Jachym.Barvinek" <jachym.barvinek@alza.cz>
Date: Mon, 12 Oct 2026 22:29:25 +0200
Subject: [PATCH 08/10] Phi_approx: OpenCL value path as a kernel-generator
 expression of the constants

---
 .../kernel_generator/elt_function_cl.hpp      |  4 --
 .../kernels/device_functions/Phi_approx.hpp   | 38 -------------------
 stan/math/opencl/prim.hpp                     |  1 +
 stan/math/opencl/prim/Phi_approx.hpp          | 38 +++++++++++++++++++
 stan/math/opencl/rev/Phi_approx.hpp           |  1 +
 stan/math/prim/fun/constants.hpp              |  3 +-
 6 files changed, 41 insertions(+), 44 deletions(-)
 delete mode 100644 stan/math/opencl/kernels/device_functions/Phi_approx.hpp
 create mode 100644 stan/math/opencl/prim/Phi_approx.hpp

diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp
index df349a7b483..e7f006eb50d 100644
--- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp
+++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp
@@ -24,7 +24,6 @@
 #include <stan/math/opencl/kernels/device_functions/logit.hpp>
 #include <stan/math/opencl/kernels/device_functions/multiply_log.hpp>
 #include <stan/math/opencl/kernels/device_functions/Phi.hpp>
-#include <stan/math/opencl/kernels/device_functions/Phi_approx.hpp>
 #include <stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp>
 #include <stan/math/opencl/kernels/device_functions/trigamma.hpp>
 #include <stan/math/opencl/matrix_cl_view.hpp>
@@ -314,9 +313,6 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit,
 ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::log1m_device_function,
                                  opencl_kernels::logit_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function)
-ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx,
-                                 opencl_kernels::inv_logit_device_function,
-                                 opencl_kernels::phi_approx_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(
     std_normal_lcdf_scaled_impl,
     opencl_kernels::std_normal_lcdf_device_function)
diff --git a/stan/math/opencl/kernels/device_functions/Phi_approx.hpp b/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
deleted file mode 100644
index cb635ebf482..00000000000
--- a/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_PHI_APPROX_HPP
-#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_PHI_APPROX_HPP
-#ifdef STAN_OPENCL
-
-#include <stan/math/opencl/stringify.hpp>
-#include <string>
-
-namespace stan {
-namespace math {
-namespace opencl_kernels {
-// \cond
-static constexpr const char* phi_approx_device_function
-    = "\n"
-      "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_PHI_APPROX\n"
-      "#define "
-      "STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_PHI_APPROX\n" STRINGIFY(
-          // \endcond
-          /** \ingroup opencl_kernels
-           * Return the Phi_approx function applied to the specified
-           * argument. The coefficients are PHI_APPROX_CUBIC_COEF and
-           * PHI_APPROX_LINEAR_COEF of stan/math/prim/fun/constants.hpp,
-           * repeated here because a kernel source string cannot use them.
-           *
-           * @return Phi_approx(x)
-           */
-          inline double Phi_approx(double x) {
-            return inv_logit(0.07056409656566716 * pow(x, 3.0)
-                             + 1.597602874240784 * x);
-          }
-          // \cond
-          ) "\n#endif\n";  // NOLINT
-// \endcond
-}  // namespace opencl_kernels
-}  // namespace math
-}  // namespace stan
-
-#endif
-#endif
diff --git a/stan/math/opencl/prim.hpp b/stan/math/opencl/prim.hpp
index 4ae8ca1649c..11f6353d250 100644
--- a/stan/math/opencl/prim.hpp
+++ b/stan/math/opencl/prim.hpp
@@ -101,6 +101,7 @@
 
 #include <stan/math/opencl/prim_constraint.hpp>
 
+#include <stan/math/opencl/prim/Phi_approx.hpp>
 #include <stan/math/opencl/prim/add_diag.hpp>
 #include <stan/math/opencl/prim/append_array.hpp>
 #include <stan/math/opencl/prim/bernoulli_cdf.hpp>
diff --git a/stan/math/opencl/prim/Phi_approx.hpp b/stan/math/opencl/prim/Phi_approx.hpp
new file mode 100644
index 00000000000..7d80972e131
--- /dev/null
+++ b/stan/math/opencl/prim/Phi_approx.hpp
@@ -0,0 +1,38 @@
+#ifndef STAN_MATH_OPENCL_PRIM_PHI_APPROX_HPP
+#define STAN_MATH_OPENCL_PRIM_PHI_APPROX_HPP
+#ifdef STAN_OPENCL
+
+#include <stan/math/opencl/kernel_generator.hpp>
+#include <stan/math/opencl/ref_type.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the elementwise logistic approximation to the standard normal CDF,
+ * `inv_logit(PHI_APPROX_CUBIC_COEF * x^3 + PHI_APPROX_LINEAR_COEF * x)`, of a
+ * kernel generator expression.
+ *
+ * @tparam T type of the kernel generator expression
+ * @param x argument
+ * @return `Phi_approx` of each element of x
+ */
+template <typename T,
+          require_all_kernel_expressions_and_none_scalar_t<T>* = nullptr>
+inline auto Phi_approx(T&& x) {
+  return make_holder_cl(
+      [](auto&& a) {
+        return inv_logit(PHI_APPROX_CUBIC_COEF * pow(a, 3.0)
+                         + PHI_APPROX_LINEAR_COEF * a);
+      },
+      to_ref(std::forward<T>(x)));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
+#endif
diff --git a/stan/math/opencl/rev/Phi_approx.hpp b/stan/math/opencl/rev/Phi_approx.hpp
index 69870016d0d..ef17e5d7697 100644
--- a/stan/math/opencl/rev/Phi_approx.hpp
+++ b/stan/math/opencl/rev/Phi_approx.hpp
@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/opencl/kernel_generator.hpp>
+#include <stan/math/opencl/prim/Phi_approx.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
 #include <stan/math/prim/fun/constants.hpp>
diff --git a/stan/math/prim/fun/constants.hpp b/stan/math/prim/fun/constants.hpp
index 38e0684430f..59d7c3724d5 100644
--- a/stan/math/prim/fun/constants.hpp
+++ b/stan/math/prim/fun/constants.hpp
@@ -199,8 +199,7 @@ const double POISSON_MAX_RATE = std::pow(2.0, 30);
  * maximum absolute deviation from \f$ \Phi \f$ over the real line, rounded
  * to double precision: the deviation is 1.40313e-4, attained at
  * \f$ x = \pm 0.4816, \pm 1.4764, \pm 2.6876 \f$. The paper prints the pair
- * rounded to 0.07056 and 1.5976 (maximum deviation 1.41448e-4). The OpenCL
- * device function repeats the values as literals.
+ * rounded to 0.07056 and 1.5976 (maximum deviation 1.41448e-4).
  */
 static constexpr double PHI_APPROX_CUBIC_COEF = 0.07056409656566716;
 static constexpr double PHI_APPROX_LINEAR_COEF = 1.597602874240784;

From 73d1a2a45ddb92270a04ea1af0bceaa1e3494ab2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= <jachymb@gmail.com>
Date: Tue, 20 Oct 2026 01:07:58 +0200
Subject: [PATCH 09/10] Compute inv_logit of Matrix<var> in one scalar pass

The Eigen::Matrix<var> overload computed exp(-|x|) with a vectorized pass
over a heap copy of the values, stored it, and divided again in the reverse
pass: 1.07-1.30 of develop on glibc. Now one loop computes the value and the
derivative e / (1 + e)^2 per element and stores the derivative, as suggested
in the review of #3446: 0.94-1.11 on glibc and 1.00-1.07 on MinGW, within the
noise of identical code; values within 2 ulp of develop, +8 arena bytes per
element as before.
---
 stan/math/rev/fun/inv_logit.hpp | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp
index e739679132e..ed2f744d427 100644
--- a/stan/math/rev/fun/inv_logit.hpp
+++ b/stan/math/rev/fun/inv_logit.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/functor/apply_scalar_binary.hpp>
+#include <cmath>
 
 namespace stan {
 namespace math {
@@ -50,7 +51,8 @@ inline auto inv_logit(T&& a) {
  * The derivative of inverse logit is
  *
  * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \frac{e}{(1 + e)^2}\f$
- * with \f$e = \exp(-|x|)\f$.
+ * with \f$e = \exp(-|x|)\f$. It is computed in the same pass as the value
+ * and stored for the reverse pass.
  *
  * @tparam T type of Eigen expression
  * @param x Eigen expression
@@ -59,11 +61,17 @@ inline auto inv_logit(T&& a) {
 template <typename T, require_eigen_vt<is_var, T>* = nullptr>
 inline auto inv_logit(T&& x) {
   auto x_arena = to_arena(std::forward<T>(x));
-  const auto x_val = x_arena.val().array().eval();
-  arena_t<plain_type_t<decltype(x_val)>> e = (-x_val.abs()).exp();
-  arena_t<T> ret = ((x_val >= 0).select(1.0, e) / (1.0 + e)).matrix();
-  reverse_pass_callback([x_arena, ret, e]() mutable {
-    x_arena.adj().array() += ret.adj().array() * e / (1.0 + e).square();
+  arena_t<T> ret(x_arena.rows(), x_arena.cols());
+  arena_t<plain_type_t<decltype(x_arena.val())>> d(x_arena.rows(),
+                                                   x_arena.cols());
+  for (Eigen::Index i = 0; i < x_arena.size(); ++i) {
+    const double xi = x_arena.coeff(i).val();
+    const double e = std::exp(-std::fabs(xi));
+    ret.coeffRef(i) = (xi >= 0 ? 1.0 : e) / (1.0 + e);
+    d.coeffRef(i) = e / ((1.0 + e) * (1.0 + e));
+  }
+  reverse_pass_callback([x_arena, ret, d]() mutable {
+    x_arena.adj().array() += ret.adj().array() * d.array();
   });
   return ret;
 }

From de0c943cbb8f0d8c1abcac9b3c16152a84c79b3b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= <jachymb@gmail.com>
Date: Tue, 20 Oct 2026 09:30:39 +0200
Subject: [PATCH 10/10] Re-run CI: Windows runner lost communication in mix
 tests 5