diff --git a/stan/math/fwd/scal/fun/rising_factorial.hpp b/stan/math/fwd/scal/fun/rising_factorial.hpp
index 2bbba9611d2..f6bf1f575eb 100644
--- a/stan/math/fwd/scal/fun/rising_factorial.hpp
+++ b/stan/math/fwd/scal/fun/rising_factorial.hpp
@@ -4,7 +4,6 @@
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/scal/fun/rising_factorial.hpp>
 #include <stan/math/prim/scal/fun/digamma.hpp>
-#include <iostream>
 
 namespace stan {
   namespace math {
diff --git a/stan/math/prim/mat/fun/cholesky_corr_constrain.hpp b/stan/math/prim/mat/fun/cholesky_corr_constrain.hpp
index 7ff8d90372d..bca942ba573 100644
--- a/stan/math/prim/mat/fun/cholesky_corr_constrain.hpp
+++ b/stan/math/prim/mat/fun/cholesky_corr_constrain.hpp
@@ -7,7 +7,6 @@
 #include <stan/math/prim/scal/fun/square.hpp>
 #include <stan/math/prim/scal/fun/corr_constrain.hpp>
 #include <cmath>
-#include <iostream>
 
 namespace stan {
   namespace math {
diff --git a/stan/math/prim/mat/fun/factor_U.hpp b/stan/math/prim/mat/fun/factor_U.hpp
index f4b9f5268d5..80fbcb649aa 100644
--- a/stan/math/prim/mat/fun/factor_U.hpp
+++ b/stan/math/prim/mat/fun/factor_U.hpp
@@ -6,10 +6,8 @@
 
 #include <cmath>
 #include <cstddef>
-#include <iostream>
 #include <limits>
 #include <stdexcept>
-#include <sstream>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/mat/fun/read_corr_L.hpp b/stan/math/prim/mat/fun/read_corr_L.hpp
index 2716bd5b022..6328e4cd7d2 100644
--- a/stan/math/prim/mat/fun/read_corr_L.hpp
+++ b/stan/math/prim/mat/fun/read_corr_L.hpp
@@ -6,7 +6,6 @@
 #include <stan/math/prim/scal/fun/square.hpp>
 #include <stan/math/prim/mat/fun/sum.hpp>
 #include <cstddef>
-#include <iostream>
 
 namespace stan {
   namespace math {
diff --git a/stan/math/prim/mat/fun/sort_indices.hpp b/stan/math/prim/mat/fun/sort_indices.hpp
index e9e83ed6917..658e83b0333 100644
--- a/stan/math/prim/mat/fun/sort_indices.hpp
+++ b/stan/math/prim/mat/fun/sort_indices.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/mat/meta/index_type.hpp>
 #include <stan/math/prim/arr/meta/index_type.hpp>
 #include <algorithm>
-#include <iostream>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/mat/fun/sort_indices_asc.hpp b/stan/math/prim/mat/fun/sort_indices_asc.hpp
index 3430a0a81b4..cbe64c17b2e 100644
--- a/stan/math/prim/mat/fun/sort_indices_asc.hpp
+++ b/stan/math/prim/mat/fun/sort_indices_asc.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/mat/meta/index_type.hpp>
 #include <stan/math/prim/mat/fun/sort_indices.hpp>
 #include <algorithm>    // std::sort
-#include <iostream>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/mat/fun/sort_indices_desc.hpp b/stan/math/prim/mat/fun/sort_indices_desc.hpp
index c9a594195e0..823d48db7c5 100644
--- a/stan/math/prim/mat/fun/sort_indices_desc.hpp
+++ b/stan/math/prim/mat/fun/sort_indices_desc.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/mat/meta/index_type.hpp>
 #include <stan/math/prim/mat/fun/sort_indices.hpp>
 #include <algorithm>    // std::sort
-#include <iostream>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/scal.hpp b/stan/math/prim/scal.hpp
index e241eada37b..9845930b871 100644
--- a/stan/math/prim/scal.hpp
+++ b/stan/math/prim/scal.hpp
@@ -37,6 +37,8 @@
 #include <stan/math/prim/scal/meta/VectorBuilder.hpp>
 #include <stan/math/prim/scal/meta/VectorView.hpp>
 
+#include <stan/math/prim/scal/err/check_2F1_converges.hpp>
+#include <stan/math/prim/scal/err/check_3F2_converges.hpp>
 #include <stan/math/prim/scal/err/check_bounded.hpp>
 #include <stan/math/prim/scal/err/check_consistent_size.hpp>
 #include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
diff --git a/stan/math/prim/scal/err/check_2F1_converges.hpp b/stan/math/prim/scal/err/check_2F1_converges.hpp
new file mode 100644
index 00000000000..fe180ba0101
--- /dev/null
+++ b/stan/math/prim/scal/err/check_2F1_converges.hpp
@@ -0,0 +1,70 @@
+#ifndef STAN_MATH_PRIM_SCAL_ERR_CHECK_2F1_CONVERGES_HPP
+#define STAN_MATH_PRIM_SCAL_ERR_CHECK_2F1_CONVERGES_HPP
+
+#include <stan/math/prim/scal/fun/value_of_rec.hpp>
+#include <cmath>
+#include <stdexcept>
+#include <sstream>
+#include <limits>
+
+namespace stan {
+  namespace math {
+
+    /**
+     * Check if the hypergeometric function (2F1) called with
+     * supplied arguments will converge, assuming arguments are
+     * finite values.
+     *
+     * @tparam T_a1 Type of a1
+     * @tparam T_a2 Type of a2
+     * @tparam T_b1 Type of b1
+     * @tparam T_z Type of z
+     *
+     * @param function Name of function ultimately relying on 2F1 (for error
+     &   messages)
+     * @param a1 Variable to check
+     * @param a2 Variable to check
+     * @param b1 Variable to check
+     * @param z Variable to check
+     *
+     * @throw <code>domain_error</code> if 2F1(a1, a2, b1, z)
+     *   does not meet convergence conditions, or if any coefficient is NaN.
+     */
+    template <typename T_a1, typename T_a2, typename T_b1, typename T_z>
+    inline void check_2F1_converges(const char* function,
+      const T_a1& a1, const T_a2& a2, const T_b1& b1, const T_z& z
+    ) {
+      using std::floor;
+      using std::fabs;
+
+      int num_terms = 0;
+      bool is_polynomial;
+      is_polynomial = (a1 <= 0.0 && floor(a1) == a1) ||
+                      (a2 <= 0.0 && floor(a2) == a2);
+      if (is_polynomial) {
+        if (a1 < 0.0 && floor(a1) == a1 && fabs(a1) > num_terms)
+          num_terms = floor(fabs(value_of_rec(a1)));
+        if (a2 < 0.0 && floor(a2) == a2 && fabs(a2) > num_terms)
+          num_terms = floor(fabs(value_of_rec(a2)));
+      }
+
+      bool is_undefined;
+      is_undefined = (b1 <= 0.0 && floor(b1) == b1 && fabs(b1) <= num_terms);
+
+      if (is_polynomial && !is_undefined) return;
+      if (fabs(z) < 1.0 && !is_undefined) return;
+      if (fabs(z) == 1.0 && !is_undefined) {
+        if (b1 > a1 + a2) return;
+      }
+      std::stringstream msg;
+      msg << "called from function '" << function << "', "
+          << "hypergeometric function 2F1 does not meet convergence "
+          << "conditions with given arguments. "
+          << "a1: " << a1 << ", a2: " << a2 << ", "
+          << "b1: " << b1 << ", z: " << z;
+      throw std::domain_error(msg.str());
+    }
+
+  }
+}
+#endif
diff --git a/stan/math/prim/scal/err/check_3F2_converges.hpp b/stan/math/prim/scal/err/check_3F2_converges.hpp
new file mode 100644
index 00000000000..6a476bfe9ed
--- /dev/null
+++ b/stan/math/prim/scal/err/check_3F2_converges.hpp
@@ -0,0 +1,80 @@
+#ifndef STAN_MATH_PRIM_SCAL_ERR_CHECK_3F2_CONVERGES_HPP
+#define STAN_MATH_PRIM_SCAL_ERR_CHECK_3F2_CONVERGES_HPP
+
+#include <stan/math/prim/scal/fun/value_of_rec.hpp>
+#include <cmath>
+#include <stdexcept>
+#include <sstream>
+#include <limits>
+
+namespace stan {
+  namespace math {
+
+    /**
+     * Check if the hypergeometric function (3F2) called with
+     * supplied arguments will converge, assuming arguments are
+     * finite values.
+     *
+     * @tparam T_a1 Type of a1
+     * @tparam T_a2 Type of a2
+     * @tparam T_a3 Type of a3
+     * @tparam T_b1 Type of b1
+     * @tparam T_b2 Type of b2
+     * @tparam T_z Type of z
+     *
+     * @param function Name of function ultimately relying on 3F2 (for error
+     &   messages)
+     * @param a1 Variable to check
+     * @param a2 Variable to check
+     * @param a3 Variable to check
+     * @param b1 Variable to check
+     * @param b2 Variable to check
+     * @param z Variable to check
+     *
+     * @throw <code>domain_error</code> if 3F2(a1, a2, a3, b1, b2, z)
+     *   does not meet convergence conditions, or if any coefficient is NaN.
+     */
+    template <typename T_a1, typename T_a2, typename T_a3, typename T_b1,
+      typename T_b2, typename T_z>
+    inline void check_3F2_converges(const char* function,
+      const T_a1& a1, const T_a2& a2, const T_a3& a3, const T_b1& b1,
+      const T_b2& b2, const T_z& z
+    ) {
+      using std::floor;
+      using std::fabs;
+
+      int num_terms = 0;
+      bool is_polynomial;
+      is_polynomial = (a1 <= 0.0 && floor(a1) == a1) ||
+                      (a2 <= 0.0 && floor(a2) == a2) ||
+                      (a3 <= 0.0 && floor(a3) == a3);
+      if (is_polynomial) {
+        if (a1 <= 0.0 && floor(a1) == a1 && fabs(a1) > num_terms)
+          num_terms = floor(fabs(value_of_rec(a1)));
+        if (a2 <= 0.0 && floor(a2) == a2 && fabs(a2) > num_terms)
+          num_terms = floor(fabs(value_of_rec(a2)));
+        if (a3 <= 0.0 && floor(a3) == a3 && fabs(a3) > num_terms)
+          num_terms = floor(fabs(value_of_rec(a3)));
+      }
+
+      bool is_undefined;
+      is_undefined = (b1 <= 0.0 && floor(b1) == b1 && fabs(b1) <= num_terms) ||
+                     (b2 <= 0.0 && floor(b2) == b2 && fabs(b2) <= num_terms);
+
+      if (is_polynomial && !is_undefined) return;
+      if (fabs(z) < 1.0 && !is_undefined) return;
+      if (fabs(z) == 1.0 && !is_undefined) {
+        if (b1 + b2 > a1 + a2 + a3) return;
+      }
+      std::stringstream msg;
+      msg << "called from function '" << function << "', "
+          << "hypergeometric function 3F2 does not meet convergence "
+          << "conditions with given arguments. "
+          << "a1: " << a1 << ", a2: " << a2 << ", a3: " << a3
+          << ", b1: " << b1 << ", b2: " << b2 << ", z: " << z;
+      throw std::domain_error(msg.str());
+    }
+
+  }
+}
+#endif
diff --git a/stan/math/prim/scal/fun/F32.hpp b/stan/math/prim/scal/fun/F32.hpp
index bf8b696586d..94dec1b8b8c 100644
--- a/stan/math/prim/scal/fun/F32.hpp
+++ b/stan/math/prim/scal/fun/F32.hpp
@@ -1,42 +1,88 @@
 #ifndef STAN_MATH_PRIM_SCAL_FUN_F32_HPP
 #define STAN_MATH_PRIM_SCAL_FUN_F32_HPP
 
+#include <stan/math/prim/scal/fun/sign.hpp>
+#include <stan/math/prim/scal/err/domain_error.hpp>
+#include <stan/math/prim/scal/fun/is_nan.hpp>
+#include <stan/math/prim/scal/fun/is_inf.hpp>
+#include <stan/math/prim/scal/err/check_3F2_converges.hpp>
 #include <cmath>
 
 namespace stan {
   namespace math {
 
+    /**
+     * Hypergeometric function, 3F2.
+     *
+     * Calculate the hypergeometric function (3F2) as the power series
+     * directly to within <code>precision</code> or until
+     * <code>max_steps</code> terms.
+     *
+     * This function does not have a closed form but will converge if:
+     *   - <code>|z|</code> is less than 1
+     *   - <code>|z|</code> is equal to one and <code>b1 + b2 < a1 + a2 + a3</code>
+     * This function is a rational polynomial if
+     *   - <code>a1</code>, <code>a2</code>, or <code>a3</code> is a
+     *     non-positive integer
+     * This function can be treated as a rational polynomial if
+     *   - <code>b1</code> or <code>b2</code> is a non-positive integer
+     *     and the series is terminated prior to the final term.
+     *
+     * @tparam T type of arguments and result
+     * @param[in] a1 a1 (always called with 1 from beta binomial cdfs)
+     * @param[in] a2 a2 (always called with a2 > 1)
+     * @param[in] a3 a3 (always called with int a3 <= 0)
+     * @param[in] b1 b1 (always called with int b1 < |a3|)
+     * @param[in] b2 b2 (always <= 1)
+     * @param[in] z z (is always called with 1 from beta binomial cdfs)
+     * @param[in] precision precision of the infinite sum. defaults to 1e-6
+     * @param[in] max_steps number of steps to take. defaults to 10000
+     */
     template<typename T>
-    T F32(T a, T b, T c, T d, T e, T z, T precision = 1e-6) {
+    T F32(const T& a1, const T& a2, const T& a3, const T& b1, const T& b2,
+          const T& z, double precision = 1e-6, int max_steps = 1e5) {
+      check_3F2_converges("F32", a1, a2, a3, b1, b2, z);
+
       using std::exp;
       using std::log;
       using std::fabs;
+      using stan::math::is_nan;
 
       T F = 1.0;
-
       T tNew = 0.0;
-
       T logT = 0.0;
-
       T logZ = log(z);
 
-      int k = 0.0;
-
-      while (fabs(tNew) > precision || k == 0) {
-        T p = (a + k) * (b + k) * (c + k) / ( (d + k) * (e + k) * (k + 1) );
-
-        // If a, b, or c is a negative integer then the series terminates
-        // after a finite number of interations
-        if (p == 0) break;
-
-        logT +=  (p > 0 ? 1.0 : -1.0) * log(fabs(p)) + logZ;
-
-        tNew = exp(logT);
+      int k = 0;
+      bool T_is_negative = false;
+      T p = 0.0;
+      do {
+        p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1));
+        if (is_nan(p) || p == 0)
+          break;
 
+        logT += log(fabs(p)) + logZ;
+        if (p < 0 && T_is_negative) {
+          T_is_negative = false;
+        } else if (p < 0 && !T_is_negative) {
+          T_is_negative = true;
+        }
+        if (T_is_negative)
+          tNew = -1 * exp(logT);
+        else
+          tNew = exp(logT);
         F += tNew;
 
         ++k;
-      }
+        if (k >= max_steps) {
+          domain_error("F32", "k (internal counter)", max_steps, "exceeded ",
+            " iterations, hypergeometric function did not converge.");
+        }
+        if (is_inf(F)) {
+          domain_error("F32", "F (output)", F,
+            "overflow ", " hypergeometric function did not converge.");
+        }
+      } while (fabs(tNew) > precision);
       return F;
     }
 
diff --git a/stan/math/prim/scal/fun/grad_2F1.hpp b/stan/math/prim/scal/fun/grad_2F1.hpp
index 748a4d03b1d..bc09c1951e2 100644
--- a/stan/math/prim/scal/fun/grad_2F1.hpp
+++ b/stan/math/prim/scal/fun/grad_2F1.hpp
@@ -1,42 +1,121 @@
 #ifndef STAN_MATH_PRIM_SCAL_FUN_GRAD_2F1_HPP
 #define STAN_MATH_PRIM_SCAL_FUN_GRAD_2F1_HPP
 
+#include <stan/math/prim/scal/fun/sign.hpp>
+#include <stan/math/prim/scal/err/domain_error.hpp>
+#include <stan/math/prim/scal/fun/is_nan.hpp>
+#include <stan/math/prim/scal/err/check_2F1_converges.hpp>
 #include <cmath>
+#include <limits>
 
 namespace stan {
   namespace math {
 
-    // Gradient of the hypergeometric function 2F1(a, b | c | z)
-    // with respect to a and c
+    /**
+     * Gradients of the hypergeometric function, 2F1.
+     *
+     * Calculate the gradients of the hypergeometric function (2F1)
+     * as the power series stopping when the series converges
+     * to within <code>precision</code> or throwing when the
+     * function takes <code>max_steps</code> steps.
+     *
+     * This power-series representation converges for all gradients
+     * under the same conditions as the 2F1 function itself.
+     *
+     * @tparam T type of arguments and result
+     * @param[out] g g pointer to array of six values of type T, result.
+     * @param[in] a1 a1 see generalized hypergeometric function definition.
+     * @param[in] a2 a2 see generalized hypergeometric function definition.
+     * @param[in] b1 b1 see generalized hypergeometric function definition.
+     * @param[in] z z see generalized hypergeometric function definition.
+     * @param[in] precision precision of the infinite sum. defaults to 1e-6
+     * @param[in] max_steps number of steps to take. defaults to 10,000
+     */
     template<typename T>
-    void grad_2F1(T& gradA, T& gradC, T a, T b, T c, T z, T precision = 1e-6) {
+    void grad_2F1(T& g_a1, T& g_b1, const T& a1, const T& a2, const T& b1,
+        const T& z, const T& precision = 1e-10, int max_steps = 1e5) {
+      check_2F1_converges("grad_2F1", a1, a2, b1, z);
+
+      using std::log;
       using std::fabs;
+      using std::exp;
+      using stan::math::is_nan;
 
-      gradA = 0;
-      gradC = 0;
+      g_a1 = 0.0;
+      g_b1 = 0.0;
 
-      T gradAold = 0;
-      T gradCold = 0;
+      T log_g_old[4];
+      for (T *q = log_g_old; q != log_g_old + 4; ++q)
+        *q = -1.0 * std::numeric_limits<double>::infinity();
 
-      int k = 0;
-      T tDak = 1.0 / (a - 1);
+      T log_t_old = 0.0;
+      T log_t_new = 0.0;
+
+      T log_z = log(z);
+
+      T p = 0.0;
 
-      while (fabs(tDak * (a + (k - 1)) ) > precision || k == 0) {
-          const T r = ( (a + k) / (c + k) ) * ( (b + k) / (T)(k + 1) ) * z;
-          tDak = r * tDak * (a + (k - 1)) / (a + k);
+      double log_t_new_sign = 1.0;
+      double log_t_old_sign = 1.0;
+      double log_g_old_sign[4];
+      for (double *q = log_g_old_sign; q != log_g_old_sign + 4; ++q)
+        *q = 1.0;
 
-          if (r == 0) break;
+      int k = 0;
+      T term;
+      while (true) {
+        p = (a1 + k) * (a2 + k) / ((b1 + k) * (1 + k));
+        if (p == 0)
+          break;
+
+        log_t_new += log(fabs(p)) + log_z;
+        if (p < 0 && log_t_new_sign < 0.0) {
+          log_t_new_sign = 1.0;
+        } else if (p < 0 && log_t_new_sign > 0.0) {
+          log_t_new_sign = -1.0;
+        }
+//      g_old[0] = t_new * (g_old[0] / t_old + 1.0 / (a1 + k));
+        term = log_g_old_sign[0] * log_t_old_sign *
+          exp(log_g_old[0] - log_t_old) + 1/(a1 + k);
+        log_g_old[0] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[0] = log_t_new_sign;
+        else
+          log_g_old_sign[0] = -1.0 * log_t_new_sign;
 
-          gradAold = r * gradAold + tDak;
-          gradCold = r * gradCold - tDak * ((a + k) / (c + k));
+//      g_old[1] = t_new * (g_old[1] / t_old + 1.0 / (a2 + k));
+        term = log_g_old_sign[1] * log_t_old_sign *
+          exp(log_g_old[1] - log_t_old) + 1/(a2 + k);
+        log_g_old[1] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[1] = log_t_new_sign;
+        else
+          log_g_old_sign[1] = -1.0 * log_t_new_sign;
+//      g_old[3] = t_new * (g_old[3] / t_old - 1.0 / (b1 + k));
+        term = log_g_old_sign[2] * log_t_old_sign *
+          exp(log_g_old[2] - log_t_old) - 1/(b1 + k);
+        log_g_old[2] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[2] = log_t_new_sign;
+        else
+          log_g_old_sign[2] = -1.0 * log_t_new_sign;
 
-          gradA += gradAold;
-          gradC += gradCold;
+        g_a1 += log_g_old_sign[0] * exp(log_g_old[0]);
+        g_b1 += log_g_old_sign[2] * exp(log_g_old[2]);
 
-          ++k;
+        if (exp(log_t_new) <= precision)
+          break;  // implicit abs
 
-          if (k > 200) break;
+        if (k >= max_steps) {
+          domain_error("grad_2F1", "k (internal counter)", max_steps,
+            "exceeded ", " iterations, hypergeometric function gradient "
+            "did not converge.");
         }
+
+        log_t_old = log_t_new;
+        log_t_old_sign = log_t_new_sign;
+        ++k;
+      }
     }
 
   }
diff --git a/stan/math/prim/scal/fun/grad_F32.hpp b/stan/math/prim/scal/fun/grad_F32.hpp
index 2c6f37d3b3b..1498b7a5bfb 100644
--- a/stan/math/prim/scal/fun/grad_F32.hpp
+++ b/stan/math/prim/scal/fun/grad_F32.hpp
@@ -1,57 +1,148 @@
 #ifndef STAN_MATH_PRIM_SCAL_FUN_GRAD_F32_HPP
 #define STAN_MATH_PRIM_SCAL_FUN_GRAD_F32_HPP
 
+#include <stan/math/prim/scal/fun/sign.hpp>
+#include <stan/math/prim/scal/err/domain_error.hpp>
+#include <stan/math/prim/scal/fun/is_nan.hpp>
+#include <stan/math/prim/scal/err/check_3F2_converges.hpp>
 #include <cmath>
+#include <limits>
 
 namespace stan {
   namespace math {
 
+    /**
+     * Gradients of the hypergeometric function, 3F2.
+     *
+     * Calculate the gradients of the hypergeometric function (3F2)
+     * as the power series stopping when the series converges
+     * to within <code>precision</code> or throwing when the
+     * function takes <code>max_steps</code> steps.
+     *
+     * This power-series representation converges for all gradients
+     * under the same conditions as the 3F2 function itself.
+     *
+     * @tparam T type of arguments and result
+     * @param[out] g g pointer to array of six values of type T, result.
+     * @param[in] a1 a1 see generalized hypergeometric function definition.
+     * @param[in] a2 a2 see generalized hypergeometric function definition.
+     * @param[in] a3 a3 see generalized hypergeometric function definition.
+     * @param[in] b1 b1 see generalized hypergeometric function definition.
+     * @param[in] b2 b2 see generalized hypergeometric function definition.
+     * @param[in] z z see generalized hypergeometric function definition.
+     * @param[in] precision precision of the infinite sum. defaults to 1e-6
+     * @param[in] max_steps number of steps to take. defaults to 10,000
+     */
     template<typename T>
-    void grad_F32(T* g, T a, T b, T c, T d, T e, T z, T precision = 1e-6) {
+    void grad_F32(T* g, const T& a1, const T& a2, const T& a3, const T& b1,
+        const T& b2, const T& z, const T& precision = 1e-6,
+        int max_steps = 1e5) {
+      check_3F2_converges("grad_F32", a1, a2, a3, b1, b2, z);
+
       using std::log;
       using std::fabs;
       using std::exp;
+      using stan::math::is_nan;
 
-      T gOld[6];
-
-      for (T *p = g; p != g + 6; ++p) *p = 0;
-      for (T *p = gOld; p != gOld + 6; ++p) *p = 0;
-
-      T tOld = 1;
-      T tNew = 0;
-
-      T logT = 0;
-
-      T logZ = log(z);
-
-      int k = 0;
-
-      while (fabs(tNew) > precision || k == 0) {
-        T C = (a + k) / (d + k);
-        C *= (b + k) / (e + k);
-        C *= (c + k) / (1 + k);
-
-        // If a, b, or c is a negative integer then the series terminates
-        // after a finite number of interations
-        if (C == 0) break;
+      for (T *q = g; q != g + 6; ++q) *q = 0.0;
 
-        logT += (C > 0 ? 1 : -1) * log(fabs(C)) + logZ;
+      T log_g_old[6];
+      for (T *q = log_g_old; q != log_g_old + 6; ++q)
+        *q = -1.0 * std::numeric_limits<double>::infinity();
 
-        tNew = exp(logT);
+      T log_t_old = 0.0;
+      T log_t_new = 0.0;
 
-        gOld[0] = tNew * (gOld[0] / tOld + 1.0 / (a + k));
-        gOld[1] = tNew * (gOld[1] / tOld + 1.0 / (b + k));
-        gOld[2] = tNew * (gOld[2] / tOld + 1.0 / (c + k));
+      T log_z = log(z);
 
-        gOld[3] = tNew * (gOld[3] / tOld - 1.0 / (d + k));
-        gOld[4] = tNew * (gOld[4] / tOld - 1.0 / (e + k));
+      T p = 0.0;
 
-        gOld[5] = tNew * (gOld[5] / tOld + 1.0 / z);
-
-        for (int i = 0; i < 6; ++i) g[i] += gOld[i];
-
-        tOld = tNew;
+      double log_t_new_sign = 1.0;
+      double log_t_old_sign = 1.0;
+      double log_g_old_sign[6];
+      for (double *q = log_g_old_sign; q != log_g_old_sign + 6; ++q)
+        *q = 1.0;
 
+      int k = 0;
+      T term;
+      while (true) {
+        p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (1 + k));
+        if (p == 0)
+          break;
+
+        log_t_new += log(fabs(p)) + log_z;
+        if (p < 0 && log_t_new_sign < 0.0) {
+          log_t_new_sign = 1.0;
+        } else if (p < 0 && log_t_new_sign > 0.0) {
+          log_t_new_sign = -1.0;
+        }
+//        g_old[0] = t_new * (g_old[0] / t_old + 1.0 / (a1 + k));
+        term = log_g_old_sign[0] * log_t_old_sign *
+          exp(log_g_old[0] - log_t_old) + 1/(a1 + k);
+        log_g_old[0] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[0] = log_t_new_sign;
+        else
+          log_g_old_sign[0] = -1.0 * log_t_new_sign;
+
+//        g_old[1] = t_new * (g_old[1] / t_old + 1.0 / (a2 + k));
+        term = log_g_old_sign[1] * log_t_old_sign *
+          exp(log_g_old[1] - log_t_old) + 1/(a2 + k);
+        log_g_old[1] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[1] = log_t_new_sign;
+        else
+          log_g_old_sign[1] = -1.0 * log_t_new_sign;
+//        g_old[2] = t_new * (g_old[2] / t_old + 1.0 / (a3 + k));
+        term = log_g_old_sign[2] * log_t_old_sign *
+          exp(log_g_old[2] - log_t_old) + 1/(a3 + k);
+        log_g_old[2] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[2] = log_t_new_sign;
+        else
+          log_g_old_sign[2] = -1.0 * log_t_new_sign;
+//
+//        g_old[3] = t_new * (g_old[3] / t_old - 1.0 / (b1 + k));
+        term = log_g_old_sign[3] * log_t_old_sign *
+          exp(log_g_old[3] - log_t_old) - 1/(b1 + k);
+        log_g_old[3] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[3] = log_t_new_sign;
+        else
+          log_g_old_sign[3] = -1.0 * log_t_new_sign;
+//        g_old[4] = t_new * (g_old[4] / t_old - 1.0 / (b2 + k));
+        term = log_g_old_sign[4] * log_t_old_sign *
+          exp(log_g_old[4] - log_t_old) - 1/(b2 + k);
+        log_g_old[4] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[4] = log_t_new_sign;
+        else
+          log_g_old_sign[4] = -1.0 * log_t_new_sign;
+//
+//        g_old[5] = t_new * (g_old[5] / t_old + 1.0 / z);
+        term = log_g_old_sign[5] * log_t_old_sign *
+          exp(log_g_old[5] - log_t_old) + 1/z;
+        log_g_old[5] = log_t_new + log(fabs(term));
+        if (term >= 0.0)
+          log_g_old_sign[5] = log_t_new_sign;
+        else
+          log_g_old_sign[5] = -1.0 * log_t_new_sign;
+
+        for (int i = 0; i < 6; ++i) {
+          g[i] += log_g_old_sign[i] * exp(log_g_old[i]);
+        }
+
+        if (exp(log_t_new) <= precision)
+          break;  // implicit abs
+
+        if (k >= max_steps) {
+          domain_error("grad_F32", "k (internal counter)", max_steps,
+            "exceeded ", " iterations, hypergeometric function gradient "
+            "did not converge.");
+        }
+
+        log_t_old = log_t_new;
+        log_t_old_sign = log_t_new_sign;
         ++k;
       }
     }
diff --git a/stan/math/prim/scal/fun/grad_reg_inc_gamma.hpp b/stan/math/prim/scal/fun/grad_reg_inc_gamma.hpp
index 15e549e56db..4f7741aa746 100644
--- a/stan/math/prim/scal/fun/grad_reg_inc_gamma.hpp
+++ b/stan/math/prim/scal/fun/grad_reg_inc_gamma.hpp
@@ -23,6 +23,9 @@ namespace stan {
      * @param g   boost::math::tgamma(a) (precomputed value)
      * @param dig boost::math::digamma(a) (precomputed value)
      * @param precision required precision; applies to series expansion only
+     * @param max_steps number of steps to take. defaults to 10000
+     * @throw throws std::domain_error if not converged after max_steps
+     *   or increment overflows to inf.
      *
      * For the asymptotic expansion, the gradient is given by:
        \f[
@@ -37,7 +40,8 @@ namespace stan {
        \f]
      */
     template<typename T>
-    T grad_reg_inc_gamma(T a, T z, T g, T dig, double precision = 1e-6) {
+    T grad_reg_inc_gamma(T a, T z, T g, T dig, double precision = 1e-6,
+        int max_steps = 1e5) {
       using std::domain_error;
       using std::exp;
       using std::fabs;
@@ -81,6 +85,11 @@ namespace stan {
           ++k;
           s *= - z / k;
           delta = s / square(k + a);
+          if (k >= max_steps)
+            stan::math::domain_error("grad_reg_inc_gamma",
+              "k (internal counter)",
+              max_steps, "exceeded ",
+              " iterations, gamma function gradient did not converge.");
           if (is_inf(delta))
             stan::math::domain_error("grad_reg_inc_gamma",
                                      "is not converging", "", "");
diff --git a/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp b/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
index 6a945e780b8..6f1495b75ee 100644
--- a/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
+++ b/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
@@ -103,10 +103,9 @@ namespace stan {
         const T_partials_return mu = alpha_dbl + n_dbl + 1;
         const T_partials_return nu = beta_dbl + N_dbl - n_dbl - 1;
 
-        const T_partials_return F = F32((T_partials_return)1, mu,
-                                                    -N_dbl + n_dbl + 1,
-                                                    n_dbl + 2, 1 - nu,
-                                                    (T_partials_return)1);
+        T_partials_return F;
+        F = F32((T_partials_return)1, mu, -N_dbl + n_dbl + 1, n_dbl + 2,
+          1 - nu, (T_partials_return)1);
 
         T_partials_return C = lgamma(nu) - lgamma(N_dbl - n_dbl);
         C += lgamma(mu) - lgamma(n_dbl + 2);
diff --git a/stan/math/prim/scal/prob/beta_binomial_rng.hpp b/stan/math/prim/scal/prob/beta_binomial_rng.hpp
index 7dbe0342006..5958d8e5125 100644
--- a/stan/math/prim/scal/prob/beta_binomial_rng.hpp
+++ b/stan/math/prim/scal/prob/beta_binomial_rng.hpp
@@ -13,8 +13,6 @@
 #include <stan/math/prim/scal/meta/include_summand.hpp>
 #include <stan/math/prim/scal/prob/binomial_rng.hpp>
 #include <stan/math/prim/scal/prob/beta_rng.hpp>
-#include <stan/math/prim/scal/fun/F32.hpp>
-#include <stan/math/prim/scal/fun/grad_F32.hpp>
 
 namespace stan {
   namespace math {
diff --git a/stan/math/rev/scal/fun/grad_inc_beta.hpp b/stan/math/rev/scal/fun/grad_inc_beta.hpp
index eadcb3d5e84..e4026094d8f 100644
--- a/stan/math/rev/scal/fun/grad_inc_beta.hpp
+++ b/stan/math/rev/scal/fun/grad_inc_beta.hpp
@@ -1,17 +1,18 @@
 #ifndef STAN_MATH_REV_SCAL_FUN_GRAD_INC_BETA_HPP
 #define STAN_MATH_REV_SCAL_FUN_GRAD_INC_BETA_HPP
 
-#include <stan/math/prim/scal/fun/grad_2F1.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/prim/scal/fun/lbeta.hpp>
+#include <stan/math/prim/scal/fun/grad_2F1.hpp>
+#include <stan/math/prim/scal/fun/value_of.hpp>
 #include <stan/math/rev/scal/fun/exp.hpp>
 #include <stan/math/rev/scal/fun/fabs.hpp>
-#include <stan/math/prim/scal/fun/lbeta.hpp>
+#include <stan/math/rev/scal/fun/inc_beta.hpp>
+#include <stan/math/rev/scal/fun/is_nan.hpp>
 #include <stan/math/rev/scal/fun/lgamma.hpp>
 #include <stan/math/rev/scal/fun/log.hpp>
 #include <stan/math/rev/scal/fun/log1m.hpp>
 #include <stan/math/rev/scal/fun/value_of.hpp>
-#include <stan/math/prim/scal/fun/value_of.hpp>
-#include <stan/math/rev/scal/fun/inc_beta.hpp>
 #include <cmath>
 
 namespace stan {
diff --git a/test/unit/math/mix/scal/fun/grad_2F1_test.cpp b/test/unit/math/mix/scal/fun/grad_2F1_test.cpp
index b8955d3aaf1..a3762b18151 100644
--- a/test/unit/math/mix/scal/fun/grad_2F1_test.cpp
+++ b/test/unit/math/mix/scal/fun/grad_2F1_test.cpp
@@ -14,9 +14,9 @@ TEST(ProbInternalMath, grad2F1_fd1) {
   fvar<double> gradC;
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
 
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242,gradA.val_,1e-6);
-  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_,1e-5);
-  EXPECT_NEAR(-0.461773435230326182245722531773361592054302268779753796048,gradC.val_,1e-6);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242,gradA.val_,1e-8);
+  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_,1e-8);
+  EXPECT_NEAR(-0.461773435230326182245722531773361592054302268779753796048,gradC.val_,1e-8);
 }
 TEST(ProbInternalMath, grad2F1_fd2) {
   using stan::math::fvar;
@@ -30,9 +30,9 @@ TEST(ProbInternalMath, grad2F1_fd2) {
   fvar<double> gradC;
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
 
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242,gradA.val_,1e-6);
-  EXPECT_NEAR(-0.461773435230326182245722531773361592054302268779753796048,gradC.val_,1e-6);
-  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_,1e-5);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242,gradA.val_,1e-8);
+  EXPECT_NEAR(-0.461773435230326182245722531773361592054302268779753796048,gradC.val_,1e-8);
+  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_,1e-8);
 }
 TEST(ProbInternalMath, grad2F1_ffd1) {
   using stan::math::fvar;
@@ -46,9 +46,9 @@ TEST(ProbInternalMath, grad2F1_ffd1) {
   fvar<fvar<double> > gradC;
 
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928,gradA.val_.val_, 1e-6);
-  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_.val_, 1e-5);
-  EXPECT_NEAR(-0.46177343523032618224572253177336159205430226877975379604859,gradC.val_.val_, 1e-6);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928,gradA.val_.val_, 1e-8);
+  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_.val_, 1e-8);
+  EXPECT_NEAR(-0.46177343523032618224572253177336159205430226877975379604859,gradC.val_.val_, 1e-8);
 }
 TEST(ProbInternalMath, grad2F1_ffd2) {
   using stan::math::fvar;
@@ -62,9 +62,9 @@ TEST(ProbInternalMath, grad2F1_ffd2) {
   fvar<fvar<double> > gradC;
 
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928,gradA.val_.val_, 1e-6);
-  EXPECT_NEAR(-0.46177343523032618224572253177336159205430226877975379604859,gradC.val_.val_, 1e-6);
-  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_.val_, 1e-5);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928,gradA.val_.val_, 1e-8);
+  EXPECT_NEAR(-0.46177343523032618224572253177336159205430226877975379604859,gradC.val_.val_, 1e-8);
+  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_.val_, 1e-8);
 }
 
 TEST(ProbInternalMath, grad2F1_fv1) {
@@ -79,9 +79,9 @@ TEST(ProbInternalMath, grad2F1_fv1) {
   fvar<var> gradA; fvar<var> gradC;
 
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928, gradA.val_.val(),1e-6);
-  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_.val(), 1e-5);
-  EXPECT_NEAR(-0.4617734352303261822457225317733615920543022687797537960, gradC.val_.val(),1e-6);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928, gradA.val_.val(),1e-8);
+  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,gradA.d_.val(), 1e-8);
+  EXPECT_NEAR(-0.4617734352303261822457225317733615920543022687797537960, gradC.val_.val(),1e-8);
 }
 TEST(ProbInternalMath, grad2F1_fv2) {
   using stan::math::fvar;
@@ -95,9 +95,9 @@ TEST(ProbInternalMath, grad2F1_fv2) {
   fvar<var> gradA; fvar<var> gradC;
 
   stan::math::grad_2F1(gradA,gradC,a, b, c, z);
-  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928, gradA.val_.val(),1e-6);
-  EXPECT_NEAR(-0.4617734352303261822457225317733615920543022687797537960, gradC.val_.val(),1e-6);
-  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_.val(), 1e-5);
+  EXPECT_NEAR(0.4617734315397201318453321291834046302225919173588625242928, gradA.val_.val(),1e-8);
+  EXPECT_NEAR(-0.4617734352303261822457225317733615920543022687797537960, gradC.val_.val(),1e-8);
+  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,gradC.d_.val(), 1e-8);
 }
 
 TEST(ProbInternalMath, grad2F1_fv_1stderiv1) {
@@ -116,7 +116,7 @@ TEST(ProbInternalMath, grad2F1_fv_1stderiv1) {
   AVEC y1 = createAVEC(a.val_);
   VEC grad1;
   gradA.val_.grad(y1,grad1);
-  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,grad1[0],1e-5);
+  EXPECT_NEAR(0.163714876516383746459968120418298168600425943651588679302872,grad1[0],1e-8);
 }
 TEST(ProbInternalMath, grad2F1_fv_1stderiv2) {
   using stan::math::fvar;
@@ -134,7 +134,7 @@ TEST(ProbInternalMath, grad2F1_fv_1stderiv2) {
   AVEC y1 = createAVEC(c.val_);
   VEC grad1;
   gradC.val_.grad(y1,grad1);
-  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,grad1[0],1e-5);
+  EXPECT_NEAR(0.5744063304437309685867184312646717864627845936245830896889,grad1[0],1e-8);
 }
 
 TEST(ProbInternalMath, grad2F1_fv_2ndderiv1) {
@@ -153,7 +153,7 @@ TEST(ProbInternalMath, grad2F1_fv_2ndderiv1) {
   AVEC y1 = createAVEC(a.val_);
   VEC grad1;
   gradA.d_.grad(y1,grad1);
-  EXPECT_NEAR(0.06425652761307923044917291721823961650191124494852382302571,grad1[0],1e-5);
+  EXPECT_NEAR(0.06425652761307923044917291721823961650191124494852382302571,grad1[0],1e-8);
 }
 
 TEST(ProbInternalMath, grad2F1_fv_2ndderiv2) {
@@ -172,5 +172,5 @@ TEST(ProbInternalMath, grad2F1_fv_2ndderiv2) {
   AVEC y1 = createAVEC(c.val_);
   VEC grad1;
   gradC.d_.grad(y1,grad1);
-  EXPECT_NEAR(-1.000245537254470801442530432195413371212643855413220347277769,grad1[0],1e-5);
+  EXPECT_NEAR(-1.000245537254470801442530432195413371212643855413220347277769,grad1[0],1e-8);
 }
diff --git a/test/unit/math/mix/scal/fun/inc_beta_test.cpp b/test/unit/math/mix/scal/fun/inc_beta_test.cpp
index 03056ef2746..6bff26caa84 100644
--- a/test/unit/math/mix/scal/fun/inc_beta_test.cpp
+++ b/test/unit/math/mix/scal/fun/inc_beta_test.cpp
@@ -193,7 +193,7 @@ TEST(ProbInternalMath, inc_beta_ffv_3rddderiv1) {
   VEC grad1;
   z1.d_.d_.grad(y1,grad1);
   EXPECT_NEAR(0.079976746033671442, 
-              grad1[0],1e-6);
+              grad1[0],1.05e-6);
 }
 TEST(ProbInternalMath, inc_beta_ffv_3rddderiv2) {
   using stan::math::fvar;
diff --git a/test/unit/math/prim/scal/err/check_2F1_converges_test.cpp b/test/unit/math/prim/scal/err/check_2F1_converges_test.cpp
new file mode 100644
index 00000000000..b2f06d1c03e
--- /dev/null
+++ b/test/unit/math/prim/scal/err/check_2F1_converges_test.cpp
@@ -0,0 +1,75 @@
+#include <stan/math/prim/scal.hpp>
+#include <gtest/gtest.h>
+
+using stan::math::check_2F1_converges;
+
+TEST(passesOnConvergentArgs,Check2F1Converges) {
+  const char* function = "check_2F1_converges";
+  double a1 = 1.0;
+  double a2 = 1.0;
+  double b1 = 5.0;
+  double z = 0.3;
+
+  // in radius of convergence for z, other args don't matter
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 5.0;
+  z = 1.0;  // still in radius of convergence, ok
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 1.0;  // now in radius of convergences, but b1 is too small.
+  z = 1.0; 
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+
+  a1 = 10.0;
+  a2 = 1.0;
+  b1 = 10.0;  // now in radius of convergences, but b1 is too small.
+  z = 1.0; 
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 5.0;
+  z = 1.3;  // outside of radius of convergence for current implementation. 
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 1.0;
+  z = 0.99999999999; // b1 is small, but z < 1 so we're ok.
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 1.0;
+  z = -0.999999999999; // checking negative z, this is fine
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 1.0;
+  z = std::numeric_limits<double>::infinity(); // limits of range?
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = 1.0;
+  z = -1.0 * std::numeric_limits<double>::infinity(); // limits of range?
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+  EXPECT_THROW(check_2F1_converges(function, a1, a2, b1, z), std::domain_error);
+   
+  a1 = 1.0;
+  a2 = 1.0;
+  b1 = std::numeric_limits<double>::infinity(); // should be ok, underflow to zero (?) 
+  z = 0.5;
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+  EXPECT_NO_THROW(check_2F1_converges(function, a1, a2, b1, z));
+
+}
+
+ 
diff --git a/test/unit/math/prim/scal/err/check_3F2_converges_test.cpp b/test/unit/math/prim/scal/err/check_3F2_converges_test.cpp
new file mode 100644
index 00000000000..67e165fbeb5
--- /dev/null
+++ b/test/unit/math/prim/scal/err/check_3F2_converges_test.cpp
@@ -0,0 +1,103 @@
+#include <stan/math/prim/scal.hpp>
+#include <gtest/gtest.h>
+
+using stan::math::check_3F2_converges;
+
+TEST(passesOnConvergentArgs,Check3F2Converges) {
+  const char* function = "check_3F2_converges";
+  double a1 = 1.0;
+  double a2 = 1.0;
+  double a3 = 1.0;
+  double b1 = 5.0;
+  double b2 = 5.0;
+  double z = 0.3;
+
+  // in radius of convergence for z, other args don't matter
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 5.0;
+  b2 = 5.0;
+  z = 1.0;  // still in radius of convergence, ok
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 1.1;  // now in radius of convergences, but b1 is too small.
+  b2 = 1.1;  // now in radius of convergences, but b1 is too small.
+  z = 1.0; 
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+
+  a1 = 40.0; // a1 is too big 
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 10.0;  
+  b2 = 10.0;  
+  z = 1.0; 
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+
+  a1 = 5.0; 
+  a2 = 0.0;
+  a3 = 1.0;
+  b1 = 10.0;  
+  b2 = 10.0;  
+  z = 1.0; 
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 5.0;
+  b2 = 5.0;
+  z = 1.3;  // outside of radius of convergence for current implementation. 
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 1.0;
+  b2 = 1.0;
+  z = 0.99999999999; // b1 is small, but z < 1 so we're ok.
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 1.0;
+  b2 = 1.0;
+  z = -0.999999999999; // checking negative z, this is fine
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 10.0;
+  b2 = 10.0;
+  z = std::numeric_limits<double>::infinity(); // limits of range?
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = 10.0;
+  b2 = 10.0;
+  z = -1.0 * std::numeric_limits<double>::infinity(); // limits of range?
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+  EXPECT_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z), std::domain_error);
+   
+  a1 = 1.0;
+  a2 = 1.0;
+  a3 = 1.0;
+  b1 = std::numeric_limits<double>::infinity(); // should be ok, underflow to zero (?) 
+  b2 = std::numeric_limits<double>::infinity(); // should be ok, underflow to zero (?) 
+  z = 0.5;
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+  EXPECT_NO_THROW(check_3F2_converges(function, a1, a2, a3, b1, b2, z));
+
+}
+
+ 
diff --git a/test/unit/math/prim/scal/fun/F32_test.cpp b/test/unit/math/prim/scal/fun/F32_test.cpp
new file mode 100644
index 00000000000..b42adca2145
--- /dev/null
+++ b/test/unit/math/prim/scal/fun/F32_test.cpp
@@ -0,0 +1,51 @@
+#include <gtest/gtest.h>
+#include <stan/math/prim/scal/fun/F32.hpp>
+
+TEST(MathPrimScalFun, F32_converges_by_z) { // converge
+  EXPECT_NEAR(2.5,
+              stan::math::F32(1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10),
+              1e-8);
+}
+
+TEST(MathPrimScalFun, F32_polynomial) { // terminate by zero numerator, no sign-flip
+  EXPECT_NEAR(11.28855722705942,
+              stan::math::F32(1.0, 31.0, -27.0, 19.0, -41.0, .99999, 1e-10),
+              1e-8);
+}
+
+TEST(MathPrimScalFun, F32_short_polynomial) { // terminate by zero numerator, single-step, no sign-flip
+  EXPECT_NEAR(-0.08000000000000007,
+              stan::math::F32(1.0, 12.0, -1.0, 10.0, 1.0, .9),
+              1e-8);
+}
+
+TEST(MathPrimScalFun, F32_short_polynomial_undef) { // at pole, should throw 
+  EXPECT_THROW(
+              stan::math::F32(1.0, 12.0, -1.0, 10.0, -1.0, 1.0),
+              std::domain_error);
+}
+
+TEST(MathPrimScalFun, F32_sign_flip_numerator) { // converge, single sign flip via numerator
+  EXPECT_NEAR(0.96935324630667443905,
+              stan::math::F32(1.0, -.5, 2.0, 10.0, 1.0, 0.3, 1e-10),
+              1e-8);
+}
+
+TEST(MathPrimScalFun, F32_diverge_by_z) {
+  // This should throw (Mathematica claims the answer is -10 but... ?
+  EXPECT_THROW(
+              stan::math::F32(1.0, 12.0, 1.0, 10.0, 1.0, 1.1),
+              std::domain_error);
+}
+
+TEST(MathPrimScalFun, F32_double_sign_flip) { // convergence, double sign flip 
+  EXPECT_NEAR(1.03711889198028226149,
+              stan::math::F32(1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10),
+              1e-8);
+  EXPECT_NEAR(1.06593846110441323674,
+              stan::math::F32(1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10),
+              1e-8);
+}
+
+
+
diff --git a/test/unit/math/prim/scal/fun/grad_2F1_test.cpp b/test/unit/math/prim/scal/fun/grad_2F1_test.cpp
new file mode 100644
index 00000000000..3ffe08028e5
--- /dev/null
+++ b/test/unit/math/prim/scal/fun/grad_2F1_test.cpp
@@ -0,0 +1,136 @@
+#include <gtest/gtest.h>
+#include <vector>
+#include <stan/math/prim/scal.hpp>
+
+TEST(MathPrimScalFun, grad2F1_1) {
+  double a1 = 1.0;
+  double a2 = 1.0;
+  double b1 = 1.0;
+  double z = 0.6;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR( 2.290726829685388, grad_a1, 1e-8);
+  EXPECT_NEAR(-2.290726829685388, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_2) {
+  double a1 = 1.0;
+  double a2 = 31.0;
+  double b1 = 41.0;
+  double z = 1.0;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z, 1e-11);
+  EXPECT_NEAR( 6.825270649241036, grad_a1, 1e-8);
+  EXPECT_NEAR(-0.382716049382716, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_3) {
+  double a1 = 1.0;
+  double a2 = -2.1;
+  double b1 = 41.0;
+  double z = 1.0;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR(-0.04921317604093563, grad_a1, 1e-8);
+  EXPECT_NEAR( 0.00118482743834665, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_4) {
+  double a1 = 1.0;
+  double a2 = 12.0;
+  double b1 = 10.0;
+  double z = 1.0;
+
+  double grad_a1;
+  double grad_b1;
+  EXPECT_THROW(
+    stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z),
+    std::domain_error);
+}
+
+TEST(MathPrimScalFun, grad2F1_5) {
+  double a1 = 1.0;
+  double a2 = 12.0;
+  double b1 = 20.0;
+  double z = 1.2;
+
+  double grad_a1;
+  double grad_b1;
+  EXPECT_THROW(
+    stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z),
+    std::domain_error);
+}
+
+TEST(MathPrimScalFun, grad2F1_6) {
+  double a1 = 1.0;
+  double a2 = -0.5;
+  double b1 = 10.6;
+  double z = 0.3;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR(-0.01443822031245647, grad_a1, 1e-8);
+  EXPECT_NEAR( 0.00136986255602642, grad_b1, 1e-8);
+}
+
+
+TEST(MathPrimScalFun, grad2F1_7) {
+  double a1 = 1.0;  
+  double a2 = -0.5;
+  double b1 = 10.0;
+  double z = 0.3;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR(-0.0153218866216130, grad_a1, 1e-8);
+  EXPECT_NEAR( 0.0015413242328729, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_8) {
+  double a1 = -.5;  
+  double a2 = -4.5;
+  double b1 = 11.0;
+  double z = 0.3;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR(-0.1227022810085707, grad_a1, 1e-8);
+  EXPECT_NEAR(-0.0053540982315572, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_9) {
+  double a1 = -.5;  
+  double a2 = -4.5;
+  double b1 = -3.2;
+  double z = 0.9;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR( 0.85880025358111, grad_a1, 1e-8);
+  EXPECT_NEAR(-4.19010422485256, grad_b1, 1e-8);
+}
+
+TEST(MathPrimScalFun, grad2F1_10) {
+  double a1 = 2;
+  double a2 = 1.0;
+  double b1 = 2;
+  double z = 0.4;
+
+  double grad_a1;
+  double grad_b1;
+  stan::math::grad_2F1(grad_a1, grad_b1, a1, a2, b1, z);
+  EXPECT_NEAR( 0.4617734323582945, grad_a1, 1e-8);
+  EXPECT_NEAR(-0.4617734323582945, grad_b1, 1e-8);
+}
+
+
diff --git a/test/unit/math/prim/scal/fun/grad_F32_test.cpp b/test/unit/math/prim/scal/fun/grad_F32_test.cpp
new file mode 100644
index 00000000000..3822ccd8fe1
--- /dev/null
+++ b/test/unit/math/prim/scal/fun/grad_F32_test.cpp
@@ -0,0 +1,123 @@
+#include <gtest/gtest.h>
+#include <vector>
+#include <stan/math/prim/scal/fun/grad_F32.hpp>
+
+TEST(MathPrimScalFun, grad_F32_converges_by_z) { // converge
+  std::vector<double> g(6);
+  g[0] =  2.290726829685388;
+  g[1] =  2.290726829685388;
+  g[2] =  2.290726829685388;
+  g[3] = -2.290726829685388;
+  g[4] = -2.290726829685388;
+  g[5] =  6.249999999999999;
+  double g_calc[6];
+  stan::math::grad_F32(g_calc, 1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10);
+  for (int i=0; i < 6; ++i) 
+    EXPECT_NEAR(g[i], g_calc[i], 1e-8);
+}
+
+// FIXME: Can't run till we use a continuation, crosses pole in digamma
+//TEST(MathPrimScalFun, grad_F32_polynomial) { // terminate by zero numerator, no sign-flip
+//  EXPECT_NEAR(11.28855722705942,
+//              stan::math::grad_F32(1.0, 31.0, -27.0, 19.0, -41.0, .99999, 1e-10),
+//              1e-8);
+//}
+
+// FIXME: Can't run till we use a continuation, crosses pole in digamma
+//TEST(MathPrimScalFun, grad_F32_short_polynomial) { // terminate by zero numerator, single-step, no sign-flip
+//  std::vector<double> g(6);
+//  g[0] = -1.08;
+//  g[1] = -0.09;
+//  g[2] = -0.07784317107588729;
+//  g[3] =  0.108;
+//  g[4] =  1.08;
+//  g[5] = -1.2;
+//  double g_calc[6];
+//  stan::math::grad_F32(g_calc, 1.0, 12.0, -1.0, 10.0, 1.0, .9, 1e-20);
+//  for (int i=0; i < 6; ++i) {
+//    std::cout << i << std::endl;
+//    EXPECT_NEAR(g[i], g_calc[i], 1e-8);
+//  }
+//         
+//}
+
+TEST(MathPrimScalFun, grad_F32_short_polynomial_undef) { // at pole, should throw 
+  double g_calc[6];
+  EXPECT_THROW(
+    stan::math::grad_F32(g_calc, 1.0, 12.0, -1.0, 10.0, -1.0, 1.0),
+    std::domain_error);
+}
+
+// FIXME: Can't run until we use a continuation, crosses pole in digamma
+//TEST(MathPrimScalFun, grad_F32_sign_flip_numerator) { // converge, single sign flip via numerator
+//  std::vector<double> g(6);
+//  g[0] = -1.08;
+//  g[1] = -0.09;
+//  g[2] = -0.07784317107588729;
+//  g[3] =  0.108;
+//  g[4] =  1.08;
+//  g[5] = -1.2;
+//  double g_calc[6];
+//  stan::math::grad_F32(g_calc, 1.0, -.5, 2.0, 10.0, 1.0, 0.3, 1e-10);
+//  for (int i=0; i < 6; ++i) 
+//    EXPECT_NEAR(g[i], g_calc[i], 1e-8);
+//         
+//}
+
+TEST(MathPrimScalFun, grad_F32_diverge_by_z) {
+  // This should throw (Mathematica claims the answer is -10 but it's by continuation
+  double g_calc[6];
+  EXPECT_THROW(
+    stan::math::grad_F32(g_calc, 1.0, 12.0, 1.0, 10.0, 1.0, 1.1),
+    std::domain_error);
+}
+
+TEST(MathPrimScalFun, grad_F32_double_sign_flip_1) { // convergence, double sign flip 
+  double g_calc[6];
+  std::vector<double> g(6);
+  g[0] =  0.03692914737912187;
+  g[1] = -0.0749983829996620;
+  g[2] = -0.0145982785371077;
+  g[3] = -0.00367744563081578;
+  g[4] = -0.0369291473791218;
+  g[5] =  0.1224673892640317;
+  stan::math::grad_F32(g_calc, 1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10);
+  for (int i=0; i < 6; ++i) 
+    EXPECT_NEAR(g[i], g_calc[i], 1e-8);
+}
+
+TEST(MathPrimScalFun, grad_F32_double_sign_flip_2) { // convergence, double sign flip 
+  double g_calc[6];
+  std::vector<double> g(6);
+  g[0] =  0.06517384196658483;
+  g[1] = -0.1349675922478992;
+  g[2] = -0.01422583581177123;
+  g[3] = -0.00645591019603955;
+  g[4] = -0.0651738419665848;
+  g[5] =  0.2147503542109778;
+  stan::math::grad_F32(g_calc, 1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10);
+  for (int i=0; i < 6; ++i) 
+    EXPECT_NEAR(g[i], g_calc[i], 1e-8);
+}
+
+//
+//m = {
+//  {1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10},
+//  {1.0, 31.0, -27.0, 19.0, -41.0, .99999, 1e-10},
+//  {1.0, 12.0, -1.0, 10.0, 1.0, .9},
+//  {1.0, 12.0, -1.0, 10.0, -1.0, 1.0},
+//  {1.0, -.5, 2.0, 10.0, 1.0, 0.3, 1e-10},
+//  {1.0, 12.0, 1.0, 10.0, 1.0, 1.1},
+//  {1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10},
+//  {1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10}
+//}
+//
+//{2.290726829685388,2.290726829685388,2.290726829685388,-2.290726829685388,-2.290726829685388,6.249999999999999},
+//IGNORE NO GRADIENT => {22.95673741899624,1.73991856855972,Indeterminate,-2.692768015657851,1.606385109522632,59.65366361215116},
+//{-1.08,-0.09,-0.07784317107588729,0.108,1.08,-1.2},
+//THROW => {1.2,0.1,Indeterminate,-0.12,1.2,1.2},
+//{-0.03098182751890103,0.05997687803893869,-0.01554776878589169,0.003126435989341982,0.03098182751890103,-0.1044310337234595},
+//THROW => {-34.02585092994046+31.41592653589793 I,-5.371146427869409+24.22439931425321 I,-34.02585092994046+31.41592653589793 I,7.371146427869409-24.22439931425321 I,34.02585092994046-31.41592653589793 I,0.},
+//{0.03692914737912187,-0.07499838299966203,-0.01459827853710772,-0.003677445630815789,-0.03692914737912187,0.1224673892640317},
+//{0.06517384196658483,-0.1349675922478992,-0.01422583581177123,-0.006455910196039558,-0.06517384196658483,0.2147503542109778}
+//
diff --git a/test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp b/test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp
index de4952f2b13..92e747be3f7 100644
--- a/test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp
+++ b/test/unit/math/prim/scal/prob/beta_binomial_cdf_log_test.cpp
@@ -1,5 +1,6 @@
 #include <stan/math/prim/scal.hpp>
 #include <gtest/gtest.h>
+#include <cmath>
 
 TEST(ProbBetaBinomial, cdf_log_matches_lcdf) {
   int n = 2;
@@ -7,8 +8,51 @@ TEST(ProbBetaBinomial, cdf_log_matches_lcdf) {
   double alpha = 1.1;
   double beta = 0.3;
 
-  EXPECT_FLOAT_EQ((stan::math::beta_binomial_lcdf(n, N, alpha, beta)),
-                  (stan::math::beta_binomial_cdf_log(n, N, alpha, beta)));
-  EXPECT_FLOAT_EQ((stan::math::beta_binomial_lcdf<int, int, double, double>(n, N, alpha, beta)),
-                  (stan::math::beta_binomial_cdf_log<int, int, double, double>(n, N, alpha, beta)));
+  EXPECT_NEAR(stan::math::beta_binomial_lcdf(n, N, alpha, beta),
+    stan::math::beta_binomial_cdf_log(n, N, alpha, beta), 1e-8);
+  EXPECT_NEAR((stan::math::beta_binomial_lcdf(n, N, alpha, beta)),
+    stan::math::beta_binomial_cdf_log(n, N, alpha, beta), 1e-8);
 }
+
+
+TEST(ProbBetaBinomial, lcdf_like_lcdf) {
+  int n = 10;
+  int N = 10;
+  double alpha = 3.0;
+  double beta = 1.0;
+
+  EXPECT_NEAR(0.0, stan::math::beta_binomial_lcdf(n, N, alpha, beta), 1e-8);
+  EXPECT_NEAR(0.0, std::exp(stan::math::beta_binomial_lcdf(0.0, N, alpha, beta)), 1e-8);
+
+}
+
+TEST(ProbBetaBinomial, lcdf_matches_lpmf) {
+  int n = 9;
+  int N = 10;
+  double alpha = 3.0;
+  double beta = 2.1;
+
+  double pmf_sum = 0.0;
+  for (int i = 0; i <= n; ++i) 
+    pmf_sum += std::exp(stan::math::beta_binomial_lpmf(i, N, alpha, beta));
+
+  EXPECT_NEAR(
+    pmf_sum,
+    std::exp(stan::math::beta_binomial_lcdf(n, N, alpha, beta)),
+    1e-8);
+}
+
+TEST(ProbBetaBinomial, lcdf_matches_mathematica) {
+  int n = 8;
+  int N = 10;
+  double alpha = 3.0;
+  double beta = 1.0;
+
+  //  EXPECT_NEAR(-0.5500463, (stan::math::beta_binomial_lcdf(n, N, alpha, beta)), 1e-8);
+  // FIXME: this point _should_ be defined for the beta_binomial_lcdf to be defined 
+  // over its full parameter range but the power-series is not defined.  Leaving the test
+  // in place with the current behavior.
+  EXPECT_THROW(stan::math::beta_binomial_lcdf(n, N, alpha, beta), std::domain_error);
+}
+
+