From 2f40496bfbcb78f63b6d74a1c5d339c12d5632ea Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Fri, 23 Oct 2026 09:55:21 +0300
Subject: [PATCH] Compute the cube in Phi_approx by multiplication instead of
 pow

---
 stan/math/fwd/fun/Phi_approx.hpp                   |  4 ++--
 .../opencl/kernels/device_functions/Phi_approx.hpp |  3 ++-
 stan/math/prim/fun/Phi_approx.hpp                  |  4 ++--
 test/unit/math/mix/fun/Phi_approx_test.cpp         | 14 ++++++++++++++
 4 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/stan/math/fwd/fun/Phi_approx.hpp b/stan/math/fwd/fun/Phi_approx.hpp
index 2509a44348a..59f20f0dd74 100644
--- a/stan/math/fwd/fun/Phi_approx.hpp
+++ b/stan/math/fwd/fun/Phi_approx.hpp
@@ -4,7 +4,6 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/fun/inv_logit.hpp>
-#include <stan/math/fwd/fun/pow.hpp>
 #include <stan/math/prim/fun/Phi_approx.hpp>
 #include <cmath>
 
@@ -23,7 +22,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);
+  const fvar<T> x_squared = x * x;
+  return inv_logit(0.07056 * x * x_squared + 1.5976 * 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..aff105302d4 100644
--- a/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
+++ b/stan/math/opencl/kernels/device_functions/Phi_approx.hpp
@@ -22,7 +22,8 @@ static constexpr const char* phi_approx_device_function
            * @return Phi_approx(x)
            */
           inline double Phi_approx(double x) {
-            return inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
+            const double x_squared = x * x;
+            return inv_logit(0.07056 * x * x_squared + 1.5976 * x);
           }
           // \cond
           ) "\n#endif\n";  // NOLINT
diff --git a/stan/math/prim/fun/Phi_approx.hpp b/stan/math/prim/fun/Phi_approx.hpp
index 8351cb13e49..fda884411c7 100644
--- a/stan/math/prim/fun/Phi_approx.hpp
+++ b/stan/math/prim/fun/Phi_approx.hpp
@@ -21,8 +21,8 @@ namespace math {
  * @return Probability random sample is less than or equal to argument.
  */
 inline double Phi_approx(double x) {
-  using std::pow;
-  return inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
+  const double x_squared = x * x;
+  return inv_logit(0.07056 * x * x_squared + 1.5976 * x);
 }
 
 /**
diff --git a/test/unit/math/mix/fun/Phi_approx_test.cpp b/test/unit/math/mix/fun/Phi_approx_test.cpp
index b89df1935b7..0843656c189 100644
--- a/test/unit/math/mix/fun/Phi_approx_test.cpp
+++ b/test/unit/math/mix/fun/Phi_approx_test.cpp
@@ -6,6 +6,20 @@ TEST(mathMixMatFun, PhiApprox) {
   stan::test::expect_unary_vectorized(f, -3.0, 1, 1.3, 3);
 }
 
+TEST(mathMixMatFun, PhiApprox_all_modes_give_the_same_value) {
+  using stan::math::fvar;
+  using stan::math::Phi_approx;
+  using stan::math::var;
+  for (int i = -2500; i <= 1000; ++i) {
+    const double x = 0.01 * i;
+    const double expected = Phi_approx(x);
+    EXPECT_EQ(expected, Phi_approx(var(x)).val()) << "x = " << x;
+    EXPECT_EQ(expected, Phi_approx(fvar<double>(x, 1)).val()) << "x = " << x;
+    EXPECT_EQ(expected, Phi_approx(fvar<var>(x, 1)).val().val()) << "x = " << x;
+  }
+  stan::math::recover_memory();
+}
+
 TEST(mathMixMatFun, PhiApprox_varmat) {
   using stan::math::vec_concat;
   using stan::test::expect_ad_vector_matvar;