Fold grad_inc_beta into grad_reg_inc_beta

stan/math/fwd/fun.hpp

@@ -40,7 +40,6 @@
 #include <stan/math/fwd/fun/fmod.hpp>
 #include <stan/math/fwd/fun/gamma_p.hpp>
 #include <stan/math/fwd/fun/gamma_q.hpp>
-#include <stan/math/fwd/fun/grad_inc_beta.hpp>
 #include <stan/math/fwd/fun/hypergeometric_1F0.hpp>
 #include <stan/math/fwd/fun/hypergeometric_2F1.hpp>
 #include <stan/math/fwd/fun/hypergeometric_pFq.hpp>

stan/math/prim/fun.hpp

@@ -116,7 +116,6 @@
 #include <stan/math/prim/fun/gp_periodic_cov.hpp>
 #include <stan/math/prim/fun/grad_2F1.hpp>
 #include <stan/math/prim/fun/grad_F32.hpp>
-#include <stan/math/prim/fun/grad_inc_beta.hpp>
 #include <stan/math/prim/fun/grad_pFq.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_beta.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>

stan/math/prim/fun/grad_reg_inc_beta.hpp

@@ -2,9 +2,8 @@
 #define STAN_MATH_PRIM_FUN_GRAD_REG_INC_BETA_HPP
 
 #include <stan/math/prim/meta.hpp>
-#include <stan/math/prim/fun/beta.hpp>
-#include <stan/math/prim/fun/grad_inc_beta.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
+#include <stan/math/prim/fun/grad_2F1.hpp>
 #include <cmath>
 
 namespace stan {
@@ -16,6 +15,9 @@ namespace math {
  * <code>ibeta(a, b, z)</code>, with respect to the arguments
  * <code>a</code> and <code>b</code>.
  *
+ * Uses the equivalence to a hypergeometric function. See
+ * http://dlmf.nist.gov/8.17#ii
+ *
  * @tparam T type of arguments
  * @param[out] g1 partial derivative of <code>ibeta(a, b, z)</code>
  * with respect to <code>a</code>
@@ -34,12 +36,25 @@ void grad_reg_inc_beta(T& g1, T& g2, const T& a, const T& b, const T& z,
                        const T& digammaA, const T& digammaB,
                        const T& digammaSum, const T& betaAB) {
   using std::exp;
-  T dBda = 0;
-  T dBdb = 0;
-  grad_inc_beta(dBda, dBdb, a, b, z);
-  T b1 = beta(a, b) * inc_beta(a, b, z);
-  g1 = (dBda - b1 * (digammaA - digammaSum)) / betaAB;
-  g2 = (dBdb - b1 * (digammaB - digammaSum)) / betaAB;
+
+  T c1 = log(z);
+  T c2 = log1m(z);
+  T c3 = betaAB * inc_beta(a, b, z);
+  T C = exp(a * c1 + b * c2) / a;
+  T dF1 = 0;
+  T dF2 = 0;
+  T dF3 = 0;
+  T dFz = 0;
+  if (value_of_rec(C)) {
+    std::forward_as_tuple(dF1, dF2, dF3, dFz)
+        = grad_2F1<true>(a + b, 1.0, a + 1, z);
+  }
+
+  T dBda = (c1 - 1.0 / a) * c3 + C * (dF1 + dF3);
+  T dBdb = c2 * c3 + C * dF1;
+
+  g1 = (dBda - c3 * (digammaA - digammaSum)) / betaAB;
+  g2 = (dBdb - c3 * (digammaB - digammaSum)) / betaAB;
 }
 
 }  // namespace math

stan/math/rev/fun.hpp

@@ -70,7 +70,6 @@
 #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_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>