1 #ifndef STAN_MATH_PRIM_FUN_F32_HPP
2 #define STAN_MATH_PRIM_FUN_F32_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/exp.hpp>
7 #include <stan/math/prim/fun/fabs.hpp>
8 #include <stan/math/prim/fun/is_inf.hpp>
9 #include <stan/math/prim/fun/log.hpp>
10 #include <cmath>
11
12 namespace stan {
13 namespace math {
14
15 /**
16 * Hypergeometric function (3F2).
17 *
18 * Function reference: http://dlmf.nist.gov/16.2
19 *
20 * \f[
21 * _3F_2 \left(
22 * \begin{matrix}a_1 a_2 a3 \\ b_1 b_2\end{matrix}; z
23 * \right) = \sum_k=0^\infty
24 * \frac{(a_1)_k(a_2)_k(a_3)_k}{(b_1)_k(b_2)_k}\frac{z^k}{k!} \f]
25 *
26 * Where $(a_1)_k$ is an upper shifted factorial.
27 *
28 * Calculate the hypergeometric function (3F2) as the power series
29 * directly to within <code>precision</code> or until
30 * <code>max_steps</code> terms.
31 *
32 * This function does not have a closed form but will converge if:
33 * - <code>|z|</code> is less than 1
34 * - <code>|z|</code> is equal to one and <code>b1 + b2 < a1 + a2 + a3</code>
35 * This function is a rational polynomial if
36 * - <code>a1</code>, <code>a2</code>, or <code>a3</code> is a
37 * non-positive integer
38 * This function can be treated as a rational polynomial if
39 * - <code>b1</code> or <code>b2</code> is a non-positive integer
40 * and the series is terminated prior to the final term.
41 *
42 * @tparam T type of arguments and result
43 * @param[in] a1 a1 (always called with 1 from beta binomial cdfs)
44 * @param[in] a2 a2 (always called with a2 > 1)
45 * @param[in] a3 a3 (always called with int a3 <= 0)
46 * @param[in] b1 b1 (always called with int b1 < |a3|)
47 * @param[in] b2 b2 (always <= 1)
48 * @param[in] z z (is always called with 1 from beta binomial cdfs)
49 * @param[in] precision precision of the infinite sum. defaults to 1e-6
50 * @param[in] max_steps number of steps to take. defaults to 1e5
51 */
52 template <typename 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)53 T F32(const T& a1, const T& a2, const T& a3, const T& b1, const T& b2,
54 const T& z, double precision = 1e-6, int max_steps = 1e5) {
55 check_3F2_converges("F32", a1, a2, a3, b1, b2, z);
56
57 using std::exp;
58 using std::fabs;
59 using std::log;
60
61 T t_acc = 1.0;
62 T log_t = 0.0;
63 T log_z = log(z);
64 double t_sign = 1.0;
65
66 for (int k = 0; k <= max_steps; ++k) {
67 T p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (k + 1));
68 if (p == 0.0) {
69 return t_acc;
70 }
71
72 log_t += log(fabs(p)) + log_z;
73 t_sign = p >= 0.0 ? t_sign : -t_sign;
74 T t_new = t_sign > 0.0 ? exp(log_t) : -exp(log_t);
75 t_acc += t_new;
76
77 if (fabs(t_new) <= precision) {
78 return t_acc;
79 }
80
81 if (is_inf(t_acc)) {
82 throw_domain_error("F32", "sum (output)", t_acc, "overflow ",
83 " hypergeometric function did not converge.");
84 }
85 }
86 throw_domain_error("F32", "k (internal counter)", max_steps, "exceeded ",
87 " iterations, hypergeometric function did not converge.");
88 return t_acc; // to silence warning.
89 }
90
91 } // namespace math
92 } // namespace stan
93 #endif
94