diff --git a/stan/math/prim/fun/grad_reg_inc_beta.hpp b/stan/math/prim/fun/grad_reg_inc_beta.hpp
index 92fab1fc056..c34e463f968 100644
--- a/stan/math/prim/fun/grad_reg_inc_beta.hpp
+++ b/stan/math/prim/fun/grad_reg_inc_beta.hpp
@@ -420,7 +420,8 @@ inline auto inc_beta_continued_fraction_kernel(const T1& a, const T2& b,
  *   the pair (d/da, d/db) for <code>BOTH</code>
  */
 template <GRAD_RIB grads_requested, typename T1, typename T2, typename T3,
-          typename T4, typename T5, typename T6, typename T7>
+          typename T4, typename T5, typename T6, typename T7,
+          require_all_not_var_t<T1, T2, T3, T4, T5, T6, T7>* = nullptr>
 inline auto grad_reg_inc_beta(const T1& a, const T2& b, const T3& z,
                               const T4& digammaA, const T5& digammaB,
                               const T6& digammaSum, const T7& /* betaAB */) {
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index e5ef938a302..efb00634da1 100644
--- a/stan/math/rev/fun.hpp
+++ b/stan/math/rev/fun.hpp
@@ -71,6 +71,7 @@
 #include <stan/math/rev/fun/generalized_inverse.hpp>
 #include <stan/math/rev/fun/gp_periodic_cov.hpp>
 #include <stan/math/rev/fun/grad.hpp>
+#include <stan/math/rev/fun/grad_reg_inc_beta.hpp>
 #include <stan/math/rev/fun/hypergeometric_1F0.hpp>
 #include <stan/math/rev/fun/hypergeometric_2F1.hpp>
 #include <stan/math/rev/fun/hypergeometric_pFq.hpp>
diff --git a/stan/math/rev/fun/grad_reg_inc_beta.hpp b/stan/math/rev/fun/grad_reg_inc_beta.hpp
new file mode 100644
index 00000000000..daf151d46f0
--- /dev/null
+++ b/stan/math/rev/fun/grad_reg_inc_beta.hpp
@@ -0,0 +1,833 @@
+#ifndef STAN_MATH_REV_FUN_GRAD_REG_INC_BETA_HPP
+#define STAN_MATH_REV_FUN_GRAD_REG_INC_BETA_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/value_of.hpp>
+#include <stan/math/prim/fun/boost_policy.hpp>
+#include <stan/math/prim/fun/digamma.hpp>
+#include <stan/math/prim/fun/grad_reg_inc_beta.hpp>
+#include <stan/math/prim/fun/inc_beta_ddz.hpp>
+#include <stan/math/prim/fun/lbeta.hpp>
+#include <stan/math/prim/fun/log1p.hpp>
+#include <boost/math/special_functions/trigamma.hpp>
+#include <cmath>
+#include <limits>
+#include <utility>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+
+/**
+ * Hessian-vector product H v of the regularized incomplete beta function
+ * I_x(a, b), where H is the Hessian of I_x(a, b) in (a, b), computed from its
+ * power series.
+ *
+ * As in internal::inc_beta_power_series_kernel, the series is
+ *
+ *   I_x(a, b) = N S,  N = x^a / B(a, b),  S = sum_{n >= 0} t_n,
+ *   t_n = P_n / (a + n),  P_0 = 1,  P_n = P_{n-1} (n - b) x / n.
+ *
+ * Write grad for the gradient in (a, b) and D_v for the derivative along v.
+ * Differentiating I = N S twice and applying the result to v gives
+ *
+ *   H v = N (D_v grad S + grad S D_v log N + grad log N D_v S
+ *            + S (grad log N D_v log N + D_v grad log N)).
+ *
+ * When only one shape j is requested, the product needs only the diagonal
+ * entry
+ *
+ *   H_jj = N (d2S + 2 dlog_N dS + S (dlog_N^2 + d2log_N)),
+ *
+ * where d and d2 are the first and second derivatives in that shape.
+ *
+ * The derivatives of S are the sums of those of its terms. The a derivatives
+ * of a term are closed form. Its b derivatives come from those of P_n, which
+ * follow from the recurrence for P_n:
+ *
+ *   dt/da = -t / (a + n),           d2t/da2 = 2 t / (a + n)^2,
+ *   dt/db = P'_n / (a + n),         d2t/dadb = -(dt/db) / (a + n),
+ *   d2t/db2 = P''_n / (a + n).
+ *
+ * For n = 0 the term is 1 / a, with a derivatives -1 / a^2 and 2 / a^3, and
+ * it does not depend on b.
+ *
+ * With psi the digamma function and psi1 the trigamma function, the
+ * derivatives of log N = a log x - log B(a, b) are
+ *
+ *   grad log N = (log x - psi(a) + psi(a + b), psi(a + b) - psi(b)),
+ *   Hess log N = [[psi1(a + b) - psi1(a), psi1(a + b)],
+ *                 [psi1(a + b), psi1(a + b) - psi1(b)]].
+ *
+ * psi1 is taken from Boost. stan::math::trigamma(double) is accurate only to
+ * about 1e-9 for arguments below 5, and psi1(a + b) appears alone in H_ab,
+ * where that error does not cancel.
+ *
+ * The loop stops once the terms are below eps / a in magnitude, or after 600
+ * terms. For integer b, t_n and its a derivatives are 0 for n >= b, but its b
+ * derivatives are not. <code>G2</code> and <code>BOTH</code> therefore also
+ * wait for the b derivatives of the terms to be small.
+ *
+ * The variable names follow two conventions. For <code>G1</code> and
+ * <code>G2</code> only the requested shape varies, and dX and d2X are the
+ * first and second derivatives of X in that shape. For <code>BOTH</code>,
+ * X_a and X_b are the partials of X, and Dv_X is the derivative of X along v.
+ *
+ * @tparam grads <code>G1</code> for H_aa v_a, <code>G2</code> for H_bb v_b,
+ *   <code>BOTH</code> for H v
+ * @param a first shape
+ * @param b second shape
+ * @param x upper limit, in (0, 1)
+ * @param v_a a component of v, not used for <code>G2</code>
+ * @param v_b b component of v, not used for <code>G1</code>
+ * @param digamma_a digamma(a), not used for <code>G2</code>
+ * @param digamma_b digamma(b), not used for <code>G1</code>
+ * @param digamma_ab digamma(a + b)
+ * @return H_aa v_a for <code>G1</code>, H_bb v_b for <code>G2</code>,
+ *   ((H v)_a, (H v)_b) for <code>BOTH</code>
+ */
+template <GRAD_RIB grads>
+inline auto inc_beta_hvp_power_series(double a, double b, double x, double v_a,
+                                      double v_b, double digamma_a,
+                                      double digamma_b, double digamma_ab) {
+  constexpr int max_iterations = 600;
+  const double tolerance = std::numeric_limits<double>::epsilon() / a;
+  const double log_x = std::log(x);
+  const double N = std::exp(a * log_x - lbeta(a, b));
+  const double trigamma_ab = boost::math::trigamma(a + b, boost_policy_t<>());
+  // n = 0: P_0 = 1 and the term is 1 / a
+  double product = 1.0;
+  double S = 1.0 / a;
+  if constexpr (grads == GRAD_RIB::G1) {
+    double dS = -S / a;
+    double d2S = 2.0 * S / (a * a);
+    for (int n = 1; n <= max_iterations; ++n) {
+      product *= (n - b) * x / n;
+      const double inv_apn = 1.0 / (a + n);
+      const double term = product * inv_apn;
+      S += term;
+      dS -= term * inv_apn;
+      d2S += 2.0 * term * inv_apn * inv_apn;
+      // for integer b the term and its a derivatives are 0 from n = b on
+      if (std::fabs(term) <= tolerance) {
+        break;
+      }
+    }
+    const double dlog_N = log_x - digamma_a + digamma_ab;
+    const double d2log_N
+        = trigamma_ab - boost::math::trigamma(a, boost_policy_t<>());
+    return v_a * N
+           * (d2S + 2.0 * dlog_N * dS + S * (dlog_N * dlog_N + d2log_N));
+  } else if constexpr (grads == GRAD_RIB::G2) {
+    double dproduct = 0.0;
+    double d2product = 0.0;
+    double dS = 0.0;
+    double d2S = 0.0;
+    for (int n = 1; n <= max_iterations; ++n) {
+      // P_n = f P_{n-1} with f = (n - b) x / n and df/db = -x / n
+      const double x_div_n = x / n;
+      const double factor = (n - b) * x_div_n;
+      d2product = factor * d2product - 2.0 * x_div_n * dproduct;
+      dproduct = factor * dproduct - x_div_n * product;
+      product *= factor;
+      const double inv_apn = 1.0 / (a + n);
+      const double term = product * inv_apn;
+      const double dterm = dproduct * inv_apn;
+      const double d2term = d2product * inv_apn;
+      S += term;
+      dS += dterm;
+      d2S += d2term;
+      if (std::fabs(term) <= tolerance && std::fabs(dterm) <= tolerance
+          && std::fabs(d2term) <= tolerance) {
+        break;
+      }
+    }
+    const double dlog_N = digamma_ab - digamma_b;
+    const double d2log_N
+        = trigamma_ab - boost::math::trigamma(b, boost_policy_t<>());
+    return v_b * N
+           * (d2S + 2.0 * dlog_N * dS + S * (dlog_N * dlog_N + d2log_N));
+  } else {
+    double product_b = 0.0;
+    double product_bb = 0.0;
+    double S_a = -S / a;
+    double S_b = 0.0;
+    double Dv_S_a = 2.0 * v_a * S / (a * a);
+    double Dv_S_b = 0.0;
+    for (int n = 1; n <= max_iterations; ++n) {
+      // P_n = f P_{n-1} with f = (n - b) x / n and df/db = -x / n
+      const double x_div_n = x / n;
+      const double factor = (n - b) * x_div_n;
+      product_bb = factor * product_bb - 2.0 * x_div_n * product_b;
+      product_b = factor * product_b - x_div_n * product;
+      product *= factor;
+      const double inv_apn = 1.0 / (a + n);
+      const double term = product * inv_apn;
+      const double term_b = product_b * inv_apn;
+      const double term_aa = 2.0 * term * inv_apn * inv_apn;
+      const double term_ab = -term_b * inv_apn;
+      const double term_bb = product_bb * inv_apn;
+      S += term;
+      S_a -= term * inv_apn;
+      S_b += term_b;
+      Dv_S_a += v_a * term_aa + v_b * term_ab;
+      Dv_S_b += v_a * term_ab + v_b * term_bb;
+      if (std::fabs(term) <= tolerance && std::fabs(term_b) <= tolerance
+          && std::fabs(term_bb) <= tolerance) {
+        break;
+      }
+    }
+    const double log_N_a = log_x - digamma_a + digamma_ab;
+    const double log_N_b = digamma_ab - digamma_b;
+    const double Dv_log_N = log_N_a * v_a + log_N_b * v_b;
+    const double Dv_S = S_a * v_a + S_b * v_b;
+    const double Dv_log_N_a
+        = (trigamma_ab - boost::math::trigamma(a, boost_policy_t<>())) * v_a
+          + trigamma_ab * v_b;
+    const double Dv_log_N_b
+        = trigamma_ab * v_a
+          + (trigamma_ab - boost::math::trigamma(b, boost_policy_t<>())) * v_b;
+    return std::pair{N
+                         * (Dv_S_a + S_a * Dv_log_N + log_N_a * Dv_S
+                            + S * (log_N_a * Dv_log_N + Dv_log_N_a)),
+                     N
+                         * (Dv_S_b + S_b * Dv_log_N + log_N_b * Dv_S
+                            + S * (log_N_b * Dv_log_N + Dv_log_N_b))};
+  }
+}
+
+/**
+ * Hessian-vector product H v of the regularized incomplete beta function
+ * I_x(a, b), where H is the Hessian of I_x(a, b) in (a, b), computed from its
+ * continued fraction.
+ *
+ * As in internal::inc_beta_continued_fraction_kernel, I_x(a, b) = Q h with
+ * Q = x^a (1 - x)^b / (a B(a, b)) and h the continued fraction
+ * https://dlmf.nist.gov/8.17.E22, which is evaluated with the modified Lentz
+ * method. Write grad for the gradient in (a, b), Hess for the Hessian, and
+ * D_v for the derivative along v. Differentiating I = Q h twice gives
+ *
+ *   H = I (grad log I grad log I^T + Hess log Q + Hess log h),
+ *
+ * with grad log I = grad log Q + grad log h. With psi the digamma function
+ * and psi1 the trigamma function,
+ *
+ *   grad log Q = (log x - 1 / a - psi(a) + psi(a + b),
+ *                 log(1 - x) - psi(b) + psi(a + b)),
+ *   Hess log Q = [[psi1(a + b) - psi1(a) + 1 / a^2, psi1(a + b)],
+ *                 [psi1(a + b), psi1(a + b) - psi1(b)]],
+ *
+ * with psi1 taken from Boost, as in the power series.
+ *
+ * Each Lentz step multiplies h by C D, so log h is the sum of log C and
+ * log D over the steps, and its derivatives are the sums of theirs. With d
+ * the partial numerator of the step, the step sets C' = 1 + d / C and
+ * D' = 1 / (1 + d D). Their log derivatives are
+ *
+ *   grad log C' = (grad d - d grad log C) / (C + d),
+ *   grad log D' = -(grad d + d grad log D) D D',
+ *
+ * using C C' = C + d. Differentiating these once more along v gives
+ *
+ *   D_v grad log C' = (Hess d v - D_v d grad log C - d D_v grad log C)
+ *                     / (C + d) - grad log C' D_v log(C C'),
+ *   D_v log(C C') = (C D_v log C + D_v d) / (C + d),
+ *   D_v grad log D' = -(Hess d v + D_v d grad log D + d D_v grad log D) D D'
+ *                     + grad log D' D_v log(D D'),
+ *   D_v log(D D') = D_v log D - D D' (D_v d + d D_v log D).
+ *
+ * The code applies each of the two second-derivative formulas in three parts.
+ * It first computes the part that uses only values from before the step. It
+ * then overwrites grad log C (or grad log D) with its new value, and last it
+ * adds the final term, which uses that new value.
+ *
+ * <code>G1</code> and <code>G2</code> take v = (1, 0) and v = (0, 1), so D_v
+ * is the derivative in the requested shape, and the result is scaled by v_a
+ * or v_b at the end.
+ *
+ * The partial numerators https://dlmf.nist.gov/8.17.E23 are d_n = num / den,
+ * with den = (a + n - 1)(a + n) and
+ *
+ *   n = 2m:     num = m (b - m) x,           num_b = m x,
+ *   n = 2m + 1: num = -(a + m)(a + b + m) x, num_a = -(2 (a + m) + b) x,
+ *               num_b = -(a + m) x,          num_aa = -2 x, num_ab = -x.
+ *
+ * Every other derivative of num is 0, so d_bb = 0. The derivatives of d_n
+ * follow from the quotient rule, with den_a = 2 (a + n) - 1 and den_aa = 2.
+ * n = 1 is the odd case with m = 0. The first step holds C at 1, so
+ * D = 1 / (1 + d_1).
+ *
+ * The test for stopping runs after each odd n. The loop stops once C D is
+ * within 3 eps of 1 and the second derivatives of log h have changed by at
+ * most a relative 3 eps since the previous test, or after n = 601.
+ *
+ * The variable names follow the same two conventions as in the power series.
+ * For <code>G1</code> and <code>G2</code>, dX and d2X are the first and
+ * second derivatives of X in the requested shape. For <code>BOTH</code>,
+ * X_a and X_b are the partials of X, and Dv_X is the derivative of X along v.
+ *
+ * @tparam grads <code>G1</code> for H_aa v_a, <code>G2</code> for H_bb v_b,
+ *   <code>BOTH</code> for H v
+ * @param a first shape
+ * @param b second shape
+ * @param x upper limit, in (0, 1)
+ * @param v_a a component of v, not used for <code>G2</code>
+ * @param v_b b component of v, not used for <code>G1</code>
+ * @param digamma_a digamma(a), not used for <code>G2</code>
+ * @param digamma_b digamma(b), not used for <code>G1</code>
+ * @param digamma_ab digamma(a + b)
+ * @return H_aa v_a for <code>G1</code>, H_bb v_b for <code>G2</code>,
+ *   ((H v)_a, (H v)_b) for <code>BOTH</code>
+ */
+template <GRAD_RIB grads>
+inline auto inc_beta_hvp_continued_fraction(double a, double b, double x,
+                                            double v_a, double v_b,
+                                            double digamma_a, double digamma_b,
+                                            double digamma_ab) {
+  constexpr int max_iterations = 300;
+  const double tolerance = 3.0 * std::numeric_limits<double>::epsilon();
+  const double small = std::sqrt(std::numeric_limits<double>::min());
+  const double log_x = std::log(x);
+  const double log1m_x = log1p(-x);
+  const double Q
+      = std::exp(a * log_x + b * log1m_x - std::log(a) - lbeta(a, b));
+  // boost's trigamma, as in the power series
+  const double trigamma_ab = boost::math::trigamma(a + b, boost_policy_t<>());
+  // first step, n = 1, with C held at 1
+  const double ap1 = a + 1.0;
+  const double d_1 = -(a + b) * x / ap1;
+  double D = 1.0 + d_1;
+  if (std::fabs(D) < small) {
+    D = small;
+  }
+  D = 1.0 / D;
+  double C = 1.0;
+  double h = D;
+  if constexpr (grads == GRAD_RIB::G1) {
+    const double dd_1 = -(x + d_1) / ap1;
+    const double d2d_1 = 2.0 * (x + d_1) / (ap1 * ap1);
+    double dlog_C = 0.0;
+    double d2log_C = 0.0;
+    double dlog_D = -dd_1 * D;
+    double d2log_D = -(d2d_1 + dlog_D * dd_1) * D;
+    double dlog_h = dlog_D;
+    double d2log_h = d2log_D;
+    double d2log_h_prev = d2log_h;
+    for (int n = 2; n <= 2 * max_iterations + 1; ++n) {
+      const int m = n / 2;
+      const bool odd = n % 2 == 1;
+      const double a_plus_n = a + n;
+      const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0));
+      const double dden = 2.0 * a_plus_n - 1.0;
+      double num;
+      double dnum = 0.0;
+      double d2num = 0.0;
+      if (odd) {
+        num = -(a + m) * (a + b + m) * x;
+        dnum = -(2.0 * (a + m) + b) * x;
+        d2num = -2.0 * x;
+      } else {
+        num = m * (b - m) * x;
+      }
+      const double d_n = num * inv_den;
+      const double dd_n = (dnum - d_n * dden) * inv_den;
+      const double d2d_n = (d2num - 2.0 * dd_n * dden - 2.0 * d_n) * inv_den;
+
+      // one step of the modified Lentz method
+      // C_new = 1 + d_n / C, with C C_new = C + d_n
+      double C_new = 1.0 + d_n / C;
+      double C_C_new = C + d_n;
+      if (std::fabs(C_new) < small) {
+        C_new = small;
+        C_C_new = C * small;
+      }
+      const double inv_C_C_new = 1.0 / C_C_new;
+      const double dlog_C_C_new = (C * dlog_C + dd_n) * inv_C_C_new;
+      d2log_C = (d2d_n - dd_n * dlog_C - d_n * d2log_C) * inv_C_C_new;
+      dlog_C = (dd_n - d_n * dlog_C) * inv_C_C_new;
+      d2log_C -= dlog_C * dlog_C_C_new;
+
+      // D_new = 1 / (1 + d_n D)
+      double D_new = 1.0 + d_n * D;
+      if (std::fabs(D_new) < small) {
+        D_new = small;
+      }
+      D_new = 1.0 / D_new;
+      const double D_D_new = D * D_new;
+      const double dlog_D_D_new = dlog_D - D_D_new * (dd_n + d_n * dlog_D);
+      d2log_D = -(d2d_n + dd_n * dlog_D + d_n * d2log_D) * D_D_new;
+      dlog_D = -(dd_n + d_n * dlog_D) * D_D_new;
+      d2log_D += dlog_D * dlog_D_D_new;
+
+      dlog_h += dlog_C + dlog_D;
+      d2log_h += d2log_C + d2log_D;
+      C = C_new;
+      D = D_new;
+      const double delta = C_new * D_new;
+      h *= delta;
+      if (odd) {
+        if (std::fabs(delta - 1.0) < tolerance
+            && std::fabs(d2log_h - d2log_h_prev)
+                   <= tolerance * std::fabs(d2log_h)) {
+          break;
+        }
+        d2log_h_prev = d2log_h;
+      }
+    }
+    const double dlog_I = log_x - 1.0 / a + digamma_ab - digamma_a + dlog_h;
+    const double d2log_Q = trigamma_ab
+                           - boost::math::trigamma(a, boost_policy_t<>())
+                           + 1.0 / (a * a);
+    return v_a * h * Q * (dlog_I * dlog_I + d2log_Q + d2log_h);
+  } else if constexpr (grads == GRAD_RIB::G2) {
+    // d_n is linear in b, so d2d_n = 0
+    const double dd_1 = -x / ap1;
+    double dlog_C = 0.0;
+    double d2log_C = 0.0;
+    double dlog_D = -dd_1 * D;
+    double d2log_D = -dlog_D * dd_1 * D;
+    double dlog_h = dlog_D;
+    double d2log_h = d2log_D;
+    double d2log_h_prev = d2log_h;
+    for (int n = 2; n <= 2 * max_iterations + 1; ++n) {
+      const int m = n / 2;
+      const bool odd = n % 2 == 1;
+      const double a_plus_n = a + n;
+      const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0));
+      double dnum;
+      double num;
+      if (odd) {
+        dnum = -(a + m) * x;
+        num = dnum * (a + b + m);
+      } else {
+        dnum = m * x;
+        num = dnum * (b - m);
+      }
+      const double d_n = num * inv_den;
+      const double dd_n = dnum * inv_den;
+
+      // one step of the modified Lentz method, as for G1 with d2d_n = 0.
+      // C_new = 1 + d_n / C, with C C_new = C + d_n
+      double C_new = 1.0 + d_n / C;
+      double C_C_new = C + d_n;
+      if (std::fabs(C_new) < small) {
+        C_new = small;
+        C_C_new = C * small;
+      }
+      const double inv_C_C_new = 1.0 / C_C_new;
+      const double dlog_C_C_new = (C * dlog_C + dd_n) * inv_C_C_new;
+      d2log_C = -(dd_n * dlog_C + d_n * d2log_C) * inv_C_C_new;
+      dlog_C = (dd_n - d_n * dlog_C) * inv_C_C_new;
+      d2log_C -= dlog_C * dlog_C_C_new;
+
+      // D_new = 1 / (1 + d_n D)
+      double D_new = 1.0 + d_n * D;
+      if (std::fabs(D_new) < small) {
+        D_new = small;
+      }
+      D_new = 1.0 / D_new;
+      const double D_D_new = D * D_new;
+      const double dlog_D_D_new = dlog_D - D_D_new * (dd_n + d_n * dlog_D);
+      d2log_D = -(dd_n * dlog_D + d_n * d2log_D) * D_D_new;
+      dlog_D = -(dd_n + d_n * dlog_D) * D_D_new;
+      d2log_D += dlog_D * dlog_D_D_new;
+
+      dlog_h += dlog_C + dlog_D;
+      d2log_h += d2log_C + d2log_D;
+      C = C_new;
+      D = D_new;
+      const double delta = C_new * D_new;
+      h *= delta;
+      if (odd) {
+        if (std::fabs(delta - 1.0) < tolerance
+            && std::fabs(d2log_h - d2log_h_prev)
+                   <= tolerance * std::fabs(d2log_h)) {
+          break;
+        }
+        d2log_h_prev = d2log_h;
+      }
+    }
+    const double dlog_I = log1m_x + digamma_ab - digamma_b + dlog_h;
+    const double d2log_Q
+        = trigamma_ab - boost::math::trigamma(b, boost_policy_t<>());
+    return v_b * h * Q * (dlog_I * dlog_I + d2log_Q + d2log_h);
+  } else {
+    double log_C_a = 0.0;
+    double log_C_b = 0.0;
+    double Dv_log_C_a = 0.0;
+    double Dv_log_C_b = 0.0;
+    const double d_1_a = -(x + d_1) / ap1;
+    const double d_1_b = -x / ap1;
+    double log_D_a = -d_1_a * D;
+    double log_D_b = -d_1_b * D;
+    // D_v log D = -D D_v d_1
+    const double d_1_aa = 2.0 * (x + d_1) / (ap1 * ap1);
+    const double d_1_ab = x / (ap1 * ap1);
+    const double Dv_d_1 = d_1_a * v_a + d_1_b * v_b;
+    double Dv_log_D_a
+        = -(d_1_aa * v_a + d_1_ab * v_b) * D - log_D_a * D * Dv_d_1;
+    double Dv_log_D_b = -d_1_ab * v_a * D - log_D_b * D * Dv_d_1;
+    double log_h_a = log_D_a;
+    double log_h_b = log_D_b;
+    double Dv_log_h_a = Dv_log_D_a;
+    double Dv_log_h_b = Dv_log_D_b;
+    double Dv_log_h_a_prev = Dv_log_h_a;
+    double Dv_log_h_b_prev = Dv_log_h_b;
+    for (int n = 2; n <= 2 * max_iterations + 1; ++n) {
+      const int m = n / 2;
+      const bool odd = n % 2 == 1;
+      const double a_plus_n = a + n;
+      const double inv_den = 1.0 / (a_plus_n * (a_plus_n - 1.0));
+      const double den_a = 2.0 * a_plus_n - 1.0;
+      double num_b;
+      double num;
+      double num_a = 0.0;
+      double num_aa = 0.0;
+      double num_ab = 0.0;
+      if (odd) {
+        num_b = -(a + m) * x;
+        num = num_b * (a + b + m);
+        num_a = -(2.0 * (a + m) + b) * x;
+        num_aa = -2.0 * x;
+        num_ab = -x;
+      } else {
+        num_b = m * x;
+        num = num_b * (b - m);
+      }
+      const double d_n = num * inv_den;
+      const double d_n_a = (num_a - d_n * den_a) * inv_den;
+      const double d_n_b = num_b * inv_den;
+      const double d_n_aa
+          = (num_aa - 2.0 * d_n_a * den_a - 2.0 * d_n) * inv_den;
+      const double d_n_ab = (num_ab - d_n_b * den_a) * inv_den;
+      const double Dv_d_n = d_n_a * v_a + d_n_b * v_b;
+      const double Dv_d_n_a = d_n_aa * v_a + d_n_ab * v_b;
+      const double Dv_d_n_b = d_n_ab * v_a;
+
+      // one step of the modified Lentz method
+      // C_new = 1 + d_n / C, with C C_new = C + d_n
+      double C_new = 1.0 + d_n / C;
+      double C_C_new = C + d_n;
+      if (std::fabs(C_new) < small) {
+        C_new = small;
+        C_C_new = C * small;
+      }
+      const double inv_C_C_new = 1.0 / C_C_new;
+      const double Dv_log_C = log_C_a * v_a + log_C_b * v_b;
+      const double Dv_log_C_C_new = (C * Dv_log_C + Dv_d_n) * inv_C_C_new;
+      Dv_log_C_a
+          = (Dv_d_n_a - Dv_d_n * log_C_a - d_n * Dv_log_C_a) * inv_C_C_new;
+      log_C_a = (d_n_a - d_n * log_C_a) * inv_C_C_new;
+      Dv_log_C_a -= log_C_a * Dv_log_C_C_new;
+      Dv_log_C_b
+          = (Dv_d_n_b - Dv_d_n * log_C_b - d_n * Dv_log_C_b) * inv_C_C_new;
+      log_C_b = (d_n_b - d_n * log_C_b) * inv_C_C_new;
+      Dv_log_C_b -= log_C_b * Dv_log_C_C_new;
+
+      // D_new = 1 / (1 + d_n D)
+      double D_new = 1.0 + d_n * D;
+      if (std::fabs(D_new) < small) {
+        D_new = small;
+      }
+      D_new = 1.0 / D_new;
+      const double D_D_new = D * D_new;
+      const double Dv_log_D = log_D_a * v_a + log_D_b * v_b;
+      const double Dv_log_D_D_new
+          = Dv_log_D - D_D_new * (Dv_d_n + d_n * Dv_log_D);
+      Dv_log_D_a = -(Dv_d_n_a + Dv_d_n * log_D_a + d_n * Dv_log_D_a) * D_D_new;
+      log_D_a = -(d_n_a + d_n * log_D_a) * D_D_new;
+      Dv_log_D_a += log_D_a * Dv_log_D_D_new;
+      Dv_log_D_b = -(Dv_d_n_b + Dv_d_n * log_D_b + d_n * Dv_log_D_b) * D_D_new;
+      log_D_b = -(d_n_b + d_n * log_D_b) * D_D_new;
+      Dv_log_D_b += log_D_b * Dv_log_D_D_new;
+
+      log_h_a += log_C_a + log_D_a;
+      log_h_b += log_C_b + log_D_b;
+      Dv_log_h_a += Dv_log_C_a + Dv_log_D_a;
+      Dv_log_h_b += Dv_log_C_b + Dv_log_D_b;
+      C = C_new;
+      D = D_new;
+      const double delta = C_new * D_new;
+      h *= delta;
+      if (odd) {
+        if (std::fabs(delta - 1.0) < tolerance
+            && std::fabs(Dv_log_h_a - Dv_log_h_a_prev)
+                   <= tolerance * std::fabs(Dv_log_h_a)
+            && std::fabs(Dv_log_h_b - Dv_log_h_b_prev)
+                   <= tolerance * std::fabs(Dv_log_h_b)) {
+          break;
+        }
+        Dv_log_h_a_prev = Dv_log_h_a;
+        Dv_log_h_b_prev = Dv_log_h_b;
+      }
+    }
+    const double I = h * Q;
+    const double log_I_a = log_x - 1.0 / a + digamma_ab - digamma_a + log_h_a;
+    const double log_I_b = log1m_x + digamma_ab - digamma_b + log_h_b;
+    const double Dv_log_I = log_I_a * v_a + log_I_b * v_b;
+    // the Hessian of log Q is that of -log B(a, b) plus diag(1 / a^2, 0)
+    const double Dv_log_Q_a
+        = (trigamma_ab - boost::math::trigamma(a, boost_policy_t<>())
+           + 1.0 / (a * a))
+              * v_a
+          + trigamma_ab * v_b;
+    const double Dv_log_Q_b
+        = trigamma_ab * v_a
+          + (trigamma_ab - boost::math::trigamma(b, boost_policy_t<>())) * v_b;
+    return std::pair{I * (log_I_a * Dv_log_I + Dv_log_Q_a + Dv_log_h_a),
+                     I * (log_I_b * Dv_log_I + Dv_log_Q_b + Dv_log_h_b)};
+  }
+}
+
+/**
+ * Hessian-vector product H v of the regularized incomplete beta function
+ * I_z(a, b), where H is the Hessian of I_z(a, b) in (a, b). When v holds the
+ * adjoints of (d/da I_z(a, b), d/db I_z(a, b)), H v is what the reverse pass
+ * adds to the adjoints of a and b.
+ *
+ * The power series or the continued fraction, and whether to apply the
+ * symmetry relation, are chosen by the same rule as in grad_reg_inc_beta.
+ * Under the symmetry relation I_z(a, b) = 1 - I_{1-z}(b, a), let H' be the
+ * Hessian of I_{1-z}(b, a) in its own arguments (b, a), and let P swap the
+ * two components of a vector. Then H = -P H' P, so H v = -P H' (P v). In
+ * particular H_aa = -H'_bb and H_bb = -H'_aa.
+ *
+ * @tparam grads <code>G1</code> for H_aa v_a, <code>G2</code> for H_bb v_b,
+ *   <code>BOTH</code> for H v
+ * @param a first shape
+ * @param b second shape
+ * @param z upper limit, in [0, 1]
+ * @param v_a a component of v, not used for <code>G2</code>
+ * @param v_b b component of v, not used for <code>G1</code>
+ * @param digamma_a digamma(a), not used for <code>G2</code>
+ * @param digamma_b digamma(b), not used for <code>G1</code>
+ * @param digamma_ab digamma(a + b)
+ * @return H_aa v_a for <code>G1</code>, H_bb v_b for <code>G2</code>,
+ *   ((H v)_a, (H v)_b) for <code>BOTH</code>. The result is 0 at z = 0 and
+ *   z = 1.
+ */
+template <GRAD_RIB grads>
+inline std::conditional_t<grads == GRAD_RIB::BOTH, std::pair<double, double>,
+                          double>
+inc_beta_hvp(double a, double b, double z, double v_a, double v_b,
+             double digamma_a, double digamma_b, double digamma_ab) {
+  if (z == 0.0 || z == 1.0) {
+    return {};
+  }
+  // the symmetry relation swaps which shape is requested
+  constexpr GRAD_RIB swapped_grads
+      = grads == GRAD_RIB::G1
+            ? GRAD_RIB::G2
+            : grads == GRAD_RIB::G2 ? GRAD_RIB::G1 : GRAD_RIB::BOTH;
+  const double mean = a / (a + b);
+  // the same branch rule as grad_reg_inc_beta
+  const bool use_power_series
+      = (z < mean && b * z <= 1.0 && z <= 0.5)
+        || (z >= mean && a * (1.0 - z) <= 1.0 && z >= 0.5);
+  if (use_power_series) {
+    if (z >= mean) {
+      const auto hv = inc_beta_hvp_power_series<swapped_grads>(
+          b, a, 1.0 - z, v_b, v_a, digamma_b, digamma_a, digamma_ab);
+      if constexpr (grads != GRAD_RIB::BOTH) {
+        return -hv;
+      } else {
+        return {-hv.second, -hv.first};
+      }
+    } else {
+      return inc_beta_hvp_power_series<grads>(a, b, z, v_a, v_b, digamma_a,
+                                              digamma_b, digamma_ab);
+    }
+  } else {
+    // the same threshold as grad_reg_inc_beta, so both take the same branch
+    if (z >= (a + 1.0) / (a + b + 2.0)) {
+      const auto hv = inc_beta_hvp_continued_fraction<swapped_grads>(
+          b, a, 1.0 - z, v_b, v_a, digamma_b, digamma_a, digamma_ab);
+      if constexpr (grads != GRAD_RIB::BOTH) {
+        return -hv;
+      } else {
+        return {-hv.second, -hv.first};
+      }
+    } else {
+      return inc_beta_hvp_continued_fraction<grads>(
+          a, b, z, v_a, v_b, digamma_a, digamma_b, digamma_ab);
+    }
+  }
+}
+
+}  // namespace internal
+
+/**
+ * Computes the requested partial derivatives of the regularized incomplete
+ * beta function I_z(a, b) = <code>inc_beta(a, b, z)</code> with respect to
+ * <code>a</code> and <code>b</code>, for <code>var</code> arguments.
+ *
+ * The values come from the prim overload, evaluated in double. The reverse
+ * pass does not differentiate through that series or continued fraction.
+ * Instead, with v the adjoints of the outputs, it adds H v to the adjoints of
+ * a and b, where H is the Hessian of I_z(a, b) in (a, b), computed directly
+ * by internal::inc_beta_hvp. To the adjoint of z it adds v times the
+ * closed-form mixed derivatives
+ *
+ *   d/dz d/da I_z(a, b) = inc_beta_ddz(a, b, z) (log z - psi(a) + psi(a+b)),
+ *   d/dz d/db I_z(a, b) = inc_beta_ddz(a, b, z) (log1p(-z) - psi(b) +
+ *                         psi(a+b)),
+ *
+ * where psi is the digamma function.
+ *
+ * The digamma arguments are passed only so that the values are cheap to
+ * compute, and they must equal digamma(a), digamma(b), and digamma(a + b).
+ * They get no adjoint. The a and b adjoints above are total derivatives, so
+ * they already include any dependence on a and b that runs through the
+ * digamma values.
+ *
+ * A one-sided request still needs the cross term of H when both shapes are
+ * vars. For <code>G1</code> the adjoint of b is H_ab v_a, which the
+ * <code>BOTH</code> product gives with v_b = 0. That product needs
+ * digamma(b), which a <code>G1</code> call does not receive, so the reverse
+ * pass computes it when b is a var. When only a is a var, H_aa v_a does not
+ * need digamma(b). <code>G2</code> is the same with a and b swapped.
+ *
+ * For <code>BOTH</code> there are two outputs, and one reverse pass callback
+ * reads both of their adjoints. It runs after every use of the outputs has
+ * added into those adjoints.
+ *
+ * Each argument is a <code>var</code> or arithmetic. Calls that mix in
+ * <code>fvar</code> go to the prim overload.
+ *
+ * @tparam grads_requested <code>G1</code> for the derivative with respect to
+ *   a, <code>G2</code> for the derivative with respect to b,
+ *   <code>BOTH</code> for both
+ * @tparam T1 type of a
+ * @tparam T2 type of b
+ * @tparam T3 type of z
+ * @tparam T4 type of digammaA
+ * @tparam T5 type of digammaB
+ * @tparam T6 type of digammaSum
+ * @tparam T7 type of betaAB
+ * @param a a
+ * @param b b
+ * @param z z
+ * @param digammaA the value of <code>digamma(a)</code>, not used for
+ *   <code>G2</code>
+ * @param digammaB the value of <code>digamma(b)</code>, not used for
+ *   <code>G1</code>
+ * @param digammaSum the value of <code>digamma(a + b)</code>
+ * @param betaAB accepted for interface compatibility and not used
+ * @return d/da I_z(a, b) for <code>G1</code>, d/db for <code>G2</code>, and
+ *   the pair (d/da, d/db) for <code>BOTH</code>
+ */
+template <GRAD_RIB grads_requested, typename T1, typename T2, typename T3,
+          typename T4, typename T5, typename T6, typename T7,
+          require_any_var_t<T1, T2, T3, T4, T5, T6, T7>* = nullptr>
+inline auto grad_reg_inc_beta(const T1& a, const T2& b, const T3& z,
+                              const T4& digammaA, const T5& digammaB,
+                              const T6& digammaSum, const T7& /* betaAB */) {
+  const double digamma_a = value_of(digammaA);
+  const double digamma_b = value_of(digammaB);
+  const double digamma_ab = value_of(digammaSum);
+  if constexpr (grads_requested == GRAD_RIB::G1) {
+    return make_callback_var(
+        grad_reg_inc_beta<grads_requested>(value_of(a), value_of(b),
+                                           value_of(z), digamma_a, digamma_b,
+                                           digamma_ab, 0.0),
+        [a, b, z, digamma_a, digamma_ab](auto& vi) {
+          const double a_val = value_of(a);
+          const double b_val = value_of(b);
+          const double z_val = value_of(z);
+          if constexpr (is_var_v<T1> && is_var_v<T2>) {
+            const auto [hv_a, hv_b] = internal::inc_beta_hvp<GRAD_RIB::BOTH>(
+                a_val, b_val, z_val, vi.adj(), 0.0, digamma_a, digamma(b_val),
+                digamma_ab);
+            a.adj() += hv_a;
+            b.adj() += hv_b;
+          } else if constexpr (is_var_v<T1>) {
+            // H_aa v_a does not use digamma(b)
+            a.adj() += internal::inc_beta_hvp<GRAD_RIB::G1>(
+                a_val, b_val, z_val, vi.adj(), 0.0, digamma_a, 0.0, digamma_ab);
+          } else if constexpr (is_var_v<T2>) {
+            b.adj() += internal::inc_beta_hvp<GRAD_RIB::BOTH>(
+                           a_val, b_val, z_val, vi.adj(), 0.0, digamma_a,
+                           digamma(b_val), digamma_ab)
+                           .second;
+          }
+          if constexpr (is_var_v<T3>) {
+            if (z_val > 0.0 && z_val < 1.0) {
+              const double log_term = std::log(z_val) - digamma_a + digamma_ab;
+              z.adj()
+                  += vi.adj() * inc_beta_ddz(a_val, b_val, z_val) * log_term;
+            }
+          }
+        });
+  } else if constexpr (grads_requested == GRAD_RIB::G2) {
+    return make_callback_var(
+        grad_reg_inc_beta<grads_requested>(value_of(a), value_of(b),
+                                           value_of(z), digamma_a, digamma_b,
+                                           digamma_ab, 0.0),
+        [a, b, z, digamma_b, digamma_ab](auto& vi) {
+          const double a_val = value_of(a);
+          const double b_val = value_of(b);
+          const double z_val = value_of(z);
+          if constexpr (is_var_v<T1> && is_var_v<T2>) {
+            const auto [hv_a, hv_b] = internal::inc_beta_hvp<GRAD_RIB::BOTH>(
+                a_val, b_val, z_val, 0.0, vi.adj(), digamma(a_val), digamma_b,
+                digamma_ab);
+            a.adj() += hv_a;
+            b.adj() += hv_b;
+          } else if constexpr (is_var_v<T1>) {
+            a.adj() += internal::inc_beta_hvp<GRAD_RIB::BOTH>(
+                           a_val, b_val, z_val, 0.0, vi.adj(), digamma(a_val),
+                           digamma_b, digamma_ab)
+                           .first;
+          } else if constexpr (is_var_v<T2>) {
+            // H_bb v_b does not use digamma(a)
+            b.adj() += internal::inc_beta_hvp<GRAD_RIB::G2>(
+                a_val, b_val, z_val, 0.0, vi.adj(), 0.0, digamma_b, digamma_ab);
+          }
+          if constexpr (is_var_v<T3>) {
+            if (z_val > 0.0 && z_val < 1.0) {
+              const double log_term = log1p(-z_val) - digamma_b + digamma_ab;
+              z.adj()
+                  += vi.adj() * inc_beta_ddz(a_val, b_val, z_val) * log_term;
+            }
+          }
+        });
+  } else {  // both
+    const auto [g1_val, g2_val] = grad_reg_inc_beta<GRAD_RIB::BOTH>(
+        value_of(a), value_of(b), value_of(z), digamma_a, digamma_b, digamma_ab,
+        0.0);
+    var g1 = g1_val;
+    var g2 = g2_val;
+    reverse_pass_callback([a, b, z, g1, g2, digamma_a, digamma_b,
+                           digamma_ab]() {
+      const double a_val = value_of(a);
+      const double b_val = value_of(b);
+      const double z_val = value_of(z);
+      const double v_a = g1.adj();
+      const double v_b = g2.adj();
+      if constexpr (is_var_v<T1> || is_var_v<T2>) {
+        const auto [hv_a, hv_b] = internal::inc_beta_hvp<GRAD_RIB::BOTH>(
+            a_val, b_val, z_val, v_a, v_b, digamma_a, digamma_b, digamma_ab);
+        if constexpr (is_var_v<T1>) {
+          a.adj() += hv_a;
+        }
+        if constexpr (is_var_v<T2>) {
+          b.adj() += hv_b;
+        }
+      }
+      if constexpr (is_var_v<T3>) {
+        if (z_val > 0.0 && z_val < 1.0) {
+          z.adj() += inc_beta_ddz(a_val, b_val, z_val)
+                     * (v_a * (std::log(z_val) - digamma_a + digamma_ab)
+                        + v_b * (log1p(-z_val) - digamma_b + digamma_ab));
+        }
+      }
+    });
+    return std::pair{g1, g2};
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/rev/fun/inc_beta.hpp b/stan/math/rev/fun/inc_beta.hpp
index be8a65b7b6a..037f3f07433 100644
--- a/stan/math/rev/fun/inc_beta.hpp
+++ b/stan/math/rev/fun/inc_beta.hpp
@@ -3,41 +3,60 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
-#include <stan/math/rev/fun/beta.hpp>
-#include <stan/math/rev/fun/digamma.hpp>
-#include <stan/math/rev/fun/pow.hpp>
+#include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
-#include <cmath>
+#include <stan/math/prim/fun/inc_beta_ddz.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
 
 namespace stan {
 namespace math {
 
-namespace internal {
-
-class inc_beta_vvv_vari : public op_vvv_vari {
- public:
-  inc_beta_vvv_vari(vari* avi, vari* bvi, vari* cvi)
-      : op_vvv_vari(inc_beta(avi->val_, bvi->val_, cvi->val_), avi, bvi, cvi) {}
-  void chain() {
-    double d_a;
-    double d_b;
-    const double beta_ab = beta(avi_->val_, bvi_->val_);
-    grad_reg_inc_beta(d_a, d_b, avi_->val_, bvi_->val_, cvi_->val_,
-                      digamma(avi_->val_), digamma(bvi_->val_),
-                      digamma(avi_->val_ + bvi_->val_), beta_ab);
-
-    avi_->adj_ += adj_ * d_a;
-    bvi_->adj_ += adj_ * d_b;
-    cvi_->adj_ += adj_ * std::pow(1 - cvi_->val_, bvi_->val_ - 1)
-                  * std::pow(cvi_->val_, avi_->val_ - 1) / beta_ab;
-  }
-};
-
-}  // namespace internal
-
-inline var inc_beta(const var& a, const var& b, const var& c) {
-  return var(new internal::inc_beta_vvv_vari(a.vi_, b.vi_, c.vi_));
+/**
+ * The regularized incomplete beta function I_z(a, b) for <code>var</code>
+ * arguments. Only the partials of the <code>var</code> arguments are
+ * computed in the reverse pass. Each argument is a <code>var</code> or
+ * arithmetic. Mixing with <code>fvar</code> is not supported.
+ *
+ * @tparam T1 type of a
+ * @tparam T2 type of b
+ * @tparam T3 type of z
+ * @param a first shape, positive
+ * @param b second shape, positive
+ * @param z upper limit, in [0, 1]
+ * @return I_z(a, b)
+ */
+template <typename T1, typename T2, typename T3,
+          require_all_stan_scalar_t<T1, T2, T3>* = nullptr,
+          require_all_var_or_arithmetic_t<T1, T2, T3>* = nullptr,
+          require_any_var_t<T1, T2, T3>* = nullptr>
+inline var inc_beta(const T1& a, const T2& b, const T3& z) {
+  return make_callback_var(
+      inc_beta(value_of(a), value_of(b), value_of(z)), [a, b, z](auto& vi) {
+        const double a_val = value_of(a);
+        const double b_val = value_of(b);
+        const double z_val = value_of(z);
+        if constexpr (is_autodiff_v<T1> && is_autodiff_v<T2>) {
+          const auto [d_a, d_b] = grad_reg_inc_beta<GRAD_RIB::BOTH>(
+              a_val, b_val, z_val, digamma(a_val), digamma(b_val),
+              digamma(a_val + b_val), 0.0);
+          a.adj() += vi.adj() * d_a;
+          b.adj() += vi.adj() * d_b;
+        } else if constexpr (is_autodiff_v<T1>) {
+          a.adj() += vi.adj()
+                     * grad_reg_inc_beta<GRAD_RIB::G1>(
+                         a_val, b_val, z_val, digamma(a_val), 0.0,
+                         digamma(a_val + b_val), 0.0);
+        } else if constexpr (is_autodiff_v<T2>) {
+          b.adj() += vi.adj()
+                     * grad_reg_inc_beta<GRAD_RIB::G2>(
+                         a_val, b_val, z_val, 0.0, digamma(b_val),
+                         digamma(a_val + b_val), 0.0);
+        }
+        if constexpr (is_autodiff_v<T3>) {
+          z.adj() += vi.adj() * inc_beta_ddz(a_val, b_val, z_val);
+        }
+      });
 }
 
 }  // namespace math
diff --git a/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp b/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp
index e6613d43ea1..7831ab106ee 100644
--- a/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp
+++ b/test/unit/math/mix/fun/grad_reg_inc_beta_test.cpp
@@ -1,145 +1,95 @@
-#include <stan/math/mix.hpp>
-#include <gtest/gtest.h>
-#include <test/unit/math/rev/fun/util.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <vector>
 
-TEST(ProbInternalMath, grad_reg_inc_beta_fv) {
-  using stan::math::digamma;
-  using stan::math::exp;
-  using stan::math::fvar;
-  using stan::math::lbeta;
-  using stan::math::var;
+namespace {
 
-  fvar<var> a = 1.0;
-  fvar<var> b = 1.0;
-  fvar<var> g = 0.4;
-  a.d_ = 1.0;
-  b.d_ = 1.0;
-  g.d_ = 1.0;
-  fvar<var> dig_a = digamma(a);
-  fvar<var> dig_b = digamma(b);
-  fvar<var> dig_sum = digamma(a + b);
-  fvar<var> beta_ab = exp(lbeta(a, b));
-  fvar<var> g_a;
-  fvar<var> g_b;
-
-  stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum,
-                                beta_ab);
-  EXPECT_FLOAT_EQ(-0.36651629442883944183907601651838247842001142107486495485,
-                  g_a.val_.val());
-  EXPECT_NEAR(0.306495375042422864944011633197968575202046200428315551199,
-              g_b.val_.val(), 1e-6);
+/**
+ * Runs expect_ad at points in each region of grad_reg_inc_beta:
+ *   (1, 1, 0.4), (2, 3, 0.25), (0.5, 0.5, 0.3): power series
+ *   (10, 1, 0.5): power series, integer b
+ *   (2, 3, 0.6): power series through the symmetry relation
+ *   (20, 20, 0.3), (5, 30, 0.1): continued fraction
+ *   (20, 20, 0.7), (5, 30, 0.3): continued fraction through the symmetry
+ *     relation
+ *   (10, 1, 0.7), (4, 4, 0.5): continued fraction, integer b, where the
+ *     partial numerator d_{2b} is 0 but its b derivative is not
+ */
+template <typename F>
+void expect_ad_grad_reg_inc_beta(const F& f) {
+  stan::test::expect_ad(f, 1.0, 1.0, 0.4);
+  stan::test::expect_ad(f, 2.0, 3.0, 0.25);
+  stan::test::expect_ad(f, 0.5, 0.5, 0.3);
+  stan::test::expect_ad(f, 10.0, 1.0, 0.5);
+  stan::test::expect_ad(f, 2.0, 3.0, 0.6);
+  stan::test::expect_ad(f, 20.0, 20.0, 0.3);
+  stan::test::expect_ad(f, 5.0, 30.0, 0.1);
+  stan::test::expect_ad(f, 20.0, 20.0, 0.7);
+  stan::test::expect_ad(f, 5.0, 30.0, 0.3);
+  stan::test::expect_ad(f, 10.0, 1.0, 0.7);
+  stan::test::expect_ad(f, 4.0, 4.0, 0.5);
 }
-TEST(ProbInternalMath, grad_reg_inc_beta_fv_1stDeriv1) {
-  using stan::math::digamma;
-  using stan::math::exp;
-  using stan::math::fvar;
-  using stan::math::lbeta;
-  using stan::math::var;
-
-  fvar<var> a = 1.0;
-  fvar<var> b = 1.0;
-  fvar<var> g = 0.4;
-  a.d_ = 1.0;
-  b.d_ = 0.0;
-  g.d_ = 0.0;
-  fvar<var> dig_a = digamma(a);
-  fvar<var> dig_b = digamma(b);
-  fvar<var> dig_sum = digamma(a + b);
-  fvar<var> beta_ab = exp(lbeta(a, b));
-  fvar<var> g_a;
-  fvar<var> g_b;
 
-  stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum,
-                                beta_ab);
+}  // namespace
 
-  std::vector<stan::math::var> y1{a.val_};
-  std::vector<double> grad1;
-  g_a.val_.grad(y1, grad1);
-  EXPECT_FLOAT_EQ(0.33583548212738989400284958327902414335945097838423129,
-                  grad1[0]);
+TEST(mathMixScalFun, grad_reg_inc_beta_g1) {
+  auto f = [](const auto& a, const auto& b, const auto& z) {
+    using stan::math::digamma;
+    return stan::math::grad_reg_inc_beta<stan::math::GRAD_RIB::G1>(
+        a, b, z, digamma(a), 0.0, digamma(a + b), 0.0);
+  };
+  expect_ad_grad_reg_inc_beta(f);
 }
-TEST(ProbInternalMath, grad_reg_inc_beta_fv_1stDeriv2) {
-  using stan::math::digamma;
-  using stan::math::exp;
-  using stan::math::fvar;
-  using stan::math::lbeta;
-  using stan::math::var;
-
-  fvar<var> a = 1.0;
-  fvar<var> b = 1.0;
-  fvar<var> g = 0.4;
-  a.d_ = 0.0;
-  b.d_ = 1.0;
-  g.d_ = 0.0;
-  fvar<var> dig_a = digamma(a);
-  fvar<var> dig_b = digamma(b);
-  fvar<var> dig_sum = digamma(a + b);
-  fvar<var> beta_ab = exp(lbeta(a, b));
-  fvar<var> g_a;
-  fvar<var> g_b;
-
-  stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum,
-                                beta_ab);
 
-  std::vector<stan::math::var> y1{b.val_};
-  std::vector<double> grad1;
-  g_b.val_.grad(y1, grad1);
-  EXPECT_NEAR(-0.156565690737548079304827886, grad1[0], 1e-6);
+TEST(mathMixScalFun, grad_reg_inc_beta_g2) {
+  auto f = [](const auto& a, const auto& b, const auto& z) {
+    using stan::math::digamma;
+    return stan::math::grad_reg_inc_beta<stan::math::GRAD_RIB::G2>(
+        a, b, z, 0.0, digamma(b), digamma(a + b), 0.0);
+  };
+  expect_ad_grad_reg_inc_beta(f);
 }
-TEST(ProbInternalMath, grad_reg_inc_beta_fv_2ndDeriv1) {
-  using stan::math::digamma;
-  using stan::math::exp;
-  using stan::math::fvar;
-  using stan::math::lbeta;
-  using stan::math::var;
 
-  fvar<var> a = 1.0;
-  fvar<var> b = 1.0;
-  fvar<var> g = 0.4;
-  a.d_ = 1.0;
-  b.d_ = 0.0;
-  g.d_ = 0.0;
-  fvar<var> dig_a = digamma(a);
-  fvar<var> dig_b = digamma(b);
-  fvar<var> dig_sum = digamma(a + b);
-  fvar<var> beta_ab = exp(lbeta(a, b));
-  fvar<var> g_a;
-  fvar<var> g_b;
-
-  stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum,
-                                beta_ab);
-
-  std::vector<stan::math::var> y1{a.val_};
-  std::vector<double> grad1;
-  g_a.d_.grad(y1, grad1);
-  EXPECT_FLOAT_EQ(-0.30772293970781581317390510390046098438962772318921,
-                  grad1[0]);
+TEST(mathMixScalFun, grad_reg_inc_beta_both) {
+  auto f = [](const auto& a, const auto& b, const auto& z) {
+    using stan::math::digamma;
+    auto grads = stan::math::grad_reg_inc_beta<stan::math::GRAD_RIB::BOTH>(
+        a, b, z, digamma(a), digamma(b), digamma(a + b), 0.0);
+    Eigen::Matrix<decltype(grads.first), Eigen::Dynamic, 1> ret(2);
+    ret << grads.first, grads.second;
+    return ret;
+  };
+  expect_ad_grad_reg_inc_beta(f);
 }
-TEST(ProbInternalMath, grad_reg_inc_beta_fv_2ndDeriv2) {
-  using stan::math::digamma;
-  using stan::math::exp;
-  using stan::math::fvar;
-  using stan::math::lbeta;
-  using stan::math::var;
 
-  fvar<var> a = 1.0;
-  fvar<var> b = 1.0;
-  fvar<var> g = 0.4;
-  a.d_ = 0.0;
-  b.d_ = 1.0;
-  g.d_ = 0.0;
-  fvar<var> dig_a = digamma(a);
-  fvar<var> dig_b = digamma(b);
-  fvar<var> dig_sum = digamma(a + b);
-  fvar<var> beta_ab = exp(lbeta(a, b));
-  fvar<var> g_a;
-  fvar<var> g_b;
-
-  stan::math::grad_reg_inc_beta(g_a, g_b, a, b, g, dig_a, dig_b, dig_sum,
-                                beta_ab);
+// the student_t_cdf call: a = nu / 2, b = 1 / 2 fixed, z = 1 - r
+TEST(mathMixScalFun, grad_reg_inc_beta_student) {
+  auto f = [](const auto& a, const auto& /* b */, const auto& z) {
+    using stan::math::digamma;
+    return stan::math::grad_reg_inc_beta<stan::math::GRAD_RIB::G1>(
+        0.5 * a, 0.5, 1.0 - z, digamma(0.5 * a), 0.0, digamma(0.5 + 0.5 * a),
+        0.0);
+  };
+  expect_ad_grad_reg_inc_beta(f);
+}
 
-  std::vector<stan::math::var> y1{b.val_};
-  std::vector<double> grad1;
-  g_b.d_.grad(y1, grad1);
-  EXPECT_NEAR(0.079977766631361187517939795, grad1[0], 1e-4);
+/**
+ * The overload writing both partials through out-parameters takes every
+ * argument as the same type.
+ */
+TEST(mathMixScalFun, grad_reg_inc_beta_out_params) {
+  auto f = [](const auto& a, const auto& b, const auto& z) {
+    using stan::math::digamma;
+    using T = stan::return_type_t<decltype(a), decltype(b), decltype(z)>;
+    const T a_t = a;
+    const T b_t = b;
+    const T z_t = z;
+    T g1 = 0;
+    T g2 = 0;
+    stan::math::grad_reg_inc_beta(g1, g2, a_t, b_t, z_t, digamma(a_t),
+                                  digamma(b_t), digamma(a_t + b_t), T(0));
+    Eigen::Matrix<T, Eigen::Dynamic, 1> ret(2);
+    ret << g1, g2;
+    return ret;
+  };
+  expect_ad_grad_reg_inc_beta(f);
 }
diff --git a/test/unit/math/mix/fun/inc_beta_test.cpp b/test/unit/math/mix/fun/inc_beta_test.cpp
index 4f00deab729..8805333dddd 100644
--- a/test/unit/math/mix/fun/inc_beta_test.cpp
+++ b/test/unit/math/mix/fun/inc_beta_test.cpp
@@ -28,6 +28,29 @@ TEST(mathMixScalFun, inc_beta) {
   stan::test::expect_ad(f, nan, nan, nan);
 }
 
+/**
+ * At z = 1 with b < 1 the density is infinite. The forward and reverse
+ * sweeps of fvar<var> must not throw there.
+ */
+TEST(mathMixScalFun, inc_beta_upper_bound_infinite_density) {
+  using stan::math::fvar;
+  using stan::math::var;
+  // grad() runs every node on the stack; start from an empty one
+  stan::math::recover_memory();
+  var a_v = 1.05;
+  var b_v = 0.35;
+  fvar<var> a(a_v, 1.0);
+  fvar<var> b(b_v, 0.0);
+  fvar<var> z(1.0, 0.0);
+  fvar<var> f;
+  EXPECT_NO_THROW(f = stan::math::inc_beta(a, b, z));
+  EXPECT_EQ(1.0, f.val_.val());
+  EXPECT_NO_THROW(f.val_.grad());
+  stan::math::set_zero_all_adjoints();
+  EXPECT_NO_THROW(f.d_.grad());
+  stan::math::recover_memory();
+}
+
 TEST(mathMixScalFun, inc_beta_vec) {
   auto f = [](const auto& x1, const auto& x2, const auto& x3) {
     return stan::math::inc_beta(x1, x2, x3);
diff --git a/test/unit/math/rev/fun/inc_beta_test.cpp b/test/unit/math/rev/fun/inc_beta_test.cpp
new file mode 100644
index 00000000000..a66e629482a
--- /dev/null
+++ b/test/unit/math/rev/fun/inc_beta_test.cpp
@@ -0,0 +1,62 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <vector>
+
+namespace {
+
+/**
+ * Shape gradients of inc_beta against a high-precision reference (mpmath,
+ * 80 digits), at points where a finite-difference check is not informative.
+ */
+void expect_inc_beta_grads(double a, double b, double z, double dda, double ddb,
+                           double rtol) {
+  using stan::math::var;
+  var av = a;
+  var bv = b;
+  var zv = z;
+  var f = stan::math::inc_beta(av, bv, zv);
+  std::vector<var> x{av, bv};
+  std::vector<double> g;
+  f.grad(x, g);
+  EXPECT_TRUE(std::isfinite(g[0])) << "a=" << a << " b=" << b << " z=" << z;
+  EXPECT_TRUE(std::isfinite(g[1])) << "a=" << a << " b=" << b << " z=" << z;
+  EXPECT_NEAR(dda, g[0], rtol * std::fabs(dda))
+      << "a=" << a << " b=" << b << " z=" << z;
+  EXPECT_NEAR(ddb, g[1], rtol * std::fabs(ddb))
+      << "a=" << a << " b=" << b << " z=" << z;
+}
+
+}  // namespace
+
+TEST_F(AgradRev, inc_beta_shape_gradients) {
+  expect_inc_beta_grads(2, 3, 0.25, -0.22805360232434637, 0.10692153005229139,
+                        1e-12);
+  // a + b above the beta(a, b) underflow limit
+  expect_inc_beta_grads(600, 600, 0.5, -0.01152047147304755,
+                        0.01152047147304755, 1e-12);
+  // z^a (1 - z)^b underflows while I is 1
+  expect_inc_beta_grads(352, 590, 0.758, -3.3772717797123242e-138,
+                        4.5390666074606649e-138, 1e-12);
+  // integer b
+  expect_inc_beta_grads(10, 1, 0.5, -0.00067690154351557159,
+                        0.00226574342615196, 1e-12);
+}
+
+/**
+ * At z = 1, I is 1 for all shapes, so the shape gradients are 0. With b < 1
+ * the density, which is the z gradient, is infinite.
+ */
+TEST_F(AgradRev, inc_beta_upper_bound_infinite_density) {
+  using stan::math::var;
+  var a = 1.05;
+  var b = 0.35;
+  var z = 1.0;
+  var f = stan::math::inc_beta(a, b, z);
+  EXPECT_NO_THROW(f.grad());
+  EXPECT_EQ(1.0, f.val());
+  EXPECT_EQ(0.0, a.adj());
+  EXPECT_EQ(0.0, b.adj());
+  EXPECT_TRUE(std::isinf(z.adj()) && z.adj() > 0);
+}