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 <stan/math/fwd/fun/digamma.hpp>
 #include <stan/math/fwd/fun/erf.hpp>
 #include <stan/math/fwd/fun/erfc.hpp>
+#include <stan/math/fwd/fun/erfcx.hpp>
 #include <stan/math/fwd/fun/exp.hpp>
 #include <stan/math/fwd/fun/exp2.hpp>
 #include <stan/math/fwd/fun/expm1.hpp>
diff --git a/stan/math/fwd/fun/erfcx.hpp b/stan/math/fwd/fun/erfcx.hpp
new file mode 100644
index 00000000000..06e8ae5ea56
--- /dev/null
+++ b/stan/math/fwd/fun/erfcx.hpp
@@ -0,0 +1,33 @@
+#ifndef STAN_MATH_FWD_FUN_ERFCX_HPP
+#define STAN_MATH_FWD_FUN_ERFCX_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/erfcx.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the scaled complementary error function of the argument.
+ *
+ * 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
+ * @return scaled complementary error function of the argument
+ */
+template <typename T>
+inline fvar<T> erfcx(const fvar<T>& x) {
+  T v = erfcx(x.val_);
+  return fvar<T>(v, x.d_ * internal::erfcx_derivative(x.val_, v));
+}
+
+}  // 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..df349a7b483 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 <stan/math/opencl/kernels/device_functions/binomial_coefficient_log.hpp>
 #include <stan/math/opencl/kernels/device_functions/beta.hpp>
 #include <stan/math/opencl/kernels/device_functions/digamma.hpp>
+#include <stan/math/opencl/kernels/device_functions/erfcx.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_logit.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_Phi.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_square.hpp>
@@ -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,
@@ -351,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
new file mode 100644
index 00000000000..f7188b20f6f
--- /dev/null
+++ b/stan/math/opencl/kernels/device_functions/erfcx.hpp
@@ -0,0 +1,205 @@
+#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX_HPP
+#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_ERFCX_HPP
+#ifdef STAN_OPENCL
+
+#include <stan/math/opencl/stringify.hpp>
+#include <string>
+
+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
+           *
+           * Correction factor of the Cody (1969) third-interval rational:
+           * erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x, u = 1 / x^2.
+           *
+           * 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_tail_correction(double u) {
+            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 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).
+           * 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.
+           *
+           * @param x argument
+           * @param value erfcx(x)
+           * @return derivative of erfcx at x
+           */
+          double erfcx_derivative(double x, double value) {
+            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);
+            }
+            // -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
+           *
+           * 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 branch for branch and
+           * formula for formula; the two must be kept in step. The whole
+           * 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
+           * 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) {
+              return erfcx_cody_tail(x);
+            } else if (x >= 0.46875) {
+              return erfcx_cody_middle(x);
+            } else if (x > -0.46875) {
+              return erfcx_small(x);
+            } else if (x < -27.0) {
+              return INFINITY;
+            }
+            double h = x * 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
+// \endcond
+
+}  // namespace opencl_kernels
+}  // namespace math
+}  // namespace stan
+
+#endif
+#endif
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 <stan/math/opencl/rev/elt_multiply.hpp>
 #include <stan/math/opencl/rev/erf.hpp>
 #include <stan/math/opencl/rev/erfc.hpp>
+#include <stan/math/opencl/rev/erfcx.hpp>
 #include <stan/math/opencl/rev/exp.hpp>
 #include <stan/math/opencl/rev/exp2.hpp>
 #include <stan/math/opencl/rev/expm1.hpp>
diff --git a/stan/math/opencl/rev/erfcx.hpp b/stan/math/opencl/rev/erfcx.hpp
new file mode 100644
index 00000000000..e5fdfeb097a
--- /dev/null
+++ b/stan/math/opencl/rev/erfcx.hpp
@@ -0,0 +1,40 @@
+#ifndef STAN_MATH_OPENCL_REV_ERFCX_HPP
+#define STAN_MATH_OPENCL_REV_ERFCX_HPP
+#ifdef STAN_OPENCL
+
+#include <stan/math/opencl/kernel_generator.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/value_of.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the elementwise `erfcx()` of a var_value<matrix_cl<double>>.
+ *
+ * The derivative `2 * x * erfcx(x) - 2 / sqrt(pi)` reuses the function
+ * 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.
+ */
+template <typename T,
+          require_all_kernel_expressions_and_none_scalar_t<T>* = nullptr>
+inline var_value<matrix_cl<double>> erfcx(const var_value<T>& A) {
+  return make_callback_var(
+      erfcx(A.val()), [A](vari_value<matrix_cl<double>>& res) mutable {
+        A.adj()
+            += elt_multiply(res.adj(), erfcx_derivative(A.val(), res.val()));
+      });
+}
+
+}  // 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 <stan/math/prim/fun/elt_multiply.hpp>
 #include <stan/math/prim/fun/erf.hpp>
 #include <stan/math/prim/fun/erfc.hpp>
+#include <stan/math/prim/fun/erfcx.hpp>
 #include <stan/math/prim/fun/eval.hpp>
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/exp2.hpp>
diff --git a/stan/math/prim/fun/erfcx.hpp b/stan/math/prim/fun/erfcx.hpp
new file mode 100644
index 00000000000..6d01dd8d7e4
--- /dev/null
+++ b/stan/math/prim/fun/erfcx.hpp
@@ -0,0 +1,327 @@
+#ifndef STAN_MATH_PRIM_FUN_ERFCX_HPP
+#define STAN_MATH_PRIM_FUN_ERFCX_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/functor/apply_scalar_unary.hpp>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <utility>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+
+// Cody (1969) third interval, ascending powers of `u = 1 / x^2`.
+inline constexpr std::array<double, 6> erfcx_tail_p{
+    0.000658749161529837803157, 0.0160837851487422766278,
+    0.125781726111229246204,    0.360344899949804439429,
+    0.305326634961232344035,    0.0163153871373020978498};
+inline constexpr std::array<double, 6> 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<double, 9> 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<double, 9> erfcx_middle_q{
+    1230.33935480374942, 3439.36767414372164, 4362.61909014324716,
+    3290.79923573345963, 1621.38957456669019, 537.181101862009858,
+    117.693950891312499, 15.7449261107098347, 1.0};
+
+/**
+ * The two tail polynomials, taking the coefficients in the order the
+ * iterators give them as descending powers of `t`.
+ *
+ * @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 <std::size_t N, typename T, typename It>
+inline std::pair<T, T> 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
+ * @tparam N number of coefficients
+ * @param c coefficients, ascending
+ * @param y argument
+ * @param y2 `y * y`
+ * @return polynomial value
+ */
+template <typename T, std::size_t N>
+inline T erfcx_paired_horner(const std::array<double, N>& c, const T& y,
+                             const T& y2) {
+  int j = static_cast<int>(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 r;
+}
+
+/**
+ * Cody (1969) third-interval rational, valid for `x >= 4`.
+ *
+ * `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) {
+  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.size()>(
+      erfcx_tail_p.cbegin(), erfcx_tail_q.cbegin(), s);
+  return (INV_SQRT_PI + pq.first / (s * pq.second)) / 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)`.
+ * 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)`. 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
+ * @param value `erfcx(x)`
+ * @return derivative of `erfcx` at `x`
+ */
+template <typename T>
+inline T erfcx_derivative(const T& x, const T& value) {
+  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.size()>(
+        erfcx_tail_p.cbegin(), erfcx_tail_q.cbegin(), x2);
+    return 2.0 * pq.first / (x2 * pq.second);
+  }
+  const T u = 1.0 / x2;
+  static constexpr std::array<double, 8> 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;
+}
+
+/**
+ * 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.
+ *
+ * @param y argument, `0.46875 <= y <= 4`
+ * @return scaled complementary error function
+ */
+inline double erfcx_cody_middle(double y) {
+  const double y2 = y * y;
+  return erfcx_paired_horner(erfcx_middle_p, y, y2)
+         / erfcx_paired_horner(erfcx_middle_q, y, y2);
+}
+
+/**
+ * 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.
+ *
+ * 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) {
+  static constexpr std::array<double, 19> 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;
+  const double x4 = x2 * x2;
+  const double x8 = x4 * x4;
+  std::array<double, 4> 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;
+  }
+  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
+
+/**
+ * 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.
+ *
+ * Four branches. The whole positive axis is covered without calling any
+ * 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
+ *   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.
+ * - `|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`. 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
+ *
+ * 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.
+ *
+ * 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 <typename T, require_arithmetic_t<T>* = nullptr>
+inline double erfcx(T&& xx) {
+  const double x = static_cast<double>(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);
+  } else if (x >= middle_min) {
+    return internal::erfcx_cody_middle(x);
+  } else if (x > -middle_min) {
+    return internal::erfcx_small(x);
+  } else if (x < overflow_max) {
+    return INFTY;
+  }
+  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;
+  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));
+}
+
+/**
+ * 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 <typename T>
+  static inline auto fun(T&& x) {
+    return erfcx(std::forward<T>(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<T>* = nullptr,
+    require_container_t<T>* = nullptr, require_not_var_matrix_t<T>* = nullptr>
+inline auto erfcx(T&& x) {
+  return apply_scalar_unary<erfcx_fun, T>::apply(std::forward<T>(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 <stan/math/rev/fun/elt_multiply.hpp>
 #include <stan/math/rev/fun/erf.hpp>
 #include <stan/math/rev/fun/erfc.hpp>
+#include <stan/math/rev/fun/erfcx.hpp>
 #include <stan/math/rev/fun/exp.hpp>
 #include <stan/math/rev/fun/exp2.hpp>
 #include <stan/math/rev/fun/expm1.hpp>
diff --git a/stan/math/rev/fun/erfcx.hpp b/stan/math/rev/fun/erfcx.hpp
new file mode 100644
index 00000000000..066f142987a
--- /dev/null
+++ b/stan/math/rev/fun/erfcx.hpp
@@ -0,0 +1,72 @@
+#ifndef STAN_MATH_REV_FUN_ERFCX_HPP
+#define STAN_MATH_REV_FUN_ERFCX_HPP
+
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/erfcx.hpp>
+#include <stan/math/prim/functor/apply_scalar_binary.hpp>
+#include <cmath>
+
+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. 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) =
+   \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]
+ *
+ * One overload covers `var` and `var_value<Matrix>`. For a matrix the
+ * derivative stays a lazy `binaryExpr`, so the reverse pass allocates
+ * nothing.
+ *
+ * @tparam T a `var` or a `var_value` of a matrix type
+ * @param a The variable.
+ * @return Scaled complementary error function applied to the variable,
+ *   elementwise for a matrix.
+ */
+template <
+    typename T, require_var_t<T>* = nullptr,
+    require_all_not_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
+inline auto erfcx(T&& a) {
+  auto val = to_arena(erfcx(a.val()));
+  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.
+    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);
+  });
+}
+
+}  // 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 <stan/math/fwd.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+
+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<double> x(x_value, 1.0);
+    fvar<double> 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<double> x(1.5, 3.0);
+  fvar<double> 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<fvar<double> > x;
+  x.val_.val_ = x_value;
+  x.val_.d_ = 1.0;
+  x.d_.val_ = 1.0;
+
+  fvar<fvar<double> > 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_derivative_test.cpp b/test/unit/math/mix/fun/erfcx_derivative_test.cpp
new file mode 100644
index 00000000000..c117bdb6ca7
--- /dev/null
+++ b/test/unit/math/mix/fun/erfcx_derivative_test.cpp
@@ -0,0 +1,100 @@
+#include <stan/math/mix.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+
+/**
+ * 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 <typename F>
+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<double> 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);
+}
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 <test/unit/math/test_ad.hpp>
+
+// 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<double> 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/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..d1cf93dd4ed
--- /dev/null
+++ b/test/unit/math/opencl/rev/erfcx_test.cpp
@@ -0,0 +1,78 @@
+#ifdef STAN_OPENCL
+#include <stan/math/opencl/rev.hpp>
+#include <stan/math.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/math/opencl/util.hpp>
+
+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);
+}
+
+/**
+ * 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<stan::math::matrix_cl<double>> 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;
+
+  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
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..9f9d272889e
--- /dev/null
+++ b/test/unit/math/prim/fun/erfcx_test.cpp
@@ -0,0 +1,133 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+// 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);
+}
+
+// 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;
+  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)));
+  }
+}
+
+// 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) {
+  using stan::math::erfcx;
+  const double inf = std::numeric_limits<double>::infinity();
+  const double nan = std::numeric_limits<double>::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<double> xs{-1.0, 0.0, 1.0, 20.0};
+  std::vector<double> 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..62e0bc4be60
--- /dev/null
+++ b/test/unit/math/rev/fun/erfcx_test.cpp
@@ -0,0 +1,86 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <gtest/gtest.h>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <utility>
+
+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) {
+  constexpr std::array<std::pair<double, double>, 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(d_ref, x.adj(), 1e-12 * std::fabs(d_ref));
+    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<stan::math::var, Eigen::Dynamic, 1> 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<Eigen::VectorXd> 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));
+  }
+}