Stan Math Library  2.20.0
reverse mode automatic differentiation
grad_F32.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_GRAD_F32_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_GRAD_F32_HPP
3 
9 #include <cmath>
10 #include <limits>
11 
12 namespace stan {
13 namespace math {
14 
37 template <typename T>
38 void grad_F32(T* g, const T& a1, const T& a2, const T& a3, const T& b1,
39  const T& b2, const T& z, const T& precision = 1e-6,
40  int max_steps = 1e5) {
41  check_3F2_converges("grad_F32", a1, a2, a3, b1, b2, z);
42 
43  using std::exp;
44  using std::fabs;
45  using std::log;
46 
47  for (int i = 0; i < 6; ++i)
48  g[i] = 0.0;
49 
50  T log_g_old[6];
51  for (auto& x : log_g_old)
52  x = NEGATIVE_INFTY;
53 
54  T log_t_old = 0.0;
55  T log_t_new = 0.0;
56 
57  T log_z = log(z);
58 
59  double log_t_new_sign = 1.0;
60  double log_t_old_sign = 1.0;
61  double log_g_old_sign[6];
62  for (int i = 0; i < 6; ++i)
63  log_g_old_sign[i] = 1.0;
64 
65  for (int k = 0; k <= max_steps; ++k) {
66  T p = (a1 + k) * (a2 + k) * (a3 + k) / ((b1 + k) * (b2 + k) * (1 + k));
67  if (p == 0)
68  return;
69 
70  log_t_new += log(fabs(p)) + log_z;
71  log_t_new_sign = p >= 0.0 ? log_t_new_sign : -log_t_new_sign;
72 
73  // g_old[0] = t_new * (g_old[0] / t_old + 1.0 / (a1 + k));
74  T term = log_g_old_sign[0] * log_t_old_sign * exp(log_g_old[0] - log_t_old)
75  + inv(a1 + k);
76  log_g_old[0] = log_t_new + log(fabs(term));
77  log_g_old_sign[0] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
78 
79  // g_old[1] = t_new * (g_old[1] / t_old + 1.0 / (a2 + k));
80  term = log_g_old_sign[1] * log_t_old_sign * exp(log_g_old[1] - log_t_old)
81  + inv(a2 + k);
82  log_g_old[1] = log_t_new + log(fabs(term));
83  log_g_old_sign[1] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
84 
85  // g_old[2] = t_new * (g_old[2] / t_old + 1.0 / (a3 + k));
86  term = log_g_old_sign[2] * log_t_old_sign * exp(log_g_old[2] - log_t_old)
87  + inv(a3 + k);
88  log_g_old[2] = log_t_new + log(fabs(term));
89  log_g_old_sign[2] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
90 
91  // g_old[3] = t_new * (g_old[3] / t_old - 1.0 / (b1 + k));
92  term = log_g_old_sign[3] * log_t_old_sign * exp(log_g_old[3] - log_t_old)
93  - inv(b1 + k);
94  log_g_old[3] = log_t_new + log(fabs(term));
95  log_g_old_sign[3] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
96 
97  // g_old[4] = t_new * (g_old[4] / t_old - 1.0 / (b2 + k));
98  term = log_g_old_sign[4] * log_t_old_sign * exp(log_g_old[4] - log_t_old)
99  - inv(b2 + k);
100  log_g_old[4] = log_t_new + log(fabs(term));
101  log_g_old_sign[4] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
102 
103  // g_old[5] = t_new * (g_old[5] / t_old + 1.0 / z);
104  term = log_g_old_sign[5] * log_t_old_sign * exp(log_g_old[5] - log_t_old)
105  + inv(z);
106  log_g_old[5] = log_t_new + log(fabs(term));
107  log_g_old_sign[5] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign;
108 
109  for (int i = 0; i < 6; ++i) {
110  g[i] += log_g_old_sign[i] * exp(log_g_old[i]);
111  }
112 
113  if (log_t_new <= log(precision))
114  return; // implicit abs
115 
116  log_t_old = log_t_new;
117  log_t_old_sign = log_t_new_sign;
118  }
119  domain_error("grad_F32", "k (internal counter)", max_steps, "exceeded ",
120  " iterations, hypergeometric function gradient "
121  "did not converge.");
122  return;
123 }
124 
125 } // namespace math
126 } // namespace stan
127 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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)
Check if the hypergeometric function (3F2) called with supplied arguments will converge, assuming arguments are finite values.
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
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
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)
Gradients of the hypergeometric function, 3F2.
Definition: grad_F32.hpp:38
const double NEGATIVE_INFTY
Negative infinity.
Definition: constants.hpp:53
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

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