Bugfix/issue 487 hypergeometric functions

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 {

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 {

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 {

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 {

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 {

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 {

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 {

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>

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

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

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;
     }