From 056dd03fd15628a36f1ca96e7e65ed97a944a238 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Tue, 15 Sep 2026 15:56:39 +0300 Subject: [PATCH 01/24] Add erfcx scaled complementary error function --- stan/math/fwd/fun.hpp | 1 + stan/math/fwd/fun/erfcx.hpp | 31 ++++ .../kernel_generator/elt_function_cl.hpp | 2 + .../opencl/kernels/device_functions/erfcx.hpp | 76 +++++++++ stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/erfcx.hpp | 152 ++++++++++++++++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/erfcx.hpp | 70 ++++++++ test/unit/math/fwd/fun/erfcx_test.cpp | 47 ++++++ test/unit/math/mix/fun/erfcx_test.cpp | 22 +++ test/unit/math/prim/fun/erfcx_test.cpp | 110 +++++++++++++ test/unit/math/rev/fun/erfcx_test.cpp | 83 ++++++++++ 12 files changed, 596 insertions(+) create mode 100644 stan/math/fwd/fun/erfcx.hpp create mode 100644 stan/math/opencl/kernels/device_functions/erfcx.hpp create mode 100644 stan/math/prim/fun/erfcx.hpp create mode 100644 stan/math/rev/fun/erfcx.hpp create mode 100644 test/unit/math/fwd/fun/erfcx_test.cpp create mode 100644 test/unit/math/mix/fun/erfcx_test.cpp create mode 100644 test/unit/math/prim/fun/erfcx_test.cpp create mode 100644 test/unit/math/rev/fun/erfcx_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index 21c6d70eee0..ffd56a9297e 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/erfcx.hpp b/stan/math/fwd/fun/erfcx.hpp new file mode 100644 index 00000000000..4e1328e2d52 --- /dev/null +++ b/stan/math/fwd/fun/erfcx.hpp @@ -0,0 +1,31 @@ +#ifndef STAN_MATH_FWD_FUN_ERFCX_HPP +#define STAN_MATH_FWD_FUN_ERFCX_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the scaled complementary error function of the argument. + * + * The derivative `2 * x * erfcx(x) - 2 / sqrt(pi)` reuses the value, so the + * tangent is as accurate as the value itself. + * + * @tparam T inner type of the fvar + * @param x argument + * @return scaled complementary error function of the argument + */ +template +inline fvar erfcx(const fvar& x) { + T v = erfcx(x.val_); + return fvar(v, x.d_ * (2.0 * x.val_ * v - TWO_OVER_SQRT_PI)); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 96bbabbaf65..6ca3be497e6 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -297,6 +298,7 @@ ADD_UNARY_FUNCTION_PASS_ZERO(trunc) ADD_UNARY_FUNCTION_WITH_INCLUDES(digamma, opencl_kernels::digamma_device_function) +ADD_UNARY_FUNCTION_WITH_INCLUDES(erfcx, opencl_kernels::erfcx_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(log1m, opencl_kernels::log1m_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(log_inv_logit, opencl_kernels::log1p_exp_device_function, diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp new file mode 100644 index 00000000000..032d8c70eaf --- /dev/null +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -0,0 +1,76 @@ +#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX_HPP +#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX_HPP +#ifdef STAN_OPENCL + +#include +#include + +namespace stan { +namespace math { +namespace opencl_kernels { + +// \cond +static constexpr const char* erfcx_device_function + = "\n" + "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX\n" + "#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX\n" STRINGIFY( + // \endcond + /** \ingroup opencl_kernels + * + * Return the scaled complementary error function + * exp(x * x) * erfc(x) of the kernel generator expression. + * + * Mirrors stan/math/prim/fun/erfcx.hpp: the Cody (1969) rational + * approximation above 4, and below it exp(x * x) * erfc(x) with + * x * x Dekker-split so that exp() receives an exact argument. + * Keep the two in step; the split is what holds the lower branch + * to a few ulp instead of x * x * eps. + * + * @param x argument + * @return scaled complementary error function of the argument + */ + double erfcx(double x) { + if (x >= 4.0) { + double u = 1.0 / (x * x); + double p + = 0.000658749161529837803157 + + u + * (0.0160837851487422766278 + + u + * (0.125781726111229246204 + + u + * (0.360344899949804439429 + + u + * (0.305326634961232344035 + + u * 0.0163153871373020978498)))); + double q + = -0.00233520497626869185443 + + u + * (-0.0605183413124413191178 + + u + * (-0.527905102951428412248 + + u + * (-1.87295284992346047209 + + u + * (-2.56852019228982242072 + + u * -1.0)))); + return (M_2_SQRTPI * 0.5 + (p / q) * u) / x; + } + if (x < -27.0) { + return INFINITY; + } + double t = 134217729.0 * x; + double x_hi = t - (t - x); + double x_lo = x - x_hi; + return exp(x_hi * x_hi) * exp(x_lo * (x + x_hi)) * erfc(x); + } + // \cond + ) "\n#endif\n"; // NOLINT +// \endcond + +} // namespace opencl_kernels +} // namespace math +} // namespace stan + +#endif +#endif diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 8e0b01d2093..19d850587f0 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -83,6 +83,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp new file mode 100644 index 00000000000..006562f1f8b --- /dev/null +++ b/stan/math/prim/fun/erfcx.hpp @@ -0,0 +1,152 @@ +#ifndef STAN_MATH_PRIM_FUN_ERFCX_HPP +#define STAN_MATH_PRIM_FUN_ERFCX_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the scaled complementary error function of the argument. + * + * \f$\mbox{erfcx}(x) = \exp(x^2)\,\mbox{erfc}(x)\f$ + * + * The scaling removes the Gaussian factor, so `erfcx` stays in a narrow + * range (it decays like \f$1/(x\sqrt{\pi})\f$ for large positive `x`) where + * the two factors individually overflow and underflow: `erfc(x)` underflows + * to zero above `x = 27` while `exp(x*x)` overflows above `x = 26.6`, so + * their product cannot be formed directly in the upper tail at all. It is + * the natural building block for normal tail quantities, e.g. the standard + * normal log CDF is + * `LOG_HALF + log(erfcx(-x * INV_SQRT_TWO)) - x * x / 2` and the inverse + * Mills ratio is `SQRT_TWO_OVER_SQRT_PI / erfcx(-x * INV_SQRT_TWO)`, both + * without cancellation. + * + * Two branches, with the crossover and the accuracy below established by + * measurement against a `long double` reference over 40000 points per + * interval: + * + * - `x >= 4`: the Cody (1969) rational approximation, which gives `erfcx` + * directly with no `exp` and no `erfc` call. Measured at 2.0 to 2.5 ulp + * over `[4, 50]` and 5.4 ulp out to 200, and about 2.4 times faster than + * forming the product. + * - `x < 4`: `exp(x * x) * erfc(x)`, with `x * x` split so that the + * argument handed to `exp` is exact. Without the split the rounding of + * `x * x` is amplified by `exp` into a relative error of about + * `x * x * eps`, which reaches 512 ulp at `x = -26`. Splitting holds this + * branch to under 5 ulp. + * + * W. J. Cody, Math. Comp. 23(107):631-638 (1969). + * https://doi.org/10.1090/S0025-5718-1969-0247736-4 + * + * The crossover at 4 is the same one R's `pnorm` uses: it switches to the + * Cody tail form at `|y| > M_SQRT_32`, and since the argument here is the + * `erfc` argument `y / sqrt(2)`, that is `sqrt(32) / sqrt(2) = 4` exactly. + * Evaluated at 60 significant digits, this coefficient set holds to 8.6e-17 + * relative down to 4 and degrades past about 3.5, so 4 sits just inside its + * range. + * + * Measured worst case over the whole range is 5.4 ulp. + * + * The derivative is + * \f$\frac{d}{dx}\mbox{erfcx}(x) = 2x\,\mbox{erfcx}(x) - + \frac{2}{\sqrt{\pi}}\f$. + * + \f[ + \mbox{erfcx}(x) = + \begin{cases} + \exp(x^2)\operatorname{erfc}(x) & \mbox{if } -\infty\leq x \leq \infty + \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + * + * @tparam T An arithmetic type. + * @param xx argument + * @return scaled complementary error function of the argument + */ +template * = nullptr> +inline double erfcx(T&& xx) { + const double x = static_cast(xx); + // Cody (1969). x * x is infinite for x large enough, which correctly + // collapses this to the leading term INV_SQRT_PI / x and, at infinity, + // to zero. + constexpr double cody_min = 4.0; + if (x >= cody_min) { + const double u = 1.0 / (x * x); + const double p + = 0.000658749161529837803157 + + u + * (0.0160837851487422766278 + + u + * (0.125781726111229246204 + + u + * (0.360344899949804439429 + + u + * (0.305326634961232344035 + + u * 0.0163153871373020978498)))); + const double q = -0.00233520497626869185443 + + u + * (-0.0605183413124413191178 + + u + * (-0.527905102951428412248 + + u + * (-1.87295284992346047209 + + u + * (-2.56852019228982242072 + + u * -1.0)))); + return (INV_SQRT_PI + (p / q) * u) / x; + } + // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves + // the binary64 range below -26.63. Returning here also keeps the split + // below from forming inf - inf on -infinity. + constexpr double overflow_max = -27.0; + if (x < overflow_max) { + return INFTY; + } + // Dekker split: x_hi holds the leading 26 bits of x, so x_hi * x_hi is + // exact and exp() is handed an argument carrying no rounding error of + // its own. x_lo * (x + x_hi) is the remainder x * x - x_hi * x_hi. + constexpr double split = 134217729.0; // 2^27 + 1 + const double t = split * x; + const double x_hi = t - (t - x); + const double x_lo = x - x_hi; + return std::exp(x_hi * x_hi) * std::exp(x_lo * (x + x_hi)) * std::erfc(x); +} + +/** + * Structure to wrap `erfcx()` so that it can be vectorized. + * + * @tparam T type of variable + * @param x variable + * @return Scaled complementary error function applied to x. + */ +struct erfcx_fun { + template + static inline auto fun(T&& x) { + return erfcx(std::forward(x)); + } +}; + +/** + * Returns the elementwise `erfcx()` of the input, + * which may be a scalar or any Stan container of numeric scalars. + * + * @tparam T type of container + * @param x container + * @return Scaled complementary error function applied to each value in x. + */ +template < + typename T, + require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr, + require_container_t* = nullptr, require_not_var_matrix_t* = nullptr> +inline auto erfcx(T&& x) { + return apply_scalar_unary::apply(std::forward(x)); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index b9a8bfa3e2c..baaf5a8f0a3 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -51,6 +51,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/erfcx.hpp b/stan/math/rev/fun/erfcx.hpp new file mode 100644 index 00000000000..5228be7af19 --- /dev/null +++ b/stan/math/rev/fun/erfcx.hpp @@ -0,0 +1,70 @@ +#ifndef STAN_MATH_REV_FUN_ERFCX_HPP +#define STAN_MATH_REV_FUN_ERFCX_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * The scaled complementary error function for variables. + * + * The derivative is + * + * \f$\frac{d}{dx}\mbox{erfcx}(x) = 2x\,\mbox{erfcx}(x) - + \frac{2}{\sqrt{\pi}}\f$ + * + * which reuses the function value, so no extra `exp` or `erfc` evaluation is + * needed and the derivative inherits the value's accuracy in the tails. + * + \f[ + \mbox{erfcx}(x) = + \begin{cases} + \exp(x^2)\operatorname{erfc}(x) & \mbox{if } -\infty\leq x \leq \infty + \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN} + \end{cases} + \f] + + \f[ + \frac{\partial\, \mbox{erfcx}(x)}{\partial x} = + \begin{cases} + 2x\operatorname{erfcx}(x) - \frac{2}{\sqrt{\pi}} & \mbox{if } + -\infty\leq x\leq \infty \\[6pt] \textrm{NaN} & \mbox{if } + x = \textrm{NaN} + \end{cases} + \f] + * + * @param a The variable. + * @return Scaled complementary error function applied to the variable. + */ +inline var erfcx(const var& a) { + double val = erfcx(a.val()); + return make_callback_var(val, [a, val](auto& vi) mutable { + a.adj() += vi.adj() * (2.0 * a.val() * val - TWO_OVER_SQRT_PI); + }); +} + +/** + * The scaled complementary error function for matrix variables. + * + * @tparam T a matrix type + * @param a The variable. + * @return Scaled complementary error function applied elementwise. + */ +template * = nullptr> +inline auto erfcx(const var_value& a) { + auto val = to_arena(erfcx(a.val())); + return make_callback_var(val, [a, val](auto& vi) mutable { + a.adj().array() + += vi.adj().array() + * (2.0 * a.val().array() * val.array() - TWO_OVER_SQRT_PI); + }); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/fwd/fun/erfcx_test.cpp b/test/unit/math/fwd/fun/erfcx_test.cpp new file mode 100644 index 00000000000..5ddc8a0dd31 --- /dev/null +++ b/test/unit/math/fwd/fun/erfcx_test.cpp @@ -0,0 +1,47 @@ +#include +#include +#include + +TEST(AgradFwdErfcx, Fvar) { + using stan::math::erfcx; + using stan::math::fvar; + + for (double x_value : {-2.0, -1.0, 0.0, 0.5, 2.0, 12.0}) { + fvar x(x_value, 1.0); + fvar y = erfcx(x); + EXPECT_FLOAT_EQ(erfcx(x_value), y.val_); + EXPECT_FLOAT_EQ( + 2.0 * x_value * erfcx(x_value) - stan::math::TWO_OVER_SQRT_PI, y.d_); + } +} + +TEST(AgradFwdErfcx, FvarChainRule) { + using stan::math::erfcx; + using stan::math::fvar; + fvar x(1.5, 3.0); + fvar y = erfcx(x); + EXPECT_FLOAT_EQ(3.0 * (2.0 * 1.5 * erfcx(1.5) - stan::math::TWO_OVER_SQRT_PI), + y.d_); +} + +TEST(AgradFwdErfcx, FvarFvarDouble) { + using stan::math::erfcx; + using stan::math::fvar; + + const double x_value = 0.75; + fvar > x; + x.val_.val_ = x_value; + x.val_.d_ = 1.0; + x.d_.val_ = 1.0; + + fvar > y = erfcx(x); + const double v = erfcx(x_value); + const double d1 = 2.0 * x_value * v - stan::math::TWO_OVER_SQRT_PI; + // d2 = 2 * erfcx + 2 * x * d1 + const double d2 = 2.0 * v + 2.0 * x_value * d1; + + EXPECT_FLOAT_EQ(v, y.val_.val_); + EXPECT_FLOAT_EQ(d1, y.val_.d_); + EXPECT_FLOAT_EQ(d1, y.d_.val_); + EXPECT_FLOAT_EQ(d2, y.d_.d_); +} diff --git a/test/unit/math/mix/fun/erfcx_test.cpp b/test/unit/math/mix/fun/erfcx_test.cpp new file mode 100644 index 00000000000..99fb0dcc025 --- /dev/null +++ b/test/unit/math/mix/fun/erfcx_test.cpp @@ -0,0 +1,22 @@ +#include + +// Finite arguments only: erfcx(-inf) is +inf and erfcx grows like +// 2 * exp(x * x) to the left, so the common-argument set (which includes the +// infinities) is not finite-differenceable here. +TEST(mathMixMatFun, erfcx) { + auto f = [](const auto& x1) { return stan::math::erfcx(x1); }; + stan::test::expect_unary_vectorized(f, -2.6, -2, -1, -0.2, 0, 0.5, 1, 1.3, + 2.6, 3.9, 4.1, 12.0); +} + +TEST(mathMixMatFun, erfcxmatvar) { + using stan::math::vec_concat; + using stan::test::expect_ad_vector_matvar; + auto f = [](const auto& x1) { return stan::math::erfcx(x1); }; + std::vector args{-2.6, -2, -1, -0.2, 0, 0.5, 1, 1.3, 2.6, 3.9, 4.1}; + Eigen::VectorXd A(args.size()); + for (int i = 0; i < A.size(); ++i) { + A(i) = args[i]; + } + expect_ad_vector_matvar(f, A); +} diff --git a/test/unit/math/prim/fun/erfcx_test.cpp b/test/unit/math/prim/fun/erfcx_test.cpp new file mode 100644 index 00000000000..a80c5b27577 --- /dev/null +++ b/test/unit/math/prim/fun/erfcx_test.cpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include +#include + +// Reference values computed as expl(x * x) * erfcl(x) in long double +// (64-bit mantissa), which carries ~19 significant digits. +TEST(MathFunctions, erfcx) { + using stan::math::erfcx; + EXPECT_NEAR(108.94090438997797, erfcx(-2.0), 1.1e-12); + EXPECT_NEAR(35.694955906025278, erfcx(-1.7), 3.6e-13); + EXPECT_NEAR(5.0089800807622833, erfcx(-1.0), 5.1e-14); + EXPECT_NEAR(1.952360489182557, erfcx(-0.5), 2.0e-14); + EXPECT_NEAR(1.0, erfcx(0.0), 1e-14); + EXPECT_NEAR(0.97782647768353936, erfcx(0.02), 1e-14); + EXPECT_NEAR(0.6156903441929259, erfcx(0.5), 1e-14); + EXPECT_NEAR(0.427583576155807, erfcx(1.0), 1e-14); + EXPECT_NEAR(0.25539567631050575, erfcx(2.0), 1e-14); + EXPECT_NEAR(0.15126529983237388, erfcx(3.6), 1e-14); + EXPECT_NEAR(0.11070463773306863, erfcx(5.0), 1e-14); +} + +// The whole point of the scaling: both factors of exp(x * x) * erfc(x) are +// out of range here, but the product is an ordinary number. +TEST(MathFunctions, erfcxUpperTail) { + using stan::math::erfcx; + EXPECT_EQ(0.0, std::erfc(30.0)); // the unscaled factor has underflowed + EXPECT_NEAR(0.065925122499980351, erfcx(8.5), 1e-15); + EXPECT_NEAR(0.046854221014893761, erfcx(12.0), 1e-15); + EXPECT_NEAR(0.028174348741051319, erfcx(20.0), 1e-15); + EXPECT_NEAR(0.021683584850562907, erfcx(26.0), 1e-15); + EXPECT_NEAR(0.005641613782989433, erfcx(100.0), 1e-16); + // collapses to the leading term 1 / (x * sqrt(pi)) once every correction + // has fallen below the rounding of the result + EXPECT_NEAR(stan::math::INV_SQRT_PI * 1e-10, erfcx(1e10), 1e-25); +} + +TEST(MathFunctions, erfcxLowerTail) { + using stan::math::erfcx; + EXPECT_NEAR(144009798674.66104, erfcx(-5.0), 1.5e-3); + EXPECT_NEAR(5.3762342836322712e+43, erfcx(-10.0), 5.4e+29); + EXPECT_NEAR(7.6577249314905682e+293, erfcx(-26.0), 7.7e+279); + // 2 * exp(x * x) leaves the binary64 range below about -26.63 + EXPECT_TRUE(std::isinf(erfcx(-27.0))); + EXPECT_GT(erfcx(-27.0), 0.0); +} + +// The internal crossover between exp(x * x) * erfc(x) and the Cody rational +// approximation must not be visible in the output. +TEST(MathFunctions, erfcxBranchContinuity) { + using stan::math::erfcx; + const double cut = 4.0; + for (double delta : {1e-15, 1e-12, 1e-9, 1e-6}) { + const double below = erfcx(cut - cut * delta); + const double above = erfcx(cut + cut * delta); + // the function is smooth and decreasing, so the two sides differ only by + // the slope over 2 * cut * delta + const double slope = 2.0 * cut * erfcx(cut) - stan::math::TWO_OVER_SQRT_PI; + EXPECT_NEAR(below - above, -slope * 2.0 * cut * delta, + 1e-9 * delta + 1e-15); + } + EXPECT_NEAR(erfcx(std::nextafter(cut, 0.0)), erfcx(cut), 1e-16); +} + +TEST(MathFunctions, erfcxEdgeCases) { + using stan::math::erfcx; + const double inf = std::numeric_limits::infinity(); + const double nan = std::numeric_limits::quiet_NaN(); + EXPECT_EQ(1.0, erfcx(0.0)); + EXPECT_EQ(0.0, erfcx(inf)); + EXPECT_TRUE(std::isinf(erfcx(-inf))); + EXPECT_GT(erfcx(-inf), 0.0); + EXPECT_TRUE(std::isnan(erfcx(nan))); +} + +TEST(MathFunctions, erfcxVectorized) { + using stan::math::erfcx; + std::vector xs{-1.0, 0.0, 1.0, 20.0}; + std::vector ys = erfcx(xs); + ASSERT_EQ(4U, ys.size()); + for (size_t i = 0; i < xs.size(); ++i) { + EXPECT_FLOAT_EQ(erfcx(xs[i]), ys[i]); + } + + Eigen::VectorXd v(3); + v << -2.0, 0.5, 12.0; + Eigen::VectorXd w = erfcx(v); + for (int i = 0; i < v.size(); ++i) { + EXPECT_FLOAT_EQ(erfcx(v(i)), w(i)); + } +} + +// erfcx is how the normal tail quantities stay finite; check the two +// identities that motivate the function. +TEST(MathFunctions, erfcxNormalTailIdentities) { + using stan::math::erfcx; + for (double x : {-40.0, -20.0, -7.5, -1.0}) { + const double w = -x * stan::math::INV_SQRT_TWO; + const double log_cdf = stan::math::LOG_HALF + std::log(erfcx(w)) - w * w; + EXPECT_NEAR(log_cdf, stan::math::std_normal_lcdf(x), + 1e-10 * std::fabs(log_cdf)); + // inverse Mills ratio phi(x) / Phi(x) + const double mills = stan::math::SQRT_TWO_OVER_SQRT_PI / erfcx(w); + EXPECT_NEAR( + mills, + std::exp(stan::math::NEG_LOG_SQRT_TWO_PI - 0.5 * x * x - log_cdf), + 1e-9 * mills); + } +} diff --git a/test/unit/math/rev/fun/erfcx_test.cpp b/test/unit/math/rev/fun/erfcx_test.cpp new file mode 100644 index 00000000000..af69ef69a50 --- /dev/null +++ b/test/unit/math/rev/fun/erfcx_test.cpp @@ -0,0 +1,83 @@ +#include +#include +#include +#include +#include + +namespace { +// d/dx erfcx(x) = 2 * x * erfcx(x) - 2 / sqrt(pi) +double expected_deriv(double x) { + return 2.0 * x * stan::math::erfcx(x) - stan::math::TWO_OVER_SQRT_PI; +} +} // namespace + +TEST_F(AgradRev, erfcx_value_and_gradient) { + for (double x_value : {-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0}) { + stan::math::var x = x_value; + stan::math::var y = stan::math::erfcx(x); + y.grad(); + EXPECT_FLOAT_EQ(stan::math::erfcx(x_value), y.val()); + EXPECT_FLOAT_EQ(expected_deriv(x_value), x.adj()); + stan::math::recover_memory(); + } +} + +// Reference derivatives from 2 * x * erfcx(x) - 2 / sqrt(pi) evaluated in +// long double. +TEST_F(AgradRev, erfcx_gradient_reference) { + struct { + double x; + double d; + } cases[] = {{-2.0, -436.89199672700738}, {-1.0, -11.14633932862008}, + {0.0, -1.1283791670955126}, {1.0, -0.27321201478389856}, + {5.0, -0.021332789764826311}, {10.0, -0.0055593122190608565}}; + for (auto c : cases) { + stan::math::var x = c.x; + stan::math::var y = stan::math::erfcx(x); + y.grad(); + EXPECT_NEAR(c.d, x.adj(), 1e-12 * std::fabs(c.d)); + stan::math::recover_memory(); + } +} + +// The gradient is built from the value, so it stays finite exactly where the +// value does -- this is the tail behaviour that motivates the function. +TEST_F(AgradRev, erfcx_upper_tail_gradient) { + stan::math::var x = 30.0; + stan::math::var y = stan::math::erfcx(x); + y.grad(); + EXPECT_TRUE(std::isfinite(y.val())); + EXPECT_GT(y.val(), 0.0); + EXPECT_TRUE(std::isfinite(x.adj())); + EXPECT_NEAR(expected_deriv(30.0), x.adj(), 1e-15); +} + +TEST_F(AgradRev, erfcx_matrix) { + Eigen::Matrix x(3); + x << -1.0, 0.5, 12.0; + stan::math::var y = stan::math::sum(stan::math::erfcx(x)); + y.grad(); + for (int i = 0; i < x.size(); ++i) { + EXPECT_FLOAT_EQ(expected_deriv(x(i).val()), x(i).adj()); + } +} + +TEST_F(AgradRev, erfcx_var_matrix) { + Eigen::VectorXd xd(3); + xd << -1.0, 0.5, 12.0; + stan::math::var_value x(xd); + auto fx = stan::math::erfcx(x); + stan::math::var y = stan::math::sum(fx); + + double expected_sum = 0.0; + for (int i = 0; i < xd.size(); ++i) { + EXPECT_FLOAT_EQ(stan::math::erfcx(xd(i)), fx.val()(i)); + expected_sum += stan::math::erfcx(xd(i)); + } + EXPECT_FLOAT_EQ(expected_sum, y.val()); + + y.grad(); + for (int i = 0; i < xd.size(); ++i) { + EXPECT_FLOAT_EQ(expected_deriv(xd(i)), x.adj()(i)); + } +} From e0707b150d257e7cd8300072eca9cc9eeb2f2a28 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 16 Sep 2026 09:45:59 +0300 Subject: [PATCH 02/24] Flatten erfcx Horner evaluation to satisfy cpplint line length --- .../opencl/kernels/device_functions/erfcx.hpp | 34 +++++++------------ stan/math/prim/fun/erfcx.hpp | 33 +++++++----------- 2 files changed, 24 insertions(+), 43 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 032d8c70eaf..3909363e4f0 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -32,28 +32,18 @@ static constexpr const char* erfcx_device_function double erfcx(double x) { if (x >= 4.0) { double u = 1.0 / (x * x); - double p - = 0.000658749161529837803157 - + u - * (0.0160837851487422766278 - + u - * (0.125781726111229246204 - + u - * (0.360344899949804439429 - + u - * (0.305326634961232344035 - + u * 0.0163153871373020978498)))); - double q - = -0.00233520497626869185443 - + u - * (-0.0605183413124413191178 - + u - * (-0.527905102951428412248 - + u - * (-1.87295284992346047209 - + u - * (-2.56852019228982242072 - + u * -1.0)))); + double p = 0.0163153871373020978498; + p = 0.305326634961232344035 + u * p; + p = 0.360344899949804439429 + u * p; + p = 0.125781726111229246204 + u * p; + p = 0.0160837851487422766278 + u * p; + p = 0.000658749161529837803157 + u * p; + double q = -1.0; + q = -2.56852019228982242072 + u * q; + q = -1.87295284992346047209 + u * q; + q = -0.527905102951428412248 + u * q; + q = -0.0605183413124413191178 + u * q; + q = -0.00233520497626869185443 + u * q; return (M_2_SQRTPI * 0.5 + (p / q) * u) / x; } if (x < -27.0) { diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 006562f1f8b..06060e9fafe 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -76,27 +76,18 @@ inline double erfcx(T&& xx) { constexpr double cody_min = 4.0; if (x >= cody_min) { const double u = 1.0 / (x * x); - const double p - = 0.000658749161529837803157 - + u - * (0.0160837851487422766278 - + u - * (0.125781726111229246204 - + u - * (0.360344899949804439429 - + u - * (0.305326634961232344035 - + u * 0.0163153871373020978498)))); - const double q = -0.00233520497626869185443 - + u - * (-0.0605183413124413191178 - + u - * (-0.527905102951428412248 - + u - * (-1.87295284992346047209 - + u - * (-2.56852019228982242072 - + u * -1.0)))); + double p = 0.0163153871373020978498; + p = 0.305326634961232344035 + u * p; + p = 0.360344899949804439429 + u * p; + p = 0.125781726111229246204 + u * p; + p = 0.0160837851487422766278 + u * p; + p = 0.000658749161529837803157 + u * p; + double q = -1.0; + q = -2.56852019228982242072 + u * q; + q = -1.87295284992346047209 + u * q; + q = -0.527905102951428412248 + u * q; + q = -0.0605183413124413191178 + u * q; + q = -0.00233520497626869185443 + u * q; return (INV_SQRT_PI + (p / q) * u) / x; } // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves From 266fb5c7ecb2e43b75d0eb0252b0f78557d434e1 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 00:10:23 +0300 Subject: [PATCH 03/24] Use fma instead of Dekker split in OpenCL erfcx device function --- .../opencl/kernels/device_functions/erfcx.hpp | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 3909363e4f0..ff54d2e0cfd 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -22,9 +22,20 @@ static constexpr const char* erfcx_device_function * * Mirrors stan/math/prim/fun/erfcx.hpp: the Cody (1969) rational * approximation above 4, and below it exp(x * x) * erfc(x) with - * x * x Dekker-split so that exp() receives an exact argument. - * Keep the two in step; the split is what holds the lower branch - * to a few ulp instead of x * x * eps. + * the rounding of x * x corrected so that exp() effectively + * receives an exact argument. Without that correction the error + * of x * x is amplified by exp into roughly x * x * eps, which is + * 512 ulp at x = -26. + * + * The correction is written with fma() rather than the Dekker + * split used on the host. The split relies on t - (t - x) not + * being simplified to x, which is true in floating point but not + * over the reals, and the OpenCL compiler simplifies it: measured + * on a Tesla V100 (OpenCL 3.0 CUDA), the split form returned + * x_lo == 0 at every one of 4096 test points and scored 512.66 + * ulp, bit-identical to the uncorrected formula. fma() is a + * single instruction, so there is nothing to reassociate; it + * measures 3.40 ulp on the same device. * * @param x argument * @return scaled complementary error function of the argument @@ -49,10 +60,8 @@ static constexpr const char* erfcx_device_function if (x < -27.0) { return INFINITY; } - double t = 134217729.0 * x; - double x_hi = t - (t - x); - double x_lo = x - x_hi; - return exp(x_hi * x_hi) * exp(x_lo * (x + x_hi)) * erfc(x); + double h = x * x; + return exp(h) * (1.0 + fma(x, x, -h)) * erfc(x); } // \cond ) "\n#endif\n"; // NOLINT From 8c16f29fd416fba2ef71bd0b008d51cfc3103632 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 00:15:15 +0300 Subject: [PATCH 04/24] Use fma instead of Dekker split in erfcx lower branch --- stan/math/prim/fun/erfcx.hpp | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 06060e9fafe..1002e663f9c 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -33,11 +33,11 @@ namespace math { * directly with no `exp` and no `erfc` call. Measured at 2.0 to 2.5 ulp * over `[4, 50]` and 5.4 ulp out to 200, and about 2.4 times faster than * forming the product. - * - `x < 4`: `exp(x * x) * erfc(x)`, with `x * x` split so that the - * argument handed to `exp` is exact. Without the split the rounding of - * `x * x` is amplified by `exp` into a relative error of about - * `x * x * eps`, which reaches 512 ulp at `x = -26`. Splitting holds this - * branch to under 5 ulp. + * - `x < 4`: `exp(x * x) * erfc(x)`, with the rounding error of `x * x` + * recovered by `fma` and folded back in, so `exp` effectively receives an + * exact argument. Uncorrected, the rounding of `x * x` is amplified by + * `exp` into a relative error of about `x * x * eps`, which reaches + * 512 ulp at `x = -26`. The correction holds this branch to under 6 ulp. * * W. J. Cody, Math. Comp. 23(107):631-638 (1969). * https://doi.org/10.1090/S0025-5718-1969-0247736-4 @@ -91,20 +91,18 @@ inline double erfcx(T&& xx) { return (INV_SQRT_PI + (p / q) * u) / x; } // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves - // the binary64 range below -26.63. Returning here also keeps the split - // below from forming inf - inf on -infinity. + // the binary64 range below -26.63. Returning here also keeps the branch + // below from forming inf * 0 on -infinity. constexpr double overflow_max = -27.0; if (x < overflow_max) { return INFTY; } - // Dekker split: x_hi holds the leading 26 bits of x, so x_hi * x_hi is - // exact and exp() is handed an argument carrying no rounding error of - // its own. x_lo * (x + x_hi) is the remainder x * x - x_hi * x_hi. - constexpr double split = 134217729.0; // 2^27 + 1 - const double t = split * x; - const double x_hi = t - (t - x); - const double x_lo = x - x_hi; - return std::exp(x_hi * x_hi) * std::exp(x_lo * (x + x_hi)) * std::erfc(x); + // fma(x, x, -h) is the exact rounding error of h = x * x, so the factor + // (1 + that) restores what exp(h) would otherwise lose. Same correction + // as the OpenCL device function, which cannot use a Dekker split because + // the compiler simplifies t - (t - x) to x. + const double h = x * x; + return std::exp(h) * (1.0 + std::fma(x, x, -h)) * std::erfc(x); } /** From ea151bea9d58af189fa7b7716c739a66d17261a4 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 15:44:26 +0300 Subject: [PATCH 05/24] Use Cody second-interval rational for erfcx on 0.46875 <= x < 4 --- stan/math/prim/fun/erfcx.hpp | 117 +++++++++++++++++++++++++---------- 1 file changed, 84 insertions(+), 33 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 1002e663f9c..528d41f94ad 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -9,6 +9,69 @@ namespace stan { namespace math { +namespace internal { + +/** + * Cody (1969) third-interval rational, valid for `x >= 4`. + * + * Gives `erfcx` directly. `x * x` is infinite for `x` large enough, which + * correctly collapses this to the leading term `INV_SQRT_PI / x` and, at + * infinity, to zero. + * + * @param x argument, `x >= 4` + * @return scaled complementary error function + */ +inline double erfcx_cody_tail(double x) { + const double u = 1.0 / (x * x); + double p = 0.0163153871373020978498; + p = 0.305326634961232344035 + u * p; + p = 0.360344899949804439429 + u * p; + p = 0.125781726111229246204 + u * p; + p = 0.0160837851487422766278 + u * p; + p = 0.000658749161529837803157 + u * p; + double q = -1.0; + q = -2.56852019228982242072 + u * q; + q = -1.87295284992346047209 + u * q; + q = -0.527905102951428412248 + u * q; + q = -0.0605183413124413191178 + u * q; + q = -0.00233520497626869185443 + u * q; + return (INV_SQRT_PI + (p / q) * u) / x; +} + +/** + * Cody (1969) second-interval rational, valid for `0.46875 <= x <= 4`. + * + * Cody writes this interval as `erfc(x) = exp(-x*x) * R(x)`, so `erfcx` + * is `R(x)` with the exponential cancelled analytically: no `exp` and no + * `erfc` call. Measured at 3.7 to 7.0 ulp, of which only 0.28 ulp is + * approximation error; the rest is rounding in the two degree-8 Horner + * chains and the final division. + * + * @param y argument, `0.46875 <= y <= 4` + * @return scaled complementary error function + */ +inline double erfcx_cody_middle(double y) { + double p = 2.15311535474403846e-8 * y; + p = (p + 5.64188496988670089e-1) * y; + p = (p + 8.88314979438837594) * y; + p = (p + 66.1191906371416295) * y; + p = (p + 298.635138197400131) * y; + p = (p + 881.952221241769090) * y; + p = (p + 1712.04761263407058) * y; + p = (p + 2051.07837782607147) * y; + double q = y; + q = (q + 15.7449261107098347) * y; + q = (q + 117.693950891312499) * y; + q = (q + 537.181101862009858) * y; + q = (q + 1621.38957456669019) * y; + q = (q + 3290.79923573345963) * y; + q = (q + 4362.61909014324716) * y; + q = (q + 3439.36767414372164) * y; + return (p + 1230.33935479799725) / (q + 1230.33935480374942); +} + +} // namespace internal + /** * Return the scaled complementary error function of the argument. * @@ -25,19 +88,19 @@ namespace math { * Mills ratio is `SQRT_TWO_OVER_SQRT_PI / erfcx(-x * INV_SQRT_TWO)`, both * without cancellation. * - * Two branches, with the crossover and the accuracy below established by - * measurement against a `long double` reference over 40000 points per - * interval: + * Three branches. `erfc` dominates the cost of forming + * `exp(x*x) * erfc(x)`, so each interval that can be served by a direct + * approximation of `erfcx` instead is much faster. * - * - `x >= 4`: the Cody (1969) rational approximation, which gives `erfcx` - * directly with no `exp` and no `erfc` call. Measured at 2.0 to 2.5 ulp - * over `[4, 50]` and 5.4 ulp out to 200, and about 2.4 times faster than - * forming the product. - * - `x < 4`: `exp(x * x) * erfc(x)`, with the rounding error of `x * x` - * recovered by `fma` and folded back in, so `exp` effectively receives an - * exact argument. Uncorrected, the rounding of `x * x` is amplified by - * `exp` into a relative error of about `x * x * eps`, which reaches - * 512 ulp at `x = -26`. The correction holds this branch to under 6 ulp. + * - `x >= 4`: Cody (1969) third-interval rational. + * - `0.46875 <= x < 4`: Cody (1969) second-interval rational. Cody writes + * this interval as `erfc(x) = exp(-x*x) * R(x)`, so `erfcx` is `R(x)` + * and the exponential is cancelled analytically rather than computed and + * divided out. No `exp` and no `erfc` call. + * - `x < 0.46875`: `exp(x * x) * erfc(x)`, with the rounding error of + * `x * x` recovered by `fma` and folded back in, because `exp` + * amplifies it into roughly `x * x * eps` -- 512 ulp at `x = -26` if + * left alone. * * W. J. Cody, Math. Comp. 23(107):631-638 (1969). * https://doi.org/10.1090/S0025-5718-1969-0247736-4 @@ -49,7 +112,7 @@ namespace math { * relative down to 4 and degrades past about 3.5, so 4 sits just inside its * range. * - * Measured worst case over the whole range is 5.4 ulp. + * Measured worst case over the whole range is 7.0 ulp. * * The derivative is * \f$\frac{d}{dx}\mbox{erfcx}(x) = 2x\,\mbox{erfcx}(x) - @@ -70,25 +133,13 @@ namespace math { template * = nullptr> inline double erfcx(T&& xx) { const double x = static_cast(xx); - // Cody (1969). x * x is infinite for x large enough, which correctly - // collapses this to the leading term INV_SQRT_PI / x and, at infinity, - // to zero. constexpr double cody_min = 4.0; + constexpr double middle_min = 0.46875; if (x >= cody_min) { - const double u = 1.0 / (x * x); - double p = 0.0163153871373020978498; - p = 0.305326634961232344035 + u * p; - p = 0.360344899949804439429 + u * p; - p = 0.125781726111229246204 + u * p; - p = 0.0160837851487422766278 + u * p; - p = 0.000658749161529837803157 + u * p; - double q = -1.0; - q = -2.56852019228982242072 + u * q; - q = -1.87295284992346047209 + u * q; - q = -0.527905102951428412248 + u * q; - q = -0.0605183413124413191178 + u * q; - q = -0.00233520497626869185443 + u * q; - return (INV_SQRT_PI + (p / q) * u) / x; + return internal::erfcx_cody_tail(x); + } + if (x >= middle_min) { + return internal::erfcx_cody_middle(x); } // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves // the binary64 range below -26.63. Returning here also keeps the branch @@ -98,9 +149,9 @@ inline double erfcx(T&& xx) { return INFTY; } // fma(x, x, -h) is the exact rounding error of h = x * x, so the factor - // (1 + that) restores what exp(h) would otherwise lose. Same correction - // as the OpenCL device function, which cannot use a Dekker split because - // the compiler simplifies t - (t - x) to x. + // (1 + that) restores what exp(h) would otherwise lose. The OpenCL + // device function uses the same form; it cannot use a Dekker split + // because the compiler simplifies t - (t - x) to x. const double h = x * x; return std::exp(h) * (1.0 + std::fma(x, x, -h)) * std::erfc(x); } From f55de06dd2bdb0b0404ae7a3b82b44d6cc56517f Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 15:44:35 +0300 Subject: [PATCH 06/24] Replace erfcx erfc call below 0.46875 with a Chebyshev expansion and reflection --- stan/math/prim/fun/erfcx.hpp | 80 +++++++++++++++++++++++++++++------- 1 file changed, 66 insertions(+), 14 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 528d41f94ad..e689f611d2a 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -43,9 +43,7 @@ inline double erfcx_cody_tail(double x) { * * Cody writes this interval as `erfc(x) = exp(-x*x) * R(x)`, so `erfcx` * is `R(x)` with the exponential cancelled analytically: no `exp` and no - * `erfc` call. Measured at 3.7 to 7.0 ulp, of which only 0.28 ulp is - * approximation error; the rest is rounding in the two degree-8 Horner - * chains and the final division. + * `erfc` call. Measured at 3.7 to 7.0 ulp. * * @param y argument, `0.46875 <= y <= 4` * @return scaled complementary error function @@ -70,6 +68,44 @@ inline double erfcx_cody_middle(double y) { return (p + 1230.33935479799725) / (q + 1230.33935480374942); } +/** + * Degree-18 polynomial for `|x| < 0.46875`. + * + * A Chebyshev-economized expansion of `erfcx` about zero, so it needs no + * library call and no branch on the sign of `x`. Its low-order coefficients + * reproduce the Maclaurin series of `erfcx` exactly (`1`, `-2/sqrt(pi)`, + * `1`, ...), which is a useful check that the fit is right. Measured at + * 2.2 ulp over the whole interval. + * + * The plain Maclaurin series is only usable to about `|x| = 0.125`; four + * further economized terms extend it to 0.46875 at no measurable cost, + * which is what removes the last `exp` call from the positive axis. + * + * @param x argument, `|x| < 0.46875` + * @return scaled complementary error function + */ +inline double erfcx_small(double x) { + double p = 3.05977060678449757e-06; + p = -9.35890030086883823e-06 + x * p; + p = 2.46655529768908249e-05 + x * p; + p = -7.08163358203131886e-05 + x * p; + p = 1.98445679338826757e-04 + x * p; + p = -5.34506929034156810e-04 + x * p; + p = 1.38888415444527033e-03 + x * p; + p = -3.47359067853470795e-03 + x * p; + p = 8.33333374332981443e-03 + x * p; + p = -1.91048337772546720e-02 + x * p; + p = 4.16666666458337179e-02 + x * p; + p = -8.59717459974174147e-02 + x * p; + p = 1.66666666667239644e-01 + x * p; + p = -3.00901111227312890e-01 + x * p; + p = 4.99999999999992839e-01 + x * p; + p = -7.52252778063651983e-01 + x * p; + p = 1.0 + x * p; + p = -1.12837916709551256 + x * p; + return 1.0 + x * p; +} + } // namespace internal /** @@ -88,19 +124,25 @@ inline double erfcx_cody_middle(double y) { * Mills ratio is `SQRT_TWO_OVER_SQRT_PI / erfcx(-x * INV_SQRT_TWO)`, both * without cancellation. * - * Three branches. `erfc` dominates the cost of forming - * `exp(x*x) * erfc(x)`, so each interval that can be served by a direct - * approximation of `erfcx` instead is much faster. + * Four branches. The whole positive axis is covered without calling any + * library function, which is what makes this fast: a direct approximation + * of `erfcx` avoids `erfc`, and `erfc` dominates the cost of forming + * `exp(x*x) * erfc(x)`. * * - `x >= 4`: Cody (1969) third-interval rational. * - `0.46875 <= x < 4`: Cody (1969) second-interval rational. Cody writes * this interval as `erfc(x) = exp(-x*x) * R(x)`, so `erfcx` is `R(x)` * and the exponential is cancelled analytically rather than computed and - * divided out. No `exp` and no `erfc` call. - * - `x < 0.46875`: `exp(x * x) * erfc(x)`, with the rounding error of - * `x * x` recovered by `fma` and folded back in, because `exp` - * amplifies it into roughly `x * x * eps` -- 512 ulp at `x = -26` if - * left alone. + * divided out. + * - `|x| < 0.46875`: a degree-18 Chebyshev-economized expansion about + * zero, covering both signs with no branch and no library call. + * - `x <= -0.46875`: the reflection `2*exp(x*x) - erfcx(-x)`, which is the + * one place a library call is unavoidable, since `erfcx` grows like + * `2*exp(x*x)` as `x -> -infinity`. The rounding error of `x * x` is + * recovered with `fma` and folded back in, because `exp` amplifies it + * into roughly `x * x * eps` -- 512 ulp at `x = -26` if left alone. + * Below `x = -6.1`, `erfcx(-x)` is under `eps/2` of `2*exp(x*x)`, so + * the subtraction is skipped. * * W. J. Cody, Math. Comp. 23(107):631-638 (1969). * https://doi.org/10.1090/S0025-5718-1969-0247736-4 @@ -141,9 +183,12 @@ inline double erfcx(T&& xx) { if (x >= middle_min) { return internal::erfcx_cody_middle(x); } + if (x > -middle_min) { + return internal::erfcx_small(x); + } // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves - // the binary64 range below -26.63. Returning here also keeps the branch - // below from forming inf * 0 on -infinity. + // the binary64 range below -26.63. Returning here also keeps the + // reflection below from forming inf * 0 on -infinity. constexpr double overflow_max = -27.0; if (x < overflow_max) { return INFTY; @@ -153,7 +198,14 @@ inline double erfcx(T&& xx) { // device function uses the same form; it cannot use a Dekker split // because the compiler simplifies t - (t - x) to x. const double h = x * x; - return std::exp(h) * (1.0 + std::fma(x, x, -h)) * std::erfc(x); + const double two_exp_x2 = 2.0 * std::exp(h) * (1.0 + std::fma(x, x, -h)); + constexpr double reflect_min = -6.1; + if (x < reflect_min) { + return two_exp_x2; + } + const double y = -x; + return two_exp_x2 - (y >= cody_min ? internal::erfcx_cody_tail(y) + : internal::erfcx_cody_middle(y)); } /** From 7dfed733f506778d3c8e82708cc1bef99b0198d5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 15:44:51 +0300 Subject: [PATCH 07/24] Mirror the new erfcx branches in the OpenCL device function and tests --- .../opencl/kernels/device_functions/erfcx.hpp | 141 ++++++++++++++---- test/unit/math/prim/fun/erfcx_test.cpp | 47 ++++-- 2 files changed, 146 insertions(+), 42 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index ff54d2e0cfd..86481a0279a 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -15,53 +15,134 @@ static constexpr const char* erfcx_device_function "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX\n" "#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX\n" STRINGIFY( // \endcond + /** \ingroup opencl_kernels + * + * Cody (1969) third-interval rational, for x >= 4. + * + * @param x argument + * @return scaled complementary error function + */ + double erfcx_cody_tail(double x) { + double u = 1.0 / (x * x); + double p = 0.0163153871373020978498; + p = 0.305326634961232344035 + u * p; + p = 0.360344899949804439429 + u * p; + p = 0.125781726111229246204 + u * p; + p = 0.0160837851487422766278 + u * p; + p = 0.000658749161529837803157 + u * p; + double q = -1.0; + q = -2.56852019228982242072 + u * q; + q = -1.87295284992346047209 + u * q; + q = -0.527905102951428412248 + u * q; + q = -0.0605183413124413191178 + u * q; + q = -0.00233520497626869185443 + u * q; + return (M_2_SQRTPI * 0.5 + (p / q) * u) / x; + } + + /** \ingroup opencl_kernels + * + * Cody (1969) second-interval rational, for 0.46875 <= x <= 4. + * Yields erfcx directly: the exponential is cancelled + * analytically, so there is no exp and no erfc call. + * + * @param y argument + * @return scaled complementary error function + */ + double erfcx_cody_middle(double y) { + double p = 2.15311535474403846e-8 * y; + p = (p + 5.64188496988670089e-1) * y; + p = (p + 8.88314979438837594) * y; + p = (p + 66.1191906371416295) * y; + p = (p + 298.635138197400131) * y; + p = (p + 881.952221241769090) * y; + p = (p + 1712.04761263407058) * y; + p = (p + 2051.07837782607147) * y; + double q = y; + q = (q + 15.7449261107098347) * y; + q = (q + 117.693950891312499) * y; + q = (q + 537.181101862009858) * y; + q = (q + 1621.38957456669019) * y; + q = (q + 3290.79923573345963) * y; + q = (q + 4362.61909014324716) * y; + q = (q + 3439.36767414372164) * y; + return (p + 1230.33935479799725) / (q + 1230.33935480374942); + } + + /** \ingroup opencl_kernels + * + * Degree-18 Chebyshev-economized expansion of erfcx about zero, + * for |x| < 0.46875. Covers both signs with no branch and no + * library call. + * + * @param x argument + * @return scaled complementary error function + */ + double erfcx_small(double x) { + double p = 3.05977060678449757e-06; + p = -9.35890030086883823e-06 + x * p; + p = 2.46655529768908249e-05 + x * p; + p = -7.08163358203131886e-05 + x * p; + p = 1.98445679338826757e-04 + x * p; + p = -5.34506929034156810e-04 + x * p; + p = 1.38888415444527033e-03 + x * p; + p = -3.47359067853470795e-03 + x * p; + p = 8.33333374332981443e-03 + x * p; + p = -1.91048337772546720e-02 + x * p; + p = 4.16666666458337179e-02 + x * p; + p = -8.59717459974174147e-02 + x * p; + p = 1.66666666667239644e-01 + x * p; + p = -3.00901111227312890e-01 + x * p; + p = 4.99999999999992839e-01 + x * p; + p = -7.52252778063651983e-01 + x * p; + p = 1.0 + x * p; + p = -1.12837916709551256 + x * p; + return 1.0 + x * p; + } + /** \ingroup opencl_kernels * * Return the scaled complementary error function * exp(x * x) * erfc(x) of the kernel generator expression. * - * Mirrors stan/math/prim/fun/erfcx.hpp: the Cody (1969) rational - * approximation above 4, and below it exp(x * x) * erfc(x) with - * the rounding of x * x corrected so that exp() effectively - * receives an exact argument. Without that correction the error - * of x * x is amplified by exp into roughly x * x * eps, which is - * 512 ulp at x = -26. + * Mirrors stan/math/prim/fun/erfcx.hpp branch for branch and + * formula for formula; the two must be kept in step. The whole + * positive axis is covered without a library call, which is + * what makes this fast. On the negative side exp is + * unavoidable, since erfcx grows like 2*exp(x*x); its argument + * is corrected with fma, because exp amplifies the rounding of + * x * x into roughly x * x * eps, which is 512 ulp at x = -26. * - * The correction is written with fma() rather than the Dekker - * split used on the host. The split relies on t - (t - x) not - * being simplified to x, which is true in floating point but not - * over the reals, and the OpenCL compiler simplifies it: measured - * on a Tesla V100 (OpenCL 3.0 CUDA), the split form returned - * x_lo == 0 at every one of 4096 test points and scored 512.66 - * ulp, bit-identical to the uncorrected formula. fma() is a - * single instruction, so there is nothing to reassociate; it - * measures 3.40 ulp on the same device. + * The correction must be written with fma and not with a Dekker + * split. The split relies on t - (t - x) not being simplified + * to x, which holds in floating point but not over the reals, + * and the OpenCL compiler does simplify it: the split form + * silently degrades to the uncorrected product on device while + * still being correct on the host. Do not reintroduce it. * * @param x argument * @return scaled complementary error function of the argument */ double erfcx(double x) { if (x >= 4.0) { - double u = 1.0 / (x * x); - double p = 0.0163153871373020978498; - p = 0.305326634961232344035 + u * p; - p = 0.360344899949804439429 + u * p; - p = 0.125781726111229246204 + u * p; - p = 0.0160837851487422766278 + u * p; - p = 0.000658749161529837803157 + u * p; - double q = -1.0; - q = -2.56852019228982242072 + u * q; - q = -1.87295284992346047209 + u * q; - q = -0.527905102951428412248 + u * q; - q = -0.0605183413124413191178 + u * q; - q = -0.00233520497626869185443 + u * q; - return (M_2_SQRTPI * 0.5 + (p / q) * u) / x; + return erfcx_cody_tail(x); + } + if (x >= 0.46875) { + return erfcx_cody_middle(x); + } + if (x > -0.46875) { + return erfcx_small(x); } if (x < -27.0) { return INFINITY; } double h = x * x; - return exp(h) * (1.0 + fma(x, x, -h)) * erfc(x); + double two_exp_x2 = 2.0 * exp(h) * (1.0 + fma(x, x, -h)); + if (x < -6.1) { + return two_exp_x2; + } + double y = -x; + return two_exp_x2 + - (y >= 4.0 ? erfcx_cody_tail(y) : erfcx_cody_middle(y)); } // \cond ) "\n#endif\n"; // NOLINT diff --git a/test/unit/math/prim/fun/erfcx_test.cpp b/test/unit/math/prim/fun/erfcx_test.cpp index a80c5b27577..9f9d272889e 100644 --- a/test/unit/math/prim/fun/erfcx_test.cpp +++ b/test/unit/math/prim/fun/erfcx_test.cpp @@ -46,21 +46,44 @@ TEST(MathFunctions, erfcxLowerTail) { EXPECT_GT(erfcx(-27.0), 0.0); } -// The internal crossover between exp(x * x) * erfc(x) and the Cody rational -// approximation must not be visible in the output. +// None of the internal crossovers may be visible in the output. There are +// three: the small expansion to the Cody second interval at 0.46875, the +// second to the third interval at 4, and the small expansion to the +// negative-side reflection at -0.46875. TEST(MathFunctions, erfcxBranchContinuity) { using stan::math::erfcx; - const double cut = 4.0; - for (double delta : {1e-15, 1e-12, 1e-9, 1e-6}) { - const double below = erfcx(cut - cut * delta); - const double above = erfcx(cut + cut * delta); - // the function is smooth and decreasing, so the two sides differ only by - // the slope over 2 * cut * delta - const double slope = 2.0 * cut * erfcx(cut) - stan::math::TWO_OVER_SQRT_PI; - EXPECT_NEAR(below - above, -slope * 2.0 * cut * delta, - 1e-9 * delta + 1e-15); + for (double cut : {4.0, 0.46875, -0.46875}) { + for (double delta : {1e-15, 1e-12, 1e-9, 1e-6}) { + const double step = std::fabs(cut) * delta; + const double below = erfcx(cut - step); + const double above = erfcx(cut + step); + // the function is smooth, so the two sides differ only by the slope + // over 2 * step + const double slope + = 2.0 * cut * erfcx(cut) - stan::math::TWO_OVER_SQRT_PI; + EXPECT_NEAR(below - above, -slope * 2.0 * step, + 1e-9 * delta + 1e-14 * std::fabs(erfcx(cut))); + } + const double inward = (cut > 0.0) ? 0.0 : -1.0; + EXPECT_NEAR(erfcx(std::nextafter(cut, inward)), erfcx(cut), + 1e-15 * std::fabs(erfcx(cut))); } - EXPECT_NEAR(erfcx(std::nextafter(cut, 0.0)), erfcx(cut), 1e-16); +} + +// Below -6.1 the reflection skips subtracting erfcx(-x). That is only +// legitimate if the skipped term really is below half an ulp of +// 2*exp(x*x), which is what this asserts directly. A continuity sweep is +// not usable here: erfcx has slope about -3.5e17 at -6.1, so the +// function's own variation swamps any fixed relative tolerance. +TEST(MathFunctions, erfcxReflectionSkippedTerm) { + const double cut = -6.1; + const double leading = 2.0 * std::exp(cut * cut); + const double skipped = stan::math::erfcx(-cut); + EXPECT_LT(skipped / leading, 0.5 * 1.1102230246251565e-16); + // adjacent representable arguments still agree to the steepness limit + const double below = stan::math::erfcx(std::nextafter(cut, -1.0)); + const double above = stan::math::erfcx(std::nextafter(cut, 0.0)); + EXPECT_NEAR(below / above, 1.0, 1e-13); } TEST(MathFunctions, erfcxEdgeCases) { From 37d8262b3cce3d4ae03ae5b32a21b8e0dce58d3d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 17 Sep 2026 09:05:27 -0400 Subject: [PATCH 08/24] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/erfcx.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index e689f611d2a..adb04216726 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -204,8 +204,9 @@ inline double erfcx(T&& xx) { return two_exp_x2; } const double y = -x; - return two_exp_x2 - (y >= cody_min ? internal::erfcx_cody_tail(y) - : internal::erfcx_cody_middle(y)); + return two_exp_x2 + - (y >= cody_min ? internal::erfcx_cody_tail(y) + : internal::erfcx_cody_middle(y)); } /** From 6ff45419639fcee90e1f069804f005f86531f7f8 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 23:51:01 +0300 Subject: [PATCH 09/24] Add reverse-mode OpenCL support for erfcx --- stan/math/opencl/rev.hpp | 1 + stan/math/opencl/rev/erfcx.hpp | 37 ++++++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+) create mode 100644 stan/math/opencl/rev/erfcx.hpp diff --git a/stan/math/opencl/rev.hpp b/stan/math/opencl/rev.hpp index 37c5b716230..4d8d8e53124 100644 --- a/stan/math/opencl/rev.hpp +++ b/stan/math/opencl/rev.hpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/opencl/rev/erfcx.hpp b/stan/math/opencl/rev/erfcx.hpp new file mode 100644 index 00000000000..bb77a027b97 --- /dev/null +++ b/stan/math/opencl/rev/erfcx.hpp @@ -0,0 +1,37 @@ +#ifndef STAN_MATH_OPENCL_REV_ERFCX_HPP +#define STAN_MATH_OPENCL_REV_ERFCX_HPP +#ifdef STAN_OPENCL + +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the elementwise `erfcx()` of a var_value>. + * + * The derivative `2 * x * erfcx(x) - 2 / sqrt(pi)` reuses the function + * value, so no second `exp` or `erfc` evaluation is needed and the + * derivative inherits the value's accuracy in the tails. + * + * @param A argument + * @return Elementwise `erfcx()` of the input. + */ +template * = nullptr> +inline var_value> erfcx(const var_value& A) { + return make_callback_var( + erfcx(A.val()), [A](vari_value>& res) mutable { + A.adj() += elt_multiply( + res.adj(), elt_multiply(2.0, elt_multiply(A.val(), res.val())) + - TWO_OVER_SQRT_PI); + }); +} + +} // namespace math +} // namespace stan + +#endif +#endif From 66304496ac38ff60796c9480bd2a0c1e3e7bc635 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 17 Sep 2026 23:51:07 +0300 Subject: [PATCH 10/24] Add OpenCL tests for erfcx --- .../kernel_generator/elt_function_cl_test.cpp | 1 + test/unit/math/opencl/rev/erfcx_test.cpp | 29 +++++++++++++++++++ 2 files changed, 30 insertions(+) create mode 100644 test/unit/math/opencl/rev/erfcx_test.cpp diff --git a/test/unit/math/opencl/kernel_generator/elt_function_cl_test.cpp b/test/unit/math/opencl/kernel_generator/elt_function_cl_test.cpp index b97b945172f..f23bce2e8d7 100644 --- a/test/unit/math/opencl/kernel_generator/elt_function_cl_test.cpp +++ b/test/unit/math/opencl/kernel_generator/elt_function_cl_test.cpp @@ -93,6 +93,7 @@ TEST_UNARY_FUNCTION(tgamma) TEST_UNARY_FUNCTION(lgamma) TEST_UNARY_FUNCTION(erf) TEST_UNARY_FUNCTION(erfc) +TEST_UNARY_FUNCTION(erfcx) TEST_UNARY_FUNCTION(floor) TEST_UNARY_FUNCTION(round) diff --git a/test/unit/math/opencl/rev/erfcx_test.cpp b/test/unit/math/opencl/rev/erfcx_test.cpp new file mode 100644 index 00000000000..ddbba02a8c4 --- /dev/null +++ b/test/unit/math/opencl/rev/erfcx_test.cpp @@ -0,0 +1,29 @@ +#ifdef STAN_OPENCL +#include +#include +#include +#include + +auto erfcx_functor = [](const auto& a) { return stan::math::erfcx(a); }; + +TEST(OpenCLerfcx, prim_rev_values_small) { + Eigen::VectorXd a(7); + a << -2.6, -2, -1, -0.2, 1, 1.3, 2.6; + stan::math::test::compare_cpu_opencl_prim_rev(erfcx_functor, a); +} + +TEST(OpenCLerfcx, prim_rev_size_0) { + int N = 0; + + Eigen::MatrixXd a(N, N); + stan::math::test::compare_cpu_opencl_prim_rev(erfcx_functor, a); +} + +TEST(OpenCLerfcx, prim_rev_values_large) { + int N = 71; + + Eigen::MatrixXd a = Eigen::MatrixXd::Random(N, N); + stan::math::test::compare_cpu_opencl_prim_rev(erfcx_functor, a); +} + +#endif From b9c353fdb588c12cd5e186fc7e0ca31171de345e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 12:59:20 +0300 Subject: [PATCH 11/24] Add fixed-reference tests for the erfcx derivative in the tail --- .../math/mix/fun/erfcx_derivative_test.cpp | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 test/unit/math/mix/fun/erfcx_derivative_test.cpp diff --git a/test/unit/math/mix/fun/erfcx_derivative_test.cpp b/test/unit/math/mix/fun/erfcx_derivative_test.cpp new file mode 100644 index 00000000000..31bdfdc4cb3 --- /dev/null +++ b/test/unit/math/mix/fun/erfcx_derivative_test.cpp @@ -0,0 +1,102 @@ +#include +#include +#include + +/** + * Accuracy of the erfcx derivative in the tail. + * + * `d/dx erfcx(x) = 2 * x * erfcx(x) - 2 / sqrt(pi)`. For large `x` the two + * terms both approach `2 / sqrt(pi)` while their difference decays like + * `1 / (sqrt(pi) * x^2)`, so the subtraction loses the result. Measured + * against a 50-digit reference, the difference form gives 6.1 ulp at x = 4 + * and 2.55e+11 ulp at x = 1e6. + * + * For `x >= 4` the same tail rational that gives the value also gives the + * derivative with no subtraction at all, because the `2 / sqrt(pi)` cancels + * analytically: + * + * erfcx(x) = (INV_SQRT_PI + u * C(u)) / x, u = 1 / x^2 + * d/dx erfcx(x) = 2 * u * C(u) + * + * The existing fvar and var suites compare against finite differences, whose + * tolerance is far looser than this error, so they pass either way. These + * rows are fixed references instead. + * + * References were evaluated at 60 significant digits from + * `2 * x * exp(x^2) * erfc(x) - 2 / sqrt(pi)` and rounded to double. + */ + +namespace { + +struct deriv_ref { + double x; + double d1; +}; + +// The suggested form measures 0.1 to 31 ulp across these points. 1e-13 +// relative leaves room for that and still rejects the difference form, which +// is wrong by 6.8e+03 ulp at x = 100 and more beyond. +const deriv_ref TAIL[] = {{4.0, -0.032383506095021455}, + {6.0, -0.015060353489052321}, + {10.0, -0.0055593122190608565}, + {100.0, -5.6410497625993184e-05}, + {1000.0, -5.641887372654967e-07}, + {10000.0, -5.641895750849127e-09}}; + +// Below 4 the difference form has no cancellation and stays in use. +const deriv_ref INTERIOR[] = {{-1.0, -11.14633932862008}, + {-5.0, -1440097986747.7388}}; + +constexpr double TOL = 1e-13; + +template +void expect_rows(const deriv_ref* rows, int n, const F& derivative, + const char* what) { + for (int i = 0; i < n; ++i) { + const double got = derivative(rows[i].x); + EXPECT_LT(std::fabs(got / rows[i].d1 - 1.0), TOL) + << what << " derivative lost precision at x = " << rows[i].x; + } +} + +double fwd_derivative(double x) { + stan::math::fvar xv(x, 1.0); + return stan::math::erfcx(xv).d_; +} + +double rev_derivative(double x) { + stan::math::var xv = x; + stan::math::var y = stan::math::erfcx(xv); + y.grad(); + const double d = xv.adj(); + stan::math::set_zero_all_adjoints(); + return d; +} + +} // namespace + +TEST(MathFunctions, erfcx_derivative_tail_fwd) { + expect_rows(TAIL, 6, fwd_derivative, "fvar tail"); +} + +TEST(MathFunctions, erfcx_derivative_tail_rev) { + expect_rows(TAIL, 6, rev_derivative, "var tail"); +} + +TEST(MathFunctions, erfcx_derivative_interior_fwd) { + expect_rows(INTERIOR, 2, fwd_derivative, "fvar interior"); +} + +TEST(MathFunctions, erfcx_derivative_interior_rev) { + expect_rows(INTERIOR, 2, rev_derivative, "var interior"); +} + +/** + * The two forms must agree where they meet, so the branch at 4 introduces no + * step in the derivative. + */ +TEST(MathFunctions, erfcx_derivative_no_step_at_four) { + const double below = fwd_derivative(std::nextafter(4.0, 0.0)); + const double at = fwd_derivative(4.0); + EXPECT_LT(std::fabs(below / at - 1.0), 1e-12); +} From baebb7ec277f044f79e0a8f1f4034b59c6940f52 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 12:59:25 +0300 Subject: [PATCH 12/24] Take the erfcx derivative from the tail rational to avoid cancellation --- stan/math/fwd/fun/erfcx.hpp | 8 ++-- stan/math/prim/fun/erfcx.hpp | 74 +++++++++++++++++++++++++++++------- stan/math/rev/fun/erfcx.hpp | 17 ++++++--- 3 files changed, 78 insertions(+), 21 deletions(-) diff --git a/stan/math/fwd/fun/erfcx.hpp b/stan/math/fwd/fun/erfcx.hpp index 4e1328e2d52..06e8ae5ea56 100644 --- a/stan/math/fwd/fun/erfcx.hpp +++ b/stan/math/fwd/fun/erfcx.hpp @@ -13,8 +13,10 @@ namespace math { /** * Return the scaled complementary error function of the argument. * - * The derivative `2 * x * erfcx(x) - 2 / sqrt(pi)` reuses the value, so the - * tangent is as accurate as the value itself. + * The derivative comes from `internal::erfcx_derivative`. Below 4 it is + * `2 * x * erfcx(x) - 2 / sqrt(pi)`, which reuses the value. At 4 and above + * that difference cancels, so the tail rational supplies the derivative + * directly. * * @tparam T inner type of the fvar * @param x argument @@ -23,7 +25,7 @@ namespace math { template inline fvar erfcx(const fvar& x) { T v = erfcx(x.val_); - return fvar(v, x.d_ * (2.0 * x.val_ * v - TWO_OVER_SQRT_PI)); + return fvar(v, x.d_ * internal::erfcx_derivative(x.val_, v)); } } // namespace math diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index adb04216726..59e361f1a10 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -11,6 +11,36 @@ namespace math { namespace internal { +/** + * Correction factor of the Cody (1969) third-interval rational: + * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. + * + * Split out from the value so that the derivative can reuse it. See + * `erfcx_derivative`. + * + * @tparam T scalar type + * @param u inverse square of the argument, `0 <= u <= 1/16` + * @return `P(u) / Q(u)` + */ +template +inline T erfcx_tail_correction(const T& u) { + static constexpr double p[] + = {0.000658749161529837803157, 0.0160837851487422766278, + 0.125781726111229246204, 0.360344899949804439429, + 0.305326634961232344035, 0.0163153871373020978498}; + static constexpr double q[] + = {-0.00233520497626869185443, -0.0605183413124413191178, + -0.527905102951428412248, -1.87295284992346047209, + -2.56852019228982242072, -1.0}; + T numerator = p[5]; + T denominator = q[5]; + for (int i = 4; i >= 0; --i) { + numerator = p[i] + u * numerator; + denominator = q[i] + u * denominator; + } + return numerator / denominator; +} + /** * Cody (1969) third-interval rational, valid for `x >= 4`. * @@ -23,19 +53,37 @@ namespace internal { */ inline double erfcx_cody_tail(double x) { const double u = 1.0 / (x * x); - double p = 0.0163153871373020978498; - p = 0.305326634961232344035 + u * p; - p = 0.360344899949804439429 + u * p; - p = 0.125781726111229246204 + u * p; - p = 0.0160837851487422766278 + u * p; - p = 0.000658749161529837803157 + u * p; - double q = -1.0; - q = -2.56852019228982242072 + u * q; - q = -1.87295284992346047209 + u * q; - q = -0.527905102951428412248 + u * q; - q = -0.0605183413124413191178 + u * q; - q = -0.00233520497626869185443 + u * q; - return (INV_SQRT_PI + (p / q) * u) / x; + return (INV_SQRT_PI + u * erfcx_tail_correction(u)) / x; +} + +/** + * Derivative of `erfcx`, `2 * x * erfcx(x) - 2 / sqrt(pi)`. + * + * That difference cancels for large `x`: both terms approach + * `2 / sqrt(pi)` while the result decays like `1 / (sqrt(pi) * x^2)`. + * Measured against a 50-digit reference, the difference form gives 6.1 ulp + * at `x = 4` and 2.55e+11 ulp at `x = 1e6`. + * + * For `x >= 4` the tail rational gives the derivative with no subtraction, + * because the constant cancels analytically: + * + * `2 * x * (INV_SQRT_PI + u * C(u)) / x - 2 / sqrt(pi) = 2 * u * C(u)` + * + * since `2 * INV_SQRT_PI` is `2 / sqrt(pi)`. That form measures 0.1 to + * 31 ulp over the same range. + * + * @tparam T scalar type + * @param x argument + * @param value `erfcx(x)` + * @return derivative of `erfcx` at `x` + */ +template +inline T erfcx_derivative(const T& x, const T& value) { + if (x >= 4.0) { + const T u = 1.0 / (x * x); + return 2.0 * u * erfcx_tail_correction(u); + } + return 2.0 * x * value - TWO_OVER_SQRT_PI; } /** diff --git a/stan/math/rev/fun/erfcx.hpp b/stan/math/rev/fun/erfcx.hpp index 5228be7af19..ab9c407a86f 100644 --- a/stan/math/rev/fun/erfcx.hpp +++ b/stan/math/rev/fun/erfcx.hpp @@ -19,7 +19,9 @@ namespace math { \frac{2}{\sqrt{\pi}}\f$ * * which reuses the function value, so no extra `exp` or `erfc` evaluation is - * needed and the derivative inherits the value's accuracy in the tails. + * needed. That difference cancels for `x >= 4`, so `internal::erfcx_derivative` + * takes the derivative from the tail rational there instead. Without that the + * error reaches 2.55e+11 ulp at `x = 1e6`. * \f[ \mbox{erfcx}(x) = @@ -44,7 +46,7 @@ namespace math { inline var erfcx(const var& a) { double val = erfcx(a.val()); return make_callback_var(val, [a, val](auto& vi) mutable { - a.adj() += vi.adj() * (2.0 * a.val() * val - TWO_OVER_SQRT_PI); + a.adj() += vi.adj() * internal::erfcx_derivative(a.val(), val); }); } @@ -59,9 +61,14 @@ template * = nullptr> inline auto erfcx(const var_value& a) { auto val = to_arena(erfcx(a.val())); return make_callback_var(val, [a, val](auto& vi) mutable { - a.adj().array() - += vi.adj().array() - * (2.0 * a.val().array() * val.array() - TWO_OVER_SQRT_PI); + a.adj().array() += vi.adj().array() + * a.val() + .array() + .binaryExpr(val.array(), + [](double x, double v) { + return internal::erfcx_derivative( + x, v); + }); }); } From ce4dea1efcc7f4f175de09fbe460348c49c68b16 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 13:11:38 +0300 Subject: [PATCH 13/24] Split the erfcx small-branch polynomial and merge the reverse-mode overloads --- stan/math/prim/fun/erfcx.hpp | 56 ++++++++++++++++++++++-------------- stan/math/rev/fun/erfcx.hpp | 43 ++++++++++++--------------- 2 files changed, 54 insertions(+), 45 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 59e361f1a10..c6830477084 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -122,8 +122,13 @@ inline double erfcx_cody_middle(double y) { * A Chebyshev-economized expansion of `erfcx` about zero, so it needs no * library call and no branch on the sign of `x`. Its low-order coefficients * reproduce the Maclaurin series of `erfcx` exactly (`1`, `-2/sqrt(pi)`, - * `1`, ...), which is a useful check that the fit is right. Measured at - * 2.2 ulp over the whole interval. + * `1`, ...), which is a useful check that the fit is right. `even[0]`, + * `even[1]` and `odd[0]` are those three terms. + * + * Measured at 2.2 ulp over the interval. Evaluating the same coefficients + * as one degree-18 Horner chain gives 1.5 ulp instead, but runs 1.7 times + * slower, because that chain is 18 dependent operations. The whole function + * is 7.0 ulp, set by the middle interval, so the 0.7 ulp costs nothing. * * The plain Maclaurin series is only usable to about `|x| = 0.125`; four * further economized terms extend it to 0.46875 at no measurable cost, @@ -133,25 +138,34 @@ inline double erfcx_cody_middle(double y) { * @return scaled complementary error function */ inline double erfcx_small(double x) { - double p = 3.05977060678449757e-06; - p = -9.35890030086883823e-06 + x * p; - p = 2.46655529768908249e-05 + x * p; - p = -7.08163358203131886e-05 + x * p; - p = 1.98445679338826757e-04 + x * p; - p = -5.34506929034156810e-04 + x * p; - p = 1.38888415444527033e-03 + x * p; - p = -3.47359067853470795e-03 + x * p; - p = 8.33333374332981443e-03 + x * p; - p = -1.91048337772546720e-02 + x * p; - p = 4.16666666458337179e-02 + x * p; - p = -8.59717459974174147e-02 + x * p; - p = 1.66666666667239644e-01 + x * p; - p = -3.00901111227312890e-01 + x * p; - p = 4.99999999999992839e-01 + x * p; - p = -7.52252778063651983e-01 + x * p; - p = 1.0 + x * p; - p = -1.12837916709551256 + x * p; - return 1.0 + x * p; + // Split into the even and odd powers of x, so the two Horner chains run + // independently. A single degree-18 chain is 18 dependent operations; two + // chains of 9 halve that latency. + static constexpr double even[] + = {1.0, + 1.0, + 4.99999999999992839e-01, + 1.66666666667239644e-01, + 4.16666666458337179e-02, + 8.33333374332981443e-03, + 1.38888415444527033e-03, + 1.98445679338826757e-04, + 2.46655529768908249e-05, + 3.05977060678449757e-06}; + static constexpr double odd[] + = {-1.12837916709551256, -7.52252778063651983e-01, + -3.00901111227312890e-01, -8.59717459974174147e-02, + -1.91048337772546720e-02, -3.47359067853470795e-03, + -5.34506929034156810e-04, -7.08163358203131886e-05, + -9.35890030086883823e-06}; + const double x2 = x * x; + double e = even[9]; + double o = odd[8]; + for (int i = 8; i >= 1; --i) { + e = even[i] + x2 * e; + o = odd[i - 1] + x2 * o; + } + return even[0] + x2 * e + x * o; } } // namespace internal diff --git a/stan/math/rev/fun/erfcx.hpp b/stan/math/rev/fun/erfcx.hpp index ab9c407a86f..e911d84aa23 100644 --- a/stan/math/rev/fun/erfcx.hpp +++ b/stan/math/rev/fun/erfcx.hpp @@ -3,8 +3,10 @@ #include #include +#include #include #include +#include #include namespace stan { @@ -40,35 +42,28 @@ namespace math { \end{cases} \f] * - * @param a The variable. - * @return Scaled complementary error function applied to the variable. - */ -inline var erfcx(const var& a) { - double val = erfcx(a.val()); - return make_callback_var(val, [a, val](auto& vi) mutable { - a.adj() += vi.adj() * internal::erfcx_derivative(a.val(), val); - }); -} - -/** - * The scaled complementary error function for matrix variables. + * One overload covers `var` and `var_value`. For a matrix the + * derivative stays a lazy `binaryExpr`, so the reverse pass allocates + * nothing. * - * @tparam T a matrix type + * @tparam T a `var` or a `var_value` of a matrix type * @param a The variable. - * @return Scaled complementary error function applied elementwise. + * @return Scaled complementary error function applied to the variable, + * elementwise for a matrix. */ -template * = nullptr> -inline auto erfcx(const var_value& a) { +template < + typename T, require_var_t* = nullptr, + require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> +inline auto erfcx(T&& a) { auto val = to_arena(erfcx(a.val())); return make_callback_var(val, [a, val](auto& vi) mutable { - a.adj().array() += vi.adj().array() - * a.val() - .array() - .binaryExpr(val.array(), - [](double x, double v) { - return internal::erfcx_derivative( - x, v); - }); + // apply_scalar_binary returns a lazy binaryExpr for Eigen inputs, so + // this allocates nothing, and calls the functor directly for scalars. + const auto& deriv = apply_scalar_binary( + [](double x, double v) { return internal::erfcx_derivative(x, v); }, + a.val(), val); + as_array_or_scalar(a.adj()) + += as_array_or_scalar(vi.adj()) * as_array_or_scalar(deriv); }); } From b9be613dc6886915f177ef3b8dbe9591cf29ed49 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 18 Sep 2026 09:14:54 -0400 Subject: [PATCH 14/24] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/erfcx.hpp | 21 +++++++++---------- .../math/mix/fun/erfcx_derivative_test.cpp | 14 ++++++------- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index c6830477084..1331f2465bb 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -141,17 +141,16 @@ inline double erfcx_small(double x) { // Split into the even and odd powers of x, so the two Horner chains run // independently. A single degree-18 chain is 18 dependent operations; two // chains of 9 halve that latency. - static constexpr double even[] - = {1.0, - 1.0, - 4.99999999999992839e-01, - 1.66666666667239644e-01, - 4.16666666458337179e-02, - 8.33333374332981443e-03, - 1.38888415444527033e-03, - 1.98445679338826757e-04, - 2.46655529768908249e-05, - 3.05977060678449757e-06}; + static constexpr double even[] = {1.0, + 1.0, + 4.99999999999992839e-01, + 1.66666666667239644e-01, + 4.16666666458337179e-02, + 8.33333374332981443e-03, + 1.38888415444527033e-03, + 1.98445679338826757e-04, + 2.46655529768908249e-05, + 3.05977060678449757e-06}; static constexpr double odd[] = {-1.12837916709551256, -7.52252778063651983e-01, -3.00901111227312890e-01, -8.59717459974174147e-02, diff --git a/test/unit/math/mix/fun/erfcx_derivative_test.cpp b/test/unit/math/mix/fun/erfcx_derivative_test.cpp index 31bdfdc4cb3..c117bdb6ca7 100644 --- a/test/unit/math/mix/fun/erfcx_derivative_test.cpp +++ b/test/unit/math/mix/fun/erfcx_derivative_test.cpp @@ -36,16 +36,14 @@ struct deriv_ref { // The suggested form measures 0.1 to 31 ulp across these points. 1e-13 // relative leaves room for that and still rejects the difference form, which // is wrong by 6.8e+03 ulp at x = 100 and more beyond. -const deriv_ref TAIL[] = {{4.0, -0.032383506095021455}, - {6.0, -0.015060353489052321}, - {10.0, -0.0055593122190608565}, - {100.0, -5.6410497625993184e-05}, - {1000.0, -5.641887372654967e-07}, - {10000.0, -5.641895750849127e-09}}; +const deriv_ref TAIL[] + = {{4.0, -0.032383506095021455}, {6.0, -0.015060353489052321}, + {10.0, -0.0055593122190608565}, {100.0, -5.6410497625993184e-05}, + {1000.0, -5.641887372654967e-07}, {10000.0, -5.641895750849127e-09}}; // Below 4 the difference form has no cancellation and stays in use. -const deriv_ref INTERIOR[] = {{-1.0, -11.14633932862008}, - {-5.0, -1440097986747.7388}}; +const deriv_ref INTERIOR[] + = {{-1.0, -11.14633932862008}, {-5.0, -1440097986747.7388}}; constexpr double TOL = 1e-13; From dadf5374dab6e51c6daa736e702aeb6f3c026854 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 16:24:17 +0300 Subject: [PATCH 15/24] Add OpenCL tail tests for erfcx including fixed-reference derivatives --- test/unit/math/opencl/rev/erfcx_test.cpp | 49 ++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/test/unit/math/opencl/rev/erfcx_test.cpp b/test/unit/math/opencl/rev/erfcx_test.cpp index ddbba02a8c4..787bce17507 100644 --- a/test/unit/math/opencl/rev/erfcx_test.cpp +++ b/test/unit/math/opencl/rev/erfcx_test.cpp @@ -12,6 +12,55 @@ TEST(OpenCLerfcx, prim_rev_values_small) { stan::math::test::compare_cpu_opencl_prim_rev(erfcx_functor, a); } +/** + * The tail, `x >= 4`, where the value uses the Cody rational and the + * derivative uses its correction factor directly. + * + * Every other case here stays below 2.6, so before this the tail branch was + * never executed on the device, for the value or for the derivative. + */ +TEST(OpenCLerfcx, prim_rev_values_tail) { + Eigen::VectorXd a(8); + a << 3.9, 4.0, 4.1, 6.0, 10.0, 100.0, 1000.0, 10000.0; + stan::math::test::compare_cpu_opencl_prim_rev(erfcx_functor, a); +} + +/** + * Device derivative in the tail against fixed references. + * + * `compare_cpu_opencl_prim_rev` cannot check this. It uses + * `EXPECT_NEAR_REL`, whose relative tolerance floors its denominator at 1, + * so for a derivative of magnitude 5.6e-09 even a 2.4e-08 relative error is + * a 1.4e-16 absolute difference and passes. Verified: the differential test + * above passes against both the old and the new device derivative. + * + * References are `2 * x * erfcx(x) - 2 / sqrt(pi)` at 60 significant + * digits, rounded to double. They are the same values the CPU suite uses in + * `test/unit/math/mix/fun/erfcx_derivative_test.cpp`. + */ +TEST(OpenCLerfcx, rev_tail_derivative_against_references) { + const int N = 6; + Eigen::VectorXd x(N); + x << 4.0, 6.0, 10.0, 100.0, 1000.0, 10000.0; + Eigen::VectorXd expected(N); + expected << -0.032383506095021455, -0.015060353489052321, + -0.0055593122190608565, -5.6410497625993184e-05, + -5.641887372654967e-07, -5.641895750849127e-09; + + stan::math::var_value> x_cl( + stan::math::to_matrix_cl(x)); + auto y = stan::math::erfcx(x_cl); + stan::math::var total = stan::math::sum(y); + total.grad(); + + const Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj()); + for (int i = 0; i < N; ++i) { + EXPECT_LT(std::fabs(adj[i] / expected[i] - 1.0), 1e-12) + << "device derivative lost precision at x = " << x[i]; + } + stan::math::set_zero_all_adjoints(); +} + TEST(OpenCLerfcx, prim_rev_size_0) { int N = 0; From 3453e8afc76ab02bb958f7343bd2975d0028a29c Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 16:24:23 +0300 Subject: [PATCH 16/24] Take the OpenCL erfcx derivative from the tail rational to avoid cancellation --- .../kernel_generator/elt_function_cl.hpp | 2 + .../opencl/kernels/device_functions/erfcx.hpp | 52 ++++++++++++++++--- stan/math/opencl/rev/erfcx.hpp | 13 +++-- 3 files changed, 56 insertions(+), 11 deletions(-) diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 6ca3be497e6..df349a7b483 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -353,6 +353,8 @@ ADD_BINARY_FUNCTION_WITH_INCLUDES(ldexp) ADD_BINARY_FUNCTION_WITH_INCLUDES(pow) ADD_BINARY_FUNCTION_WITH_INCLUDES(copysign) +ADD_BINARY_FUNCTION_WITH_INCLUDES( + erfcx_derivative, stan::math::opencl_kernels::erfcx_device_function) ADD_BINARY_FUNCTION_WITH_INCLUDES( beta, stan::math::opencl_kernels::beta_device_function) ADD_BINARY_FUNCTION_WITH_INCLUDES( diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 86481a0279a..8c61159502a 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -17,13 +17,16 @@ static constexpr const char* erfcx_device_function // \endcond /** \ingroup opencl_kernels * - * Cody (1969) third-interval rational, for x >= 4. + * Correction factor of the Cody (1969) third-interval rational: + * erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x, u = 1 / x^2. * - * @param x argument - * @return scaled complementary error function + * Split out so the derivative can reuse it. See + * erfcx_tail_derivative. + * + * @param u inverse square of the argument, 0 <= u <= 1/16 + * @return P(u) / Q(u) */ - double erfcx_cody_tail(double x) { - double u = 1.0 / (x * x); + double erfcx_tail_correction(double u) { double p = 0.0163153871373020978498; p = 0.305326634961232344035 + u * p; p = 0.360344899949804439429 + u * p; @@ -36,7 +39,44 @@ static constexpr const char* erfcx_device_function q = -0.527905102951428412248 + u * q; q = -0.0605183413124413191178 + u * q; q = -0.00233520497626869185443 + u * q; - return (M_2_SQRTPI * 0.5 + (p / q) * u) / x; + return p / q; + } + + /** \ingroup opencl_kernels + * + * Cody (1969) third-interval rational, for x >= 4. + * + * @param x argument + * @return scaled complementary error function + */ + double erfcx_cody_tail(double x) { + double u = 1.0 / (x * x); + return (M_2_SQRTPI * 0.5 + u * erfcx_tail_correction(u)) / x; + } + + /** \ingroup opencl_kernels + * + * Derivative of erfcx, 2 * x * erfcx(x) - 2 / sqrt(pi). + * + * That difference cancels for large x: both terms approach + * 2 / sqrt(pi) while the result decays like 1 / (sqrt(pi) * x^2). + * The difference form reaches 2.55e+11 ulp at x = 1e6. For + * x >= 4 the constant cancels analytically against the leading + * term of the tail rational, leaving 2 * u * C(u). + * + * Takes the value as an argument so the reverse pass does not + * evaluate erfcx twice. + * + * @param x argument + * @param value erfcx(x) + * @return derivative of erfcx at x + */ + double erfcx_derivative(double x, double value) { + if (x >= 4.0) { + double u = 1.0 / (x * x); + return 2.0 * u * erfcx_tail_correction(u); + } + return 2.0 * x * value - M_2_SQRTPI; } /** \ingroup opencl_kernels diff --git a/stan/math/opencl/rev/erfcx.hpp b/stan/math/opencl/rev/erfcx.hpp index bb77a027b97..fc43e8c94f1 100644 --- a/stan/math/opencl/rev/erfcx.hpp +++ b/stan/math/opencl/rev/erfcx.hpp @@ -13,8 +13,12 @@ namespace math { * Returns the elementwise `erfcx()` of a var_value>. * * The derivative `2 * x * erfcx(x) - 2 / sqrt(pi)` reuses the function - * value, so no second `exp` or `erfc` evaluation is needed and the - * derivative inherits the value's accuracy in the tails. + * value, so no second `exp` or `erfc` evaluation is needed. That difference + * cancels for `x >= 4`, where both terms approach `2 / sqrt(pi)` while the + * result decays like `1 / (sqrt(pi) * x^2)`; the device function + * `erfcx_derivative` takes the derivative from the tail rational there + * instead. Without that the error reaches 2.55e+11 ulp at `x = 1e6`. The + * CPU implementation branches at the same point. * * @param A argument * @return Elementwise `erfcx()` of the input. @@ -24,9 +28,8 @@ template > erfcx(const var_value& A) { return make_callback_var( erfcx(A.val()), [A](vari_value>& res) mutable { - A.adj() += elt_multiply( - res.adj(), elt_multiply(2.0, elt_multiply(A.val(), res.val())) - - TWO_OVER_SQRT_PI); + A.adj() += elt_multiply(res.adj(), + erfcx_derivative(A.val(), res.val())); }); } From 3b1ded96147ff4275f50d99cc29fa656076f4f02 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 18 Sep 2026 09:52:09 -0400 Subject: [PATCH 17/24] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/rev/erfcx.hpp | 4 ++-- test/unit/math/opencl/rev/erfcx_test.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/opencl/rev/erfcx.hpp b/stan/math/opencl/rev/erfcx.hpp index fc43e8c94f1..e5fdfeb097a 100644 --- a/stan/math/opencl/rev/erfcx.hpp +++ b/stan/math/opencl/rev/erfcx.hpp @@ -28,8 +28,8 @@ template > erfcx(const var_value& A) { return make_callback_var( erfcx(A.val()), [A](vari_value>& res) mutable { - A.adj() += elt_multiply(res.adj(), - erfcx_derivative(A.val(), res.val())); + A.adj() + += elt_multiply(res.adj(), erfcx_derivative(A.val(), res.val())); }); } diff --git a/test/unit/math/opencl/rev/erfcx_test.cpp b/test/unit/math/opencl/rev/erfcx_test.cpp index 787bce17507..d1cf93dd4ed 100644 --- a/test/unit/math/opencl/rev/erfcx_test.cpp +++ b/test/unit/math/opencl/rev/erfcx_test.cpp @@ -44,8 +44,8 @@ TEST(OpenCLerfcx, rev_tail_derivative_against_references) { x << 4.0, 6.0, 10.0, 100.0, 1000.0, 10000.0; Eigen::VectorXd expected(N); expected << -0.032383506095021455, -0.015060353489052321, - -0.0055593122190608565, -5.6410497625993184e-05, - -5.641887372654967e-07, -5.641895750849127e-09; + -0.0055593122190608565, -5.6410497625993184e-05, -5.641887372654967e-07, + -5.641895750849127e-09; stan::math::var_value> x_cl( stan::math::to_matrix_cl(x)); From 8663eccb36788a3ddc2a2713d5feda9a620e3baa Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 17:25:08 +0300 Subject: [PATCH 18/24] Apply clang-format 10 to the OpenCL erfcx files --- stan/math/opencl/rev/erfcx.hpp | 4 ++-- test/unit/math/opencl/rev/erfcx_test.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/opencl/rev/erfcx.hpp b/stan/math/opencl/rev/erfcx.hpp index fc43e8c94f1..e5fdfeb097a 100644 --- a/stan/math/opencl/rev/erfcx.hpp +++ b/stan/math/opencl/rev/erfcx.hpp @@ -28,8 +28,8 @@ template > erfcx(const var_value& A) { return make_callback_var( erfcx(A.val()), [A](vari_value>& res) mutable { - A.adj() += elt_multiply(res.adj(), - erfcx_derivative(A.val(), res.val())); + A.adj() + += elt_multiply(res.adj(), erfcx_derivative(A.val(), res.val())); }); } diff --git a/test/unit/math/opencl/rev/erfcx_test.cpp b/test/unit/math/opencl/rev/erfcx_test.cpp index 787bce17507..d1cf93dd4ed 100644 --- a/test/unit/math/opencl/rev/erfcx_test.cpp +++ b/test/unit/math/opencl/rev/erfcx_test.cpp @@ -44,8 +44,8 @@ TEST(OpenCLerfcx, rev_tail_derivative_against_references) { x << 4.0, 6.0, 10.0, 100.0, 1000.0, 10000.0; Eigen::VectorXd expected(N); expected << -0.032383506095021455, -0.015060353489052321, - -0.0055593122190608565, -5.6410497625993184e-05, - -5.641887372654967e-07, -5.641895750849127e-09; + -0.0055593122190608565, -5.6410497625993184e-05, -5.641887372654967e-07, + -5.641895750849127e-09; stan::math::var_value> x_cl( stan::math::to_matrix_cl(x)); From 2933f2897879aea39f7e9ec9f1963500ea13f516 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 21:47:25 +0300 Subject: [PATCH 19/24] Apply batched suggestions from code review Co-authored-by: Steve Bronder --- .../opencl/kernels/device_functions/erfcx.hpp | 142 ++++++++++++------ stan/math/prim/fun/erfcx.hpp | 57 +++---- 2 files changed, 117 insertions(+), 82 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 8c61159502a..e094fa743c5 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -88,25 +88,53 @@ static constexpr const char* erfcx_device_function * @param y argument * @return scaled complementary error function */ - double erfcx_cody_middle(double y) { - double p = 2.15311535474403846e-8 * y; - p = (p + 5.64188496988670089e-1) * y; - p = (p + 8.88314979438837594) * y; - p = (p + 66.1191906371416295) * y; - p = (p + 298.635138197400131) * y; - p = (p + 881.952221241769090) * y; - p = (p + 1712.04761263407058) * y; - p = (p + 2051.07837782607147) * y; - double q = y; - q = (q + 15.7449261107098347) * y; - q = (q + 117.693950891312499) * y; - q = (q + 537.181101862009858) * y; - q = (q + 1621.38957456669019) * y; - q = (q + 3290.79923573345963) * y; - q = (q + 4362.61909014324716) * y; - q = (q + 3439.36767414372164) * y; - return (p + 1230.33935479799725) / (q + 1230.33935480374942); - } +inline double erfcx_cody_middle(double y) { + constexpr std::array p{ + 1230.33935479799725, + 2051.07837782607147, + 1712.04761263407058, + 881.952221241769090, + 298.635138197400131, + 66.1191906371416295, + 8.88314979438837594, + 5.64188496988670089e-1 + }; + + constexpr std::array q{ + 1230.33935480374942, + 3439.36767414372164, + 4362.61909014324716, + 3290.79923573345963, + 1621.38957456669019, + 537.181101862009858, + 117.693950891312499, + 15.7449261107098347, + }; + + const double y2 = y * y; + + const double p01 = std::fma(p[1], y, p[0]); + const double p23 = std::fma(p[3], y, p[2]); + const double p45 = std::fma(p[5], y, p[4]); + const double p67 = std::fma(p[7], y, p[6]); + + const double q01 = std::fma(q[1], y, q[0]); + const double q23 = std::fma(q[3], y, q[2]); + const double q45 = std::fma(q[5], y, q[4]); + const double q67 = std::fma(q[7], y, q[6]); + + double num = std::fma(2.15311535474403846e-8, y2, p67); + num = std::fma(num, y2, p45); + num = std::fma(num, y2, p23); + num = std::fma(num, y2, p01); + + double den = y2 + q67; + den = std::fma(den, y2, q45); + den = std::fma(den, y2, q23); + den = std::fma(den, y2, q01); + + return num / den; +} /** \ingroup opencl_kernels * @@ -117,27 +145,52 @@ static constexpr const char* erfcx_device_function * @param x argument * @return scaled complementary error function */ - double erfcx_small(double x) { - double p = 3.05977060678449757e-06; - p = -9.35890030086883823e-06 + x * p; - p = 2.46655529768908249e-05 + x * p; - p = -7.08163358203131886e-05 + x * p; - p = 1.98445679338826757e-04 + x * p; - p = -5.34506929034156810e-04 + x * p; - p = 1.38888415444527033e-03 + x * p; - p = -3.47359067853470795e-03 + x * p; - p = 8.33333374332981443e-03 + x * p; - p = -1.91048337772546720e-02 + x * p; - p = 4.16666666458337179e-02 + x * p; - p = -8.59717459974174147e-02 + x * p; - p = 1.66666666667239644e-01 + x * p; - p = -3.00901111227312890e-01 + x * p; - p = 4.99999999999992839e-01 + x * p; - p = -7.52252778063651983e-01 + x * p; - p = 1.0 + x * p; - p = -1.12837916709551256 + x * p; - return 1.0 + x * p; - } +inline double erfcx_small(double x) { + constexpr std::array p{ + 1.0, + -1.12837916709551256, + 1.0, + -7.52252778063651983e-01, + 4.99999999999992839e-01, + -3.00901111227312890e-01, + 1.66666666667239644e-01, + -8.59717459974174147e-02, + 4.16666666458337179e-02, + -1.91048337772546720e-02, + 8.33333374332981443e-03, + -3.47359067853470795e-03, + 1.38888415444527033e-03, + -5.34506929034156810e-04, + 1.98445679338826757e-04, + -7.08163358203131886e-05, + 2.46655529768908249e-05, + -9.35890030086883823e-06 + }; + + const double x2 = x * x; + + // Independent linear terms. + const double p01 = std::fma(p[1], x, p[0]); + const double p23 = std::fma(p[3], x, p[2]); + const double p45 = std::fma(p[5], x, p[4]); + const double p67 = std::fma(p[7], x, p[6]); + const double p89 = std::fma(p[9], x, p[8]); + const double p1011 = std::fma(p[11], x, p[10]); + const double p1213 = std::fma(p[13], x, p[12]); + const double p1415 = std::fma(p[15], x, p[14]); + const double p1617 = std::fma(p[17], x, p[16]); + + // Combine using Horner's method in x². + double result = std::fma(3.05977060678449757e-06, x2, p1617); + result = std::fma(result, x2, p1415); + result = std::fma(result, x2, p1213); + result = std::fma(result, x2, p1011); + result = std::fma(result, x2, p89); + result = std::fma(result, x2, p67); + result = std::fma(result, x2, p45); + result = std::fma(result, x2, p23); + return std::fma(result, x2, p01); +} /** \ingroup opencl_kernels * @@ -165,14 +218,11 @@ static constexpr const char* erfcx_device_function double erfcx(double x) { if (x >= 4.0) { return erfcx_cody_tail(x); - } - if (x >= 0.46875) { + } else if (x >= 0.46875) { return erfcx_cody_middle(x); - } - if (x > -0.46875) { + } else if (x > -0.46875) { return erfcx_small(x); - } - if (x < -27.0) { + } else if (x < -27.0) { return INFINITY; } double h = x * x; diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 1331f2465bb..f3f01e50169 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -24,19 +24,17 @@ namespace internal { */ template inline T erfcx_tail_correction(const T& u) { - static constexpr double p[] - = {0.000658749161529837803157, 0.0160837851487422766278, + static constexpr std::array p{0.000658749161529837803157, 0.0160837851487422766278, 0.125781726111229246204, 0.360344899949804439429, - 0.305326634961232344035, 0.0163153871373020978498}; - static constexpr double q[] - = {-0.00233520497626869185443, -0.0605183413124413191178, + 0.305326634961232344035}; + static constexpr std::array q{-0.00233520497626869185443, -0.0605183413124413191178, -0.527905102951428412248, -1.87295284992346047209, - -2.56852019228982242072, -1.0}; - T numerator = p[5]; - T denominator = q[5]; + -2.56852019228982242072}; + T numerator = 0.0163153871373020978498; + T denominator = -1.0; for (int i = 4; i >= 0; --i) { - numerator = p[i] + u * numerator; - denominator = q[i] + u * denominator; + numerator = numerator * u + p[i]; + denominator = denominator * u + q[i]; } return numerator / denominator; } @@ -141,30 +139,27 @@ inline double erfcx_small(double x) { // Split into the even and odd powers of x, so the two Horner chains run // independently. A single degree-18 chain is 18 dependent operations; two // chains of 9 halve that latency. - static constexpr double even[] = {1.0, - 1.0, + static constexpr std::array even = {1.0, 4.99999999999992839e-01, 1.66666666667239644e-01, 4.16666666458337179e-02, 8.33333374332981443e-03, 1.38888415444527033e-03, 1.98445679338826757e-04, - 2.46655529768908249e-05, - 3.05977060678449757e-06}; - static constexpr double odd[] + 2.46655529768908249e-05}; + static constexpr std::array odd = {-1.12837916709551256, -7.52252778063651983e-01, -3.00901111227312890e-01, -8.59717459974174147e-02, -1.91048337772546720e-02, -3.47359067853470795e-03, - -5.34506929034156810e-04, -7.08163358203131886e-05, - -9.35890030086883823e-06}; + -5.34506929034156810e-04, -7.08163358203131886e-05}; const double x2 = x * x; - double e = even[9]; - double o = odd[8]; - for (int i = 8; i >= 1; --i) { + double e = 3.05977060678449757e-06; + double o = -9.35890030086883823e-06; + for (int i = 8; i > 0; --i) { e = even[i] + x2 * e; - o = odd[i - 1] + x2 * o; + o = odd[i] + x2 * o; } - return even[0] + x2 * e + x * o; + return 1.0+ x2 * e + x * o; } } // namespace internal @@ -238,26 +233,16 @@ inline double erfcx(T&& xx) { const double x = static_cast(xx); constexpr double cody_min = 4.0; constexpr double middle_min = 0.46875; + constexpr double overflow_max = -27.0; if (x >= cody_min) { return internal::erfcx_cody_tail(x); - } - if (x >= middle_min) { + } else if (x >= middle_min) { return internal::erfcx_cody_middle(x); - } - if (x > -middle_min) { + } else if (x > -middle_min) { return internal::erfcx_small(x); - } - // erfcx(x) = 2 * exp(x * x) * (1 + o(1)) as x -> -infinity, which leaves - // the binary64 range below -26.63. Returning here also keeps the - // reflection below from forming inf * 0 on -infinity. - constexpr double overflow_max = -27.0; - if (x < overflow_max) { + } else if (x < overflow_max) { return INFTY; } - // fma(x, x, -h) is the exact rounding error of h = x * x, so the factor - // (1 + that) restores what exp(h) would otherwise lose. The OpenCL - // device function uses the same form; it cannot use a Dekker split - // because the compiler simplifies t - (t - x) to x. const double h = x * x; const double two_exp_x2 = 2.0 * std::exp(h) * (1.0 + std::fma(x, x, -h)); constexpr double reflect_min = -6.1; From 0b1094f3e2798884755334d01b26794ade09f728 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 18 Sep 2026 14:50:32 -0400 Subject: [PATCH 20/24] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../opencl/kernels/device_functions/erfcx.hpp | 173 ++++++++---------- stan/math/prim/fun/erfcx.hpp | 30 +-- 2 files changed, 96 insertions(+), 107 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index e094fa743c5..afe9b47d05a 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -88,53 +88,42 @@ static constexpr const char* erfcx_device_function * @param y argument * @return scaled complementary error function */ -inline double erfcx_cody_middle(double y) { - constexpr std::array p{ - 1230.33935479799725, - 2051.07837782607147, - 1712.04761263407058, - 881.952221241769090, - 298.635138197400131, - 66.1191906371416295, - 8.88314979438837594, - 5.64188496988670089e-1 - }; - - constexpr std::array q{ - 1230.33935480374942, - 3439.36767414372164, - 4362.61909014324716, - 3290.79923573345963, - 1621.38957456669019, - 537.181101862009858, - 117.693950891312499, - 15.7449261107098347, - }; - - const double y2 = y * y; - - const double p01 = std::fma(p[1], y, p[0]); - const double p23 = std::fma(p[3], y, p[2]); - const double p45 = std::fma(p[5], y, p[4]); - const double p67 = std::fma(p[7], y, p[6]); - - const double q01 = std::fma(q[1], y, q[0]); - const double q23 = std::fma(q[3], y, q[2]); - const double q45 = std::fma(q[5], y, q[4]); - const double q67 = std::fma(q[7], y, q[6]); - - double num = std::fma(2.15311535474403846e-8, y2, p67); - num = std::fma(num, y2, p45); - num = std::fma(num, y2, p23); - num = std::fma(num, y2, p01); - - double den = y2 + q67; - den = std::fma(den, y2, q45); - den = std::fma(den, y2, q23); - den = std::fma(den, y2, q01); - - return num / den; -} + inline double erfcx_cody_middle(double y) { + constexpr std::array p{1230.33935479799725, 2051.07837782607147, + 1712.04761263407058, 881.952221241769090, + 298.635138197400131, 66.1191906371416295, + 8.88314979438837594, 5.64188496988670089e-1}; + + constexpr std::array q{ + 1230.33935480374942, 3439.36767414372164, 4362.61909014324716, + 3290.79923573345963, 1621.38957456669019, 537.181101862009858, + 117.693950891312499, 15.7449261107098347, + }; + + const double y2 = y * y; + + const double p01 = std::fma(p[1], y, p[0]); + const double p23 = std::fma(p[3], y, p[2]); + const double p45 = std::fma(p[5], y, p[4]); + const double p67 = std::fma(p[7], y, p[6]); + + const double q01 = std::fma(q[1], y, q[0]); + const double q23 = std::fma(q[3], y, q[2]); + const double q45 = std::fma(q[5], y, q[4]); + const double q67 = std::fma(q[7], y, q[6]); + + double num = std::fma(2.15311535474403846e-8, y2, p67); + num = std::fma(num, y2, p45); + num = std::fma(num, y2, p23); + num = std::fma(num, y2, p01); + + double den = y2 + q67; + den = std::fma(den, y2, q45); + den = std::fma(den, y2, q23); + den = std::fma(den, y2, q01); + + return num / den; + } /** \ingroup opencl_kernels * @@ -145,52 +134,50 @@ inline double erfcx_cody_middle(double y) { * @param x argument * @return scaled complementary error function */ -inline double erfcx_small(double x) { - constexpr std::array p{ - 1.0, - -1.12837916709551256, - 1.0, - -7.52252778063651983e-01, - 4.99999999999992839e-01, - -3.00901111227312890e-01, - 1.66666666667239644e-01, - -8.59717459974174147e-02, - 4.16666666458337179e-02, - -1.91048337772546720e-02, - 8.33333374332981443e-03, - -3.47359067853470795e-03, - 1.38888415444527033e-03, - -5.34506929034156810e-04, - 1.98445679338826757e-04, - -7.08163358203131886e-05, - 2.46655529768908249e-05, - -9.35890030086883823e-06 - }; - - const double x2 = x * x; - - // Independent linear terms. - const double p01 = std::fma(p[1], x, p[0]); - const double p23 = std::fma(p[3], x, p[2]); - const double p45 = std::fma(p[5], x, p[4]); - const double p67 = std::fma(p[7], x, p[6]); - const double p89 = std::fma(p[9], x, p[8]); - const double p1011 = std::fma(p[11], x, p[10]); - const double p1213 = std::fma(p[13], x, p[12]); - const double p1415 = std::fma(p[15], x, p[14]); - const double p1617 = std::fma(p[17], x, p[16]); - - // Combine using Horner's method in x². - double result = std::fma(3.05977060678449757e-06, x2, p1617); - result = std::fma(result, x2, p1415); - result = std::fma(result, x2, p1213); - result = std::fma(result, x2, p1011); - result = std::fma(result, x2, p89); - result = std::fma(result, x2, p67); - result = std::fma(result, x2, p45); - result = std::fma(result, x2, p23); - return std::fma(result, x2, p01); -} + inline double erfcx_small(double x) { + constexpr std::array p{1.0, + -1.12837916709551256, + 1.0, + -7.52252778063651983e-01, + 4.99999999999992839e-01, + -3.00901111227312890e-01, + 1.66666666667239644e-01, + -8.59717459974174147e-02, + 4.16666666458337179e-02, + -1.91048337772546720e-02, + 8.33333374332981443e-03, + -3.47359067853470795e-03, + 1.38888415444527033e-03, + -5.34506929034156810e-04, + 1.98445679338826757e-04, + -7.08163358203131886e-05, + 2.46655529768908249e-05, + -9.35890030086883823e-06}; + + const double x2 = x * x; + + // Independent linear terms. + const double p01 = std::fma(p[1], x, p[0]); + const double p23 = std::fma(p[3], x, p[2]); + const double p45 = std::fma(p[5], x, p[4]); + const double p67 = std::fma(p[7], x, p[6]); + const double p89 = std::fma(p[9], x, p[8]); + const double p1011 = std::fma(p[11], x, p[10]); + const double p1213 = std::fma(p[13], x, p[12]); + const double p1415 = std::fma(p[15], x, p[14]); + const double p1617 = std::fma(p[17], x, p[16]); + + // Combine using Horner's method in x². + double result = std::fma(3.05977060678449757e-06, x2, p1617); + result = std::fma(result, x2, p1415); + result = std::fma(result, x2, p1213); + result = std::fma(result, x2, p1011); + result = std::fma(result, x2, p89); + result = std::fma(result, x2, p67); + result = std::fma(result, x2, p45); + result = std::fma(result, x2, p23); + return std::fma(result, x2, p01); + } /** \ingroup opencl_kernels * diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index f3f01e50169..71142f70afc 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -24,12 +24,14 @@ namespace internal { */ template inline T erfcx_tail_correction(const T& u) { - static constexpr std::array p{0.000658749161529837803157, 0.0160837851487422766278, - 0.125781726111229246204, 0.360344899949804439429, - 0.305326634961232344035}; - static constexpr std::array q{-0.00233520497626869185443, -0.0605183413124413191178, - -0.527905102951428412248, -1.87295284992346047209, - -2.56852019228982242072}; + static constexpr std::array p{ + 0.000658749161529837803157, 0.0160837851487422766278, + 0.125781726111229246204, 0.360344899949804439429, + 0.305326634961232344035}; + static constexpr std::array q{ + -0.00233520497626869185443, -0.0605183413124413191178, + -0.527905102951428412248, -1.87295284992346047209, + -2.56852019228982242072}; T numerator = 0.0163153871373020978498; T denominator = -1.0; for (int i = 4; i >= 0; --i) { @@ -140,13 +142,13 @@ inline double erfcx_small(double x) { // independently. A single degree-18 chain is 18 dependent operations; two // chains of 9 halve that latency. static constexpr std::array even = {1.0, - 4.99999999999992839e-01, - 1.66666666667239644e-01, - 4.16666666458337179e-02, - 8.33333374332981443e-03, - 1.38888415444527033e-03, - 1.98445679338826757e-04, - 2.46655529768908249e-05}; + 4.99999999999992839e-01, + 1.66666666667239644e-01, + 4.16666666458337179e-02, + 8.33333374332981443e-03, + 1.38888415444527033e-03, + 1.98445679338826757e-04, + 2.46655529768908249e-05}; static constexpr std::array odd = {-1.12837916709551256, -7.52252778063651983e-01, -3.00901111227312890e-01, -8.59717459974174147e-02, @@ -159,7 +161,7 @@ inline double erfcx_small(double x) { e = even[i] + x2 * e; o = odd[i] + x2 * o; } - return 1.0+ x2 * e + x * o; + return 1.0 + x2 * e + x * o; } } // namespace internal From cd7227a9c0d92f40a83912fa3b407afd7c2d7b32 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 22:07:28 +0300 Subject: [PATCH 21/24] Fix the out-of-bounds erfcx_small loop and the C++ constructs in the OpenCL source --- .../opencl/kernels/device_functions/erfcx.hpp | 112 +++++++++--------- stan/math/prim/fun/erfcx.hpp | 8 +- 2 files changed, 63 insertions(+), 57 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index afe9b47d05a..5ad28e1a2e2 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -89,12 +89,12 @@ static constexpr const char* erfcx_device_function * @return scaled complementary error function */ inline double erfcx_cody_middle(double y) { - constexpr std::array p{1230.33935479799725, 2051.07837782607147, - 1712.04761263407058, 881.952221241769090, - 298.635138197400131, 66.1191906371416295, - 8.88314979438837594, 5.64188496988670089e-1}; + const double p[] = {1230.33935479799725, 2051.07837782607147, + 1712.04761263407058, 881.952221241769090, + 298.635138197400131, 66.1191906371416295, + 8.88314979438837594, 5.64188496988670089e-1}; - constexpr std::array q{ + const double q[] = { 1230.33935480374942, 3439.36767414372164, 4362.61909014324716, 3290.79923573345963, 1621.38957456669019, 537.181101862009858, 117.693950891312499, 15.7449261107098347, @@ -102,25 +102,25 @@ static constexpr const char* erfcx_device_function const double y2 = y * y; - const double p01 = std::fma(p[1], y, p[0]); - const double p23 = std::fma(p[3], y, p[2]); - const double p45 = std::fma(p[5], y, p[4]); - const double p67 = std::fma(p[7], y, p[6]); + const double p01 = fma(p[1], y, p[0]); + const double p23 = fma(p[3], y, p[2]); + const double p45 = fma(p[5], y, p[4]); + const double p67 = fma(p[7], y, p[6]); - const double q01 = std::fma(q[1], y, q[0]); - const double q23 = std::fma(q[3], y, q[2]); - const double q45 = std::fma(q[5], y, q[4]); - const double q67 = std::fma(q[7], y, q[6]); + const double q01 = fma(q[1], y, q[0]); + const double q23 = fma(q[3], y, q[2]); + const double q45 = fma(q[5], y, q[4]); + const double q67 = fma(q[7], y, q[6]); - double num = std::fma(2.15311535474403846e-8, y2, p67); - num = std::fma(num, y2, p45); - num = std::fma(num, y2, p23); - num = std::fma(num, y2, p01); + double num = fma(2.15311535474403846e-8, y2, p67); + num = fma(num, y2, p45); + num = fma(num, y2, p23); + num = fma(num, y2, p01); double den = y2 + q67; - den = std::fma(den, y2, q45); - den = std::fma(den, y2, q23); - den = std::fma(den, y2, q01); + den = fma(den, y2, q45); + den = fma(den, y2, q23); + den = fma(den, y2, q01); return num / den; } @@ -135,48 +135,48 @@ static constexpr const char* erfcx_device_function * @return scaled complementary error function */ inline double erfcx_small(double x) { - constexpr std::array p{1.0, - -1.12837916709551256, - 1.0, - -7.52252778063651983e-01, - 4.99999999999992839e-01, - -3.00901111227312890e-01, - 1.66666666667239644e-01, - -8.59717459974174147e-02, - 4.16666666458337179e-02, - -1.91048337772546720e-02, - 8.33333374332981443e-03, - -3.47359067853470795e-03, - 1.38888415444527033e-03, - -5.34506929034156810e-04, - 1.98445679338826757e-04, - -7.08163358203131886e-05, - 2.46655529768908249e-05, - -9.35890030086883823e-06}; + const double p[] = {1.0, + -1.12837916709551256, + 1.0, + -7.52252778063651983e-01, + 4.99999999999992839e-01, + -3.00901111227312890e-01, + 1.66666666667239644e-01, + -8.59717459974174147e-02, + 4.16666666458337179e-02, + -1.91048337772546720e-02, + 8.33333374332981443e-03, + -3.47359067853470795e-03, + 1.38888415444527033e-03, + -5.34506929034156810e-04, + 1.98445679338826757e-04, + -7.08163358203131886e-05, + 2.46655529768908249e-05, + -9.35890030086883823e-06}; const double x2 = x * x; // Independent linear terms. - const double p01 = std::fma(p[1], x, p[0]); - const double p23 = std::fma(p[3], x, p[2]); - const double p45 = std::fma(p[5], x, p[4]); - const double p67 = std::fma(p[7], x, p[6]); - const double p89 = std::fma(p[9], x, p[8]); - const double p1011 = std::fma(p[11], x, p[10]); - const double p1213 = std::fma(p[13], x, p[12]); - const double p1415 = std::fma(p[15], x, p[14]); - const double p1617 = std::fma(p[17], x, p[16]); + const double p01 = fma(p[1], x, p[0]); + const double p23 = fma(p[3], x, p[2]); + const double p45 = fma(p[5], x, p[4]); + const double p67 = fma(p[7], x, p[6]); + const double p89 = fma(p[9], x, p[8]); + const double p1011 = fma(p[11], x, p[10]); + const double p1213 = fma(p[13], x, p[12]); + const double p1415 = fma(p[15], x, p[14]); + const double p1617 = fma(p[17], x, p[16]); // Combine using Horner's method in x². - double result = std::fma(3.05977060678449757e-06, x2, p1617); - result = std::fma(result, x2, p1415); - result = std::fma(result, x2, p1213); - result = std::fma(result, x2, p1011); - result = std::fma(result, x2, p89); - result = std::fma(result, x2, p67); - result = std::fma(result, x2, p45); - result = std::fma(result, x2, p23); - return std::fma(result, x2, p01); + double result = fma(3.05977060678449757e-06, x2, p1617); + result = fma(result, x2, p1415); + result = fma(result, x2, p1213); + result = fma(result, x2, p1011); + result = fma(result, x2, p89); + result = fma(result, x2, p67); + result = fma(result, x2, p45); + result = fma(result, x2, p23); + return fma(result, x2, p01); } /** \ingroup opencl_kernels diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 71142f70afc..64a4ed0fa23 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -154,10 +154,16 @@ inline double erfcx_small(double x) { -3.00901111227312890e-01, -8.59717459974174147e-02, -1.91048337772546720e-02, -3.47359067853470795e-03, -5.34506929034156810e-04, -7.08163358203131886e-05}; + static_assert(even.size() == odd.size(), + "the two chains are stepped by one loop"); + // The leading coefficient of each chain seeds the accumulator, and the + // trailing 1.0 is added at the end, so both arrays hold the interior + // coefficients only. Derive the bound from the array rather than writing + // it out, so shortening an array cannot leave the loop reading past it. const double x2 = x * x; double e = 3.05977060678449757e-06; double o = -9.35890030086883823e-06; - for (int i = 8; i > 0; --i) { + for (int i = static_cast(even.size()) - 1; i >= 0; --i) { e = even[i] + x2 * e; o = odd[i] + x2 * o; } From 6aff1cd2f0d9d831a155f7591f49418044900190 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 22:15:46 +0300 Subject: [PATCH 22/24] Keep the previous OpenCL erfcx device function, which is 8 to 13 percent faster --- .../opencl/kernels/device_functions/erfcx.hpp | 125 ++++++------------ 1 file changed, 44 insertions(+), 81 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 5ad28e1a2e2..8c61159502a 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -88,41 +88,24 @@ static constexpr const char* erfcx_device_function * @param y argument * @return scaled complementary error function */ - inline double erfcx_cody_middle(double y) { - const double p[] = {1230.33935479799725, 2051.07837782607147, - 1712.04761263407058, 881.952221241769090, - 298.635138197400131, 66.1191906371416295, - 8.88314979438837594, 5.64188496988670089e-1}; - - const double q[] = { - 1230.33935480374942, 3439.36767414372164, 4362.61909014324716, - 3290.79923573345963, 1621.38957456669019, 537.181101862009858, - 117.693950891312499, 15.7449261107098347, - }; - - const double y2 = y * y; - - const double p01 = fma(p[1], y, p[0]); - const double p23 = fma(p[3], y, p[2]); - const double p45 = fma(p[5], y, p[4]); - const double p67 = fma(p[7], y, p[6]); - - const double q01 = fma(q[1], y, q[0]); - const double q23 = fma(q[3], y, q[2]); - const double q45 = fma(q[5], y, q[4]); - const double q67 = fma(q[7], y, q[6]); - - double num = fma(2.15311535474403846e-8, y2, p67); - num = fma(num, y2, p45); - num = fma(num, y2, p23); - num = fma(num, y2, p01); - - double den = y2 + q67; - den = fma(den, y2, q45); - den = fma(den, y2, q23); - den = fma(den, y2, q01); - - return num / den; + double erfcx_cody_middle(double y) { + double p = 2.15311535474403846e-8 * y; + p = (p + 5.64188496988670089e-1) * y; + p = (p + 8.88314979438837594) * y; + p = (p + 66.1191906371416295) * y; + p = (p + 298.635138197400131) * y; + p = (p + 881.952221241769090) * y; + p = (p + 1712.04761263407058) * y; + p = (p + 2051.07837782607147) * y; + double q = y; + q = (q + 15.7449261107098347) * y; + q = (q + 117.693950891312499) * y; + q = (q + 537.181101862009858) * y; + q = (q + 1621.38957456669019) * y; + q = (q + 3290.79923573345963) * y; + q = (q + 4362.61909014324716) * y; + q = (q + 3439.36767414372164) * y; + return (p + 1230.33935479799725) / (q + 1230.33935480374942); } /** \ingroup opencl_kernels @@ -134,49 +117,26 @@ static constexpr const char* erfcx_device_function * @param x argument * @return scaled complementary error function */ - inline double erfcx_small(double x) { - const double p[] = {1.0, - -1.12837916709551256, - 1.0, - -7.52252778063651983e-01, - 4.99999999999992839e-01, - -3.00901111227312890e-01, - 1.66666666667239644e-01, - -8.59717459974174147e-02, - 4.16666666458337179e-02, - -1.91048337772546720e-02, - 8.33333374332981443e-03, - -3.47359067853470795e-03, - 1.38888415444527033e-03, - -5.34506929034156810e-04, - 1.98445679338826757e-04, - -7.08163358203131886e-05, - 2.46655529768908249e-05, - -9.35890030086883823e-06}; - - const double x2 = x * x; - - // Independent linear terms. - const double p01 = fma(p[1], x, p[0]); - const double p23 = fma(p[3], x, p[2]); - const double p45 = fma(p[5], x, p[4]); - const double p67 = fma(p[7], x, p[6]); - const double p89 = fma(p[9], x, p[8]); - const double p1011 = fma(p[11], x, p[10]); - const double p1213 = fma(p[13], x, p[12]); - const double p1415 = fma(p[15], x, p[14]); - const double p1617 = fma(p[17], x, p[16]); - - // Combine using Horner's method in x². - double result = fma(3.05977060678449757e-06, x2, p1617); - result = fma(result, x2, p1415); - result = fma(result, x2, p1213); - result = fma(result, x2, p1011); - result = fma(result, x2, p89); - result = fma(result, x2, p67); - result = fma(result, x2, p45); - result = fma(result, x2, p23); - return fma(result, x2, p01); + double erfcx_small(double x) { + double p = 3.05977060678449757e-06; + p = -9.35890030086883823e-06 + x * p; + p = 2.46655529768908249e-05 + x * p; + p = -7.08163358203131886e-05 + x * p; + p = 1.98445679338826757e-04 + x * p; + p = -5.34506929034156810e-04 + x * p; + p = 1.38888415444527033e-03 + x * p; + p = -3.47359067853470795e-03 + x * p; + p = 8.33333374332981443e-03 + x * p; + p = -1.91048337772546720e-02 + x * p; + p = 4.16666666458337179e-02 + x * p; + p = -8.59717459974174147e-02 + x * p; + p = 1.66666666667239644e-01 + x * p; + p = -3.00901111227312890e-01 + x * p; + p = 4.99999999999992839e-01 + x * p; + p = -7.52252778063651983e-01 + x * p; + p = 1.0 + x * p; + p = -1.12837916709551256 + x * p; + return 1.0 + x * p; } /** \ingroup opencl_kernels @@ -205,11 +165,14 @@ static constexpr const char* erfcx_device_function double erfcx(double x) { if (x >= 4.0) { return erfcx_cody_tail(x); - } else if (x >= 0.46875) { + } + if (x >= 0.46875) { return erfcx_cody_middle(x); - } else if (x > -0.46875) { + } + if (x > -0.46875) { return erfcx_small(x); - } else if (x < -27.0) { + } + if (x < -27.0) { return INFINITY; } double h = x * x; From 17fcb9d363660eeb72fb527680146fd085105902 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 18 Sep 2026 23:49:37 +0300 Subject: [PATCH 23/24] Apply the new review suggestions; pair the Cody middle coefficients without libm fma --- .../opencl/kernels/device_functions/erfcx.hpp | 9 ++-- stan/math/prim/fun/erfcx.hpp | 45 ++++++++++++------- stan/math/rev/fun/erfcx.hpp | 2 +- test/unit/math/rev/fun/erfcx_test.cpp | 21 +++++---- 4 files changed, 44 insertions(+), 33 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 8c61159502a..7c771d17d12 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -165,14 +165,11 @@ static constexpr const char* erfcx_device_function double erfcx(double x) { if (x >= 4.0) { return erfcx_cody_tail(x); - } - if (x >= 0.46875) { + } else if (x >= 0.46875) { return erfcx_cody_middle(x); - } - if (x > -0.46875) { + } else if (x > -0.46875) { return erfcx_small(x); - } - if (x < -27.0) { + } else if (x < -27.0) { return INFINITY; } double h = x * x; diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index 64a4ed0fa23..ff6dbcd3c8b 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -97,23 +97,34 @@ inline T erfcx_derivative(const T& x, const T& value) { * @return scaled complementary error function */ inline double erfcx_cody_middle(double y) { - double p = 2.15311535474403846e-8 * y; - p = (p + 5.64188496988670089e-1) * y; - p = (p + 8.88314979438837594) * y; - p = (p + 66.1191906371416295) * y; - p = (p + 298.635138197400131) * y; - p = (p + 881.952221241769090) * y; - p = (p + 1712.04761263407058) * y; - p = (p + 2051.07837782607147) * y; - double q = y; - q = (q + 15.7449261107098347) * y; - q = (q + 117.693950891312499) * y; - q = (q + 537.181101862009858) * y; - q = (q + 1621.38957456669019) * y; - q = (q + 3290.79923573345963) * y; - q = (q + 4362.61909014324716) * y; - q = (q + 3439.36767414372164) * y; - return (p + 1230.33935479799725) / (q + 1230.33935480374942); + static constexpr std::array p{1230.33935479799725, 2051.07837782607147, + 1712.04761263407058, 881.952221241769090, + 298.635138197400131, 66.1191906371416295, + 8.88314979438837594, 5.64188496988670089e-1}; + static constexpr std::array q{1230.33935480374942, 3439.36767414372164, + 4362.61909014324716, 3290.79923573345963, + 1621.38957456669019, 537.181101862009858, + 117.693950891312499, 15.7449261107098347}; + // Pair the coefficients first, so the four products are independent, then + // run a Horner chain half as long in y2. Written as arithmetic, not + // std::fma: with the Stan Math flags FP_FAST_FMA is not defined, so an + // explicit std::fma is a libm call. Measured on a Xeon E5-2680 v3, the + // std::fma form costs 31.0 ns per call against 4.46 ns for this one. The + // compiler contracts these into hardware fma wherever -march allows it. + const double y2 = y * y; + std::array p_vals; + std::array q_vals; + for (int i = 0, j = 0; i < 4; ++i, j += 2) { + p_vals[i] = p[j] + p[j + 1] * y; + q_vals[i] = q[j] + q[j + 1] * y; + } + double num = p_vals[3] + 2.15311535474403846e-8 * y2; + double den = y2 + q_vals[3]; + for (int i = 2; i >= 0; --i) { + num = p_vals[i] + num * y2; + den = q_vals[i] + den * y2; + } + return num / den; } /** diff --git a/stan/math/rev/fun/erfcx.hpp b/stan/math/rev/fun/erfcx.hpp index e911d84aa23..066f142987a 100644 --- a/stan/math/rev/fun/erfcx.hpp +++ b/stan/math/rev/fun/erfcx.hpp @@ -59,7 +59,7 @@ inline auto erfcx(T&& a) { return make_callback_var(val, [a, val](auto& vi) mutable { // apply_scalar_binary returns a lazy binaryExpr for Eigen inputs, so // this allocates nothing, and calls the functor directly for scalars. - const auto& deriv = apply_scalar_binary( + auto deriv = apply_scalar_binary( [](double x, double v) { return internal::erfcx_derivative(x, v); }, a.val(), val); as_array_or_scalar(a.adj()) diff --git a/test/unit/math/rev/fun/erfcx_test.cpp b/test/unit/math/rev/fun/erfcx_test.cpp index af69ef69a50..62e0bc4be60 100644 --- a/test/unit/math/rev/fun/erfcx_test.cpp +++ b/test/unit/math/rev/fun/erfcx_test.cpp @@ -1,8 +1,10 @@ #include #include #include +#include #include #include +#include namespace { // d/dx erfcx(x) = 2 * x * erfcx(x) - 2 / sqrt(pi) @@ -25,17 +27,18 @@ TEST_F(AgradRev, erfcx_value_and_gradient) { // Reference derivatives from 2 * x * erfcx(x) - 2 / sqrt(pi) evaluated in // long double. TEST_F(AgradRev, erfcx_gradient_reference) { - struct { - double x; - double d; - } cases[] = {{-2.0, -436.89199672700738}, {-1.0, -11.14633932862008}, - {0.0, -1.1283791670955126}, {1.0, -0.27321201478389856}, - {5.0, -0.021332789764826311}, {10.0, -0.0055593122190608565}}; - for (auto c : cases) { - stan::math::var x = c.x; + constexpr std::array, 6> cases{ + {{-2.0, -436.89199672700738}, + {-1.0, -11.14633932862008}, + {0.0, -1.1283791670955126}, + {1.0, -0.27321201478389856}, + {5.0, -0.021332789764826311}, + {10.0, -0.0055593122190608565}}}; + for (auto [x_val, d_ref] : cases) { + stan::math::var x = x_val; stan::math::var y = stan::math::erfcx(x); y.grad(); - EXPECT_NEAR(c.d, x.adj(), 1e-12 * std::fabs(c.d)); + EXPECT_NEAR(d_ref, x.adj(), 1e-12 * std::fabs(d_ref)); stan::math::recover_memory(); } } From 3612a0bda4ae76b4a9fb4413947f964d0d04877a Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 19 Sep 2026 15:44:20 +0300 Subject: [PATCH 24/24] Apply the 2026-09-19 review; keep the cheap derivative and share the tail polynomials --- .../opencl/kernels/device_functions/erfcx.hpp | 34 ++- stan/math/prim/fun/erfcx.hpp | 243 ++++++++++-------- 2 files changed, 153 insertions(+), 124 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/erfcx.hpp b/stan/math/opencl/kernels/device_functions/erfcx.hpp index 7c771d17d12..f7188b20f6f 100644 --- a/stan/math/opencl/kernels/device_functions/erfcx.hpp +++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp @@ -60,9 +60,11 @@ static constexpr const char* erfcx_device_function * * That difference cancels for large x: both terms approach * 2 / sqrt(pi) while the result decays like 1 / (sqrt(pi) * x^2). - * The difference form reaches 2.55e+11 ulp at x = 1e6. For - * x >= 4 the constant cancels analytically against the leading - * term of the tail rational, leaving 2 * u * C(u). + * For x >= 4 the constant cancels analytically against the + * leading term of the tail rational, leaving 2 * u * C(u). + * Above 30 it comes from the asymptotic expansion of + * -sqrt(pi) * x^2 * erfcx'(x) in u, whose coefficients satisfy + * c_{n+1} = -(n + 3/2) * c_n. * * Takes the value as an argument so the reverse pass does not * evaluate erfcx twice. @@ -72,11 +74,22 @@ static constexpr const char* erfcx_device_function * @return derivative of erfcx at x */ double erfcx_derivative(double x, double value) { - if (x >= 4.0) { - double u = 1.0 / (x * x); + if (x < 4.0) { + return 2.0 * x * value - M_2_SQRTPI; + } + double u = 1.0 / (x * x); + if (x < 30.0) { return 2.0 * u * erfcx_tail_correction(u); } - return 2.0 * x * value - M_2_SQRTPI; + // -sqrt(pi) * x^2 * erfcx'(x), asymptotic, ascending in u + const double s[8] + = {1.0, -1.5, 3.75, -13.125, + 59.0625, -324.84375, 2111.484375, -15836.1328125}; + double series = s[7]; + for (int i = 6; i >= 0; --i) { + series = series * u + s[i]; + } + return -(0.5 * M_2_SQRTPI) * u * series; } /** \ingroup opencl_kernels @@ -146,11 +159,10 @@ static constexpr const char* erfcx_device_function * * Mirrors stan/math/prim/fun/erfcx.hpp branch for branch and * formula for formula; the two must be kept in step. The whole - * positive axis is covered without a library call, which is - * what makes this fast. On the negative side exp is - * unavoidable, since erfcx grows like 2*exp(x*x); its argument - * is corrected with fma, because exp amplifies the rounding of - * x * x into roughly x * x * eps, which is 512 ulp at x = -26. + * positive axis is covered without a library call. On the + * negative side exp is unavoidable, since erfcx grows like + * 2*exp(x*x); its argument is corrected with fma, because exp + * amplifies the rounding of x * x into roughly x * x * eps. * * The correction must be written with fma and not with a Dekker * split. The split relies on t - (t - x) not being simplified diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp index ff6dbcd3c8b..6d01dd8d7e4 100644 --- a/stan/math/prim/fun/erfcx.hpp +++ b/stan/math/prim/fun/erfcx.hpp @@ -4,56 +4,102 @@ #include #include #include +#include #include +#include +#include namespace stan { namespace math { namespace internal { +// Cody (1969) third interval, ascending powers of `u = 1 / x^2`. +inline constexpr std::array erfcx_tail_p{ + 0.000658749161529837803157, 0.0160837851487422766278, + 0.125781726111229246204, 0.360344899949804439429, + 0.305326634961232344035, 0.0163153871373020978498}; +inline constexpr std::array erfcx_tail_q{ + -0.00233520497626869185443, -0.0605183413124413191178, + -0.527905102951428412248, -1.87295284992346047209, + -2.56852019228982242072, -1.0}; + +// above this the correction is below eps / 8 and x^12 overflows +inline constexpr double erfcx_tail_leading_only = 0x1p27; + +// Cody (1969) second interval, ascending powers of `y`. +inline constexpr std::array erfcx_middle_p{ + 1230.33935479799725, 2051.07837782607147, 1712.04761263407058, + 881.952221241769090, 298.635138197400131, 66.1191906371416295, + 8.88314979438837594, 5.64188496988670089e-1, 2.15311535474403846e-8}; +inline constexpr std::array erfcx_middle_q{ + 1230.33935480374942, 3439.36767414372164, 4362.61909014324716, + 3290.79923573345963, 1621.38957456669019, 537.181101862009858, + 117.693950891312499, 15.7449261107098347, 1.0}; + /** - * Correction factor of the Cody (1969) third-interval rational: - * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. + * The two tail polynomials, taking the coefficients in the order the + * iterators give them as descending powers of `t`. * - * Split out from the value so that the derivative can reuse it. See - * `erfcx_derivative`. + * @tparam N number of coefficients in each array + * @tparam T scalar type + * @tparam It coefficient iterator + * @param p first numerator coefficient + * @param q first denominator coefficient + * @param t argument + * @return numerator and denominator + */ +template +inline std::pair erfcx_tail_polynomials(It p, It q, const T& t) { + T numerator = *p; + T denominator = *q; + for (std::size_t i = 1; i < N; ++i) { + numerator = numerator * t + *++p; + denominator = denominator * t + *++q; + } + return {numerator, denominator}; +} + +/** + * Horner in `y^2` over adjacent coefficient pairs, so the two halves of the + * chain are independent. * * @tparam T scalar type - * @param u inverse square of the argument, `0 <= u <= 1/16` - * @return `P(u) / Q(u)` + * @tparam N number of coefficients + * @param c coefficients, ascending + * @param y argument + * @param y2 `y * y` + * @return polynomial value */ -template -inline T erfcx_tail_correction(const T& u) { - static constexpr std::array p{ - 0.000658749161529837803157, 0.0160837851487422766278, - 0.125781726111229246204, 0.360344899949804439429, - 0.305326634961232344035}; - static constexpr std::array q{ - -0.00233520497626869185443, -0.0605183413124413191178, - -0.527905102951428412248, -1.87295284992346047209, - -2.56852019228982242072}; - T numerator = 0.0163153871373020978498; - T denominator = -1.0; - for (int i = 4; i >= 0; --i) { - numerator = numerator * u + p[i]; - denominator = denominator * u + q[i]; +template +inline T erfcx_paired_horner(const std::array& c, const T& y, + const T& y2) { + int j = static_cast(N) - (N % 2 ? 3 : 4); + T r = N % 2 ? T(c[N - 1]) : T(c[N - 2] + c[N - 1] * y); + for (; j >= 0; j -= 2) { + r = r * y2 + (c[j] + c[j + 1] * y); } - return numerator / denominator; + return r; } /** * Cody (1969) third-interval rational, valid for `x >= 4`. * - * Gives `erfcx` directly. `x * x` is infinite for `x` large enough, which - * correctly collapses this to the leading term `INV_SQRT_PI / x` and, at - * infinity, to zero. + * `erfcx(x) = (INV_SQRT_PI + u * C(u)) / x` with `u = 1 / x^2` and `C` the + * Cody correction. `u * C(u)` is formed from the coefficients in `x^2`, so + * `1 / x^2` is never formed. * * @param x argument, `x >= 4` * @return scaled complementary error function */ inline double erfcx_cody_tail(double x) { - const double u = 1.0 / (x * x); - return (INV_SQRT_PI + u * erfcx_tail_correction(u)) / x; + if (x > erfcx_tail_leading_only) { + return INV_SQRT_PI / x; + } + const double s = x * x; + const auto pq = erfcx_tail_polynomials( + erfcx_tail_p.cbegin(), erfcx_tail_q.cbegin(), s); + return (INV_SQRT_PI + pq.first / (s * pq.second)) / x; } /** @@ -61,16 +107,17 @@ inline double erfcx_cody_tail(double x) { * * That difference cancels for large `x`: both terms approach * `2 / sqrt(pi)` while the result decays like `1 / (sqrt(pi) * x^2)`. - * Measured against a 50-digit reference, the difference form gives 6.1 ulp - * at `x = 4` and 2.55e+11 ulp at `x = 1e6`. - * * For `x >= 4` the tail rational gives the derivative with no subtraction, * because the constant cancels analytically: * * `2 * x * (INV_SQRT_PI + u * C(u)) / x - 2 / sqrt(pi) = 2 * u * C(u)` * - * since `2 * INV_SQRT_PI` is `2 / sqrt(pi)`. That form measures 0.1 to - * 31 ulp over the same range. + * since `2 * INV_SQRT_PI` is `2 / sqrt(pi)`. The walk in `x^2` is the one + * `erfcx_cody_tail` uses. + * + * Above 30 the derivative is taken from the asymptotic expansion of + * `-sqrt(pi) * x^2 * erfcx'(x)` in `u = 1 / x^2`, whose coefficients + * satisfy `c_{n+1} = -(n + 3/2) * c_n`. * * @tparam T scalar type * @param x argument @@ -79,11 +126,21 @@ inline double erfcx_cody_tail(double x) { */ template inline T erfcx_derivative(const T& x, const T& value) { - if (x >= 4.0) { - const T u = 1.0 / (x * x); - return 2.0 * u * erfcx_tail_correction(u); + if (x < 4.0) { + return 2.0 * x * value - TWO_OVER_SQRT_PI; + } + const T x2 = x * x; + if (x < 30.0) { + const auto pq = erfcx_tail_polynomials( + erfcx_tail_p.cbegin(), erfcx_tail_q.cbegin(), x2); + return 2.0 * pq.first / (x2 * pq.second); } - return 2.0 * x * value - TWO_OVER_SQRT_PI; + const T u = 1.0 / x2; + static constexpr std::array series_coefficients{ + 1.0, -1.5, 3.75, -13.125, + 59.0625, -324.84375, 2111.484375, -15836.1328125}; + const T series = erfcx_paired_horner(series_coefficients, u, u * u); + return -INV_SQRT_PI * u * series; } /** @@ -91,40 +148,15 @@ inline T erfcx_derivative(const T& x, const T& value) { * * Cody writes this interval as `erfc(x) = exp(-x*x) * R(x)`, so `erfcx` * is `R(x)` with the exponential cancelled analytically: no `exp` and no - * `erfc` call. Measured at 3.7 to 7.0 ulp. + * `erfc` call. * * @param y argument, `0.46875 <= y <= 4` * @return scaled complementary error function */ inline double erfcx_cody_middle(double y) { - static constexpr std::array p{1230.33935479799725, 2051.07837782607147, - 1712.04761263407058, 881.952221241769090, - 298.635138197400131, 66.1191906371416295, - 8.88314979438837594, 5.64188496988670089e-1}; - static constexpr std::array q{1230.33935480374942, 3439.36767414372164, - 4362.61909014324716, 3290.79923573345963, - 1621.38957456669019, 537.181101862009858, - 117.693950891312499, 15.7449261107098347}; - // Pair the coefficients first, so the four products are independent, then - // run a Horner chain half as long in y2. Written as arithmetic, not - // std::fma: with the Stan Math flags FP_FAST_FMA is not defined, so an - // explicit std::fma is a libm call. Measured on a Xeon E5-2680 v3, the - // std::fma form costs 31.0 ns per call against 4.46 ns for this one. The - // compiler contracts these into hardware fma wherever -march allows it. const double y2 = y * y; - std::array p_vals; - std::array q_vals; - for (int i = 0, j = 0; i < 4; ++i, j += 2) { - p_vals[i] = p[j] + p[j + 1] * y; - q_vals[i] = q[j] + q[j + 1] * y; - } - double num = p_vals[3] + 2.15311535474403846e-8 * y2; - double den = y2 + q_vals[3]; - for (int i = 2; i >= 0; --i) { - num = p_vals[i] + num * y2; - den = q_vals[i] + den * y2; - } - return num / den; + return erfcx_paired_horner(erfcx_middle_p, y, y2) + / erfcx_paired_horner(erfcx_middle_q, y, y2); } /** @@ -133,52 +165,45 @@ inline double erfcx_cody_middle(double y) { * A Chebyshev-economized expansion of `erfcx` about zero, so it needs no * library call and no branch on the sign of `x`. Its low-order coefficients * reproduce the Maclaurin series of `erfcx` exactly (`1`, `-2/sqrt(pi)`, - * `1`, ...), which is a useful check that the fit is right. `even[0]`, - * `even[1]` and `odd[0]` are those three terms. + * `1`, ...), which is a useful check that the fit is right. * - * Measured at 2.2 ulp over the interval. Evaluating the same coefficients - * as one degree-18 Horner chain gives 1.5 ulp instead, but runs 1.7 times - * slower, because that chain is 18 dependent operations. The whole function - * is 7.0 ulp, set by the middle interval, so the 0.7 ulp costs nothing. - * - * The plain Maclaurin series is only usable to about `|x| = 0.125`; four - * further economized terms extend it to 0.46875 at no measurable cost, - * which is what removes the last `exp` call from the positive axis. + * Evaluated by Estrin's method: four independent degree-3 blocks are + * combined in `x^4` and `x^8`. * * @param x argument, `|x| < 0.46875` * @return scaled complementary error function */ inline double erfcx_small(double x) { - // Split into the even and odd powers of x, so the two Horner chains run - // independently. A single degree-18 chain is 18 dependent operations; two - // chains of 9 halve that latency. - static constexpr std::array even = {1.0, - 4.99999999999992839e-01, - 1.66666666667239644e-01, - 4.16666666458337179e-02, - 8.33333374332981443e-03, - 1.38888415444527033e-03, - 1.98445679338826757e-04, - 2.46655529768908249e-05}; - static constexpr std::array odd - = {-1.12837916709551256, -7.52252778063651983e-01, - -3.00901111227312890e-01, -8.59717459974174147e-02, - -1.91048337772546720e-02, -3.47359067853470795e-03, - -5.34506929034156810e-04, -7.08163358203131886e-05}; - static_assert(even.size() == odd.size(), - "the two chains are stepped by one loop"); - // The leading coefficient of each chain seeds the accumulator, and the - // trailing 1.0 is added at the end, so both arrays hold the interior - // coefficients only. Derive the bound from the array rather than writing - // it out, so shortening an array cannot leave the loop reading past it. + static constexpr std::array c{1.0, + -1.12837916709551256, + 1.0, + -7.52252778063651983e-01, + 4.99999999999992839e-01, + -3.00901111227312890e-01, + 1.66666666667239644e-01, + -8.59717459974174147e-02, + 4.16666666458337179e-02, + -1.91048337772546720e-02, + 8.33333374332981443e-03, + -3.47359067853470795e-03, + 1.38888415444527033e-03, + -5.34506929034156810e-04, + 1.98445679338826757e-04, + -7.08163358203131886e-05, + 2.46655529768908249e-05, + -9.35890030086883823e-06, + 3.05977060678449757e-06}; const double x2 = x * x; - double e = 3.05977060678449757e-06; - double o = -9.35890030086883823e-06; - for (int i = static_cast(even.size()) - 1; i >= 0; --i) { - e = even[i] + x2 * e; - o = odd[i] + x2 * o; + const double x4 = x2 * x2; + const double x8 = x4 * x4; + std::array quad; + for (int i = 0; i < 4; ++i) { + const int j = 4 * i + 2; + quad[i] = (c[j] + c[j + 1] * x) + (c[j + 2] + c[j + 3] * x) * x2; } - return 1.0 + x2 * e + x * o; + const double rest + = (quad[0] + quad[1] * x4) + (quad[2] + quad[3] * x4 + c[18] * x8) * x8; + return (c[0] + c[1] * x) + x2 * rest; } } // namespace internal @@ -200,9 +225,7 @@ inline double erfcx_small(double x) { * without cancellation. * * Four branches. The whole positive axis is covered without calling any - * library function, which is what makes this fast: a direct approximation - * of `erfcx` avoids `erfc`, and `erfc` dominates the cost of forming - * `exp(x*x) * erfc(x)`. + * library function: a direct approximation of `erfcx` needs no `erfc`. * * - `x >= 4`: Cody (1969) third-interval rational. * - `0.46875 <= x < 4`: Cody (1969) second-interval rational. Cody writes @@ -215,9 +238,8 @@ inline double erfcx_small(double x) { * one place a library call is unavoidable, since `erfcx` grows like * `2*exp(x*x)` as `x -> -infinity`. The rounding error of `x * x` is * recovered with `fma` and folded back in, because `exp` amplifies it - * into roughly `x * x * eps` -- 512 ulp at `x = -26` if left alone. - * Below `x = -6.1`, `erfcx(-x)` is under `eps/2` of `2*exp(x*x)`, so - * the subtraction is skipped. + * into roughly `x * x * eps`. Below `x = -6.1`, `erfcx(-x)` is under + * `eps/2` of `2*exp(x*x)`, so the subtraction is skipped. * * W. J. Cody, Math. Comp. 23(107):631-638 (1969). * https://doi.org/10.1090/S0025-5718-1969-0247736-4 @@ -225,11 +247,6 @@ inline double erfcx_small(double x) { * The crossover at 4 is the same one R's `pnorm` uses: it switches to the * Cody tail form at `|y| > M_SQRT_32`, and since the argument here is the * `erfc` argument `y / sqrt(2)`, that is `sqrt(32) / sqrt(2) = 4` exactly. - * Evaluated at 60 significant digits, this coefficient set holds to 8.6e-17 - * relative down to 4 and degrades past about 3.5, so 4 sits just inside its - * range. - * - * Measured worst case over the whole range is 7.0 ulp. * * The derivative is * \f$\frac{d}{dx}\mbox{erfcx}(x) = 2x\,\mbox{erfcx}(x) -