1 #ifndef STAN_MATH_PRIM_SCAL_FUN_GRAD_REG_INC_GAMMA_HPP 2 #define STAN_MATH_PRIM_SCAL_FUN_GRAD_REG_INC_GAMMA_HPP 46 int max_steps = 1e5) {
52 if (
is_nan(a))
return std::numeric_limits<T>::quiet_NaN();
53 if (
is_nan(z))
return std::numeric_limits<T>::quiet_NaN();
54 if (
is_nan(g))
return std::numeric_limits<T>::quiet_NaN();
55 if (
is_nan(dig))
return std::numeric_limits<T>::quiet_NaN();
58 if (z >= a && z >= 8) {
61 T a_minus_one_minus_k = a - 1;
62 T fac = a_minus_one_minus_k;
65 T delta = dfac / zpow;
67 for (
int k = 1; k < 10; ++k) {
68 a_minus_one_minus_k -= 1;
73 dfac = a_minus_one_minus_k * dfac + fac;
74 fac *= a_minus_one_minus_k;
79 "is not converging",
"",
"");
82 return gamma_q(a, z) * (l - dig) +
exp(-z + (a - 1) * l) * S / g;
90 T log_delta = log_s - 2 *
log(a);
91 for (
int k = 1; k <= max_steps; ++k) {
92 S += s_sign >= 0.0 ?
exp(log_delta) : -
exp(log_delta);
93 log_s += log_z -
log(k);
95 log_delta = log_s - 2 *
log(k + a);
98 "is not converging",
"",
"");
99 if (log_delta <=
log(precision))
100 return gamma_p(a, z) * ( dig - l ) +
exp( a * l ) * S / g;
103 "k (internal counter)",
104 max_steps,
"exceeded ",
105 " iterations, gamma function gradient did not converge.");
106 return std::numeric_limits<T>::infinity();
fvar< T > fabs(const fvar< T > &x)
fvar< T > log(const fvar< T > &x)
fvar< T > exp(const fvar< T > &x)
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.
int is_inf(const fvar< T > &x)
Returns 1 if the input's value is infinite and 0 otherwise.
fvar< T > gamma_p(const fvar< T > &x1, const fvar< T > &x2)
T grad_reg_inc_gamma(T a, T z, T g, T dig, double precision=1e-6, int max_steps=1e5)
Gradient of the regularized incomplete gamma functions igamma(a, z)
double e()
Return the base of the natural logarithm.
int is_nan(const fvar< T > &x)
Returns 1 if the input's value is NaN and 0 otherwise.
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)