Stan Math Library  2.20.0
reverse mode automatic differentiation
grad_2F1.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_GRAD_2F1_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_GRAD_2F1_HPP
3 
8 #include <cmath>
9 #include <limits>
10 
11 namespace stan {
12 namespace math {
13 
35 template <typename T>
36 void grad_2F1(T& g_a1, T& g_b1, const T& a1, const T& a2, const T& b1,
37  const T& z, const T& precision = 1e-10, int max_steps = 1e5) {
38  check_2F1_converges("grad_2F1", a1, a2, b1, z);
39 
40  using std::exp;
41  using std::fabs;
42  using std::log;
43 
44  g_a1 = 0.0;
45  g_b1 = 0.0;
46 
47  T log_g_old[2];
48  for (auto& i : log_g_old)
49  i = NEGATIVE_INFTY;
50 
51  T log_t_old = 0.0;
52  T log_t_new = 0.0;
53 
54  T log_z = log(z);
55 
56  double log_t_new_sign = 1.0;
57  double log_t_old_sign = 1.0;
58  double log_g_old_sign[2];
59  for (double& x : log_g_old_sign)
60  x = 1.0;
61 
62  for (int k = 0; k <= max_steps; ++k) {
63  T p = (a1 + k) * (a2 + k) / ((b1 + k) * (1 + k));
64  if (p == 0)
65  return;
66 
67  log_t_new += log(fabs(p)) + log_z;
68  log_t_new_sign = p >= 0.0 ? log_t_new_sign : -log_t_new_sign;
69 
70  T term = log_g_old_sign[0] * log_t_old_sign * exp(log_g_old[0] - log_t_old)
71  + inv(a1 + k);
72  log_g_old[0] = log_t_new + log(fabs(term));
73  log_g_old_sign[0] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
74 
75  term = log_g_old_sign[1] * log_t_old_sign * exp(log_g_old[1] - log_t_old)
76  - inv(b1 + k);
77  log_g_old[1] = log_t_new + log(fabs(term));
78  log_g_old_sign[1] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
79 
80  g_a1 += log_g_old_sign[0] > 0 ? exp(log_g_old[0]) : -exp(log_g_old[0]);
81  g_b1 += log_g_old_sign[1] > 0 ? exp(log_g_old[1]) : -exp(log_g_old[1]);
82 
83  if (log_t_new <= log(precision))
84  return; // implicit abs
85 
86  log_t_old = log_t_new;
87  log_t_old_sign = log_t_new_sign;
88  }
89  domain_error("grad_2F1", "k (internal counter)", max_steps, "exceeded ",
90  " iterations, hypergeometric function gradient "
91  "did not converge.");
92  return;
93 }
94 
95 } // namespace math
96 } // namespace stan
97 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
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)
Gradients of the hypergeometric function, 2F1.
Definition: grad_2F1.hpp:36
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
const double NEGATIVE_INFTY
Negative infinity.
Definition: constants.hpp:53
void check_2F1_converges(const char *function, const T_a1 &a1, const T_a2 &a2, const T_b1 &b1, const T_z &z)
Check if the hypergeometric function (2F1) called with supplied arguments will converge, assuming arguments are finite values.
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

     [ Stan Home Page ] © 2011–2018, Stan Development Team.