Automatic Differentiation
 
Loading...
Searching...
No Matches
trigamma.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_TRIGAMMA_HPP
2#define STAN_MATH_PRIM_FUN_TRIGAMMA_HPP
3
12#include <cmath>
13
14// Reference:
15// BE Schneider,
16// Algorithm AS 121:
17// Trigamma Function,
18// Applied Statistics,
19// Volume 27, Number 1, pages 97-99, 1978.
20
21namespace stan {
22namespace math {
23
36template <typename T>
37inline T trigamma_impl(const T& x) {
38 using std::floor;
39 using std::sin;
40
41 double small = 0.0001;
42 double large = 5.0;
43 T value;
44 T y;
45 T z;
46
47 // bernoulli numbers
48 double b2 = inv(6.0);
49 double b4 = -inv(30.0);
50 double b6 = inv(42.0);
51 double b8 = -inv(30.0);
52
53 // negative integers and zero return positive infinity
54 // see http://mathworld.wolfram.com/PolygammaFunction.html
55 if (x <= 0.0 && floor(x) == x) {
56 value = positive_infinity();
57 return value;
58 }
59
60 // negative non-integers: use the reflection formula
61 // see http://mathworld.wolfram.com/PolygammaFunction.html
62 if (x <= 0 && floor(x) != x) {
63 value = -trigamma_impl(-x + 1.0) + square(pi() / sin(-pi() * x));
64 return value;
65 }
66
67 // small value approximation if x <= small.
68 if (x <= small) {
69 return inv_square(x);
70 }
71
72 // use recurrence relation until x >= large
73 // see http://mathworld.wolfram.com/PolygammaFunction.html
74 z = x;
75 value = 0.0;
76 while (z < large) {
77 value += inv_square(z);
78 z += 1.0;
79 }
80
81 // asymptotic expansion as a Laurent series if x >= large
82 // see http://en.wikipedia.org/wiki/Trigamma_function
83 y = inv_square(z);
84 value += 0.5 * y + (1.0 + y * (b2 + y * (b4 + y * (b6 + y * b8)))) / z;
85
86 return value;
87}
88
120inline double trigamma(double u) { return trigamma_impl(u); }
121
129inline double trigamma(int u) { return trigamma(static_cast<double>(u)); }
130
143 template <typename T>
144 static inline auto fun(const T& x) {
145 return trigamma(x);
146 }
147};
148
159template <typename T,
161inline auto trigamma(const T& x) {
163}
164
165} // namespace math
166} // namespace stan
167
168#endif
require_not_t< is_nonscalar_prim_or_rev_kernel_expression< std::decay_t< T > > > require_not_nonscalar_prim_or_rev_kernel_expression_t
Require type does not satisfy is_nonscalar_prim_or_rev_kernel_expression.
fvar< T > sin(const fvar< T > &x)
Definition sin.hpp:14
static constexpr double positive_infinity()
Return positive infinity.
fvar< T > inv_square(const fvar< T > &x)
T trigamma_impl(const T &x)
Return the trigamma function applied to the argument.
Definition trigamma.hpp:37
fvar< T > trigamma(const fvar< T > &u)
Return the value of the trigamma function at the specified argument (i.e., the second derivative of t...
Definition trigamma.hpp:21
fvar< T > floor(const fvar< T > &x)
Definition floor.hpp:12
static constexpr double pi()
Return the value of pi.
Definition constants.hpp:36
fvar< T > inv(const fvar< T > &x)
Definition inv.hpp:12
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Base template class for vectorization of unary scalar functions defined by a template class F to a sc...
static auto fun(const T &x)
Return the trigamma() function applied to the argument.
Definition trigamma.hpp:144
Structure to wrap trigamma() so it can be vectorized.
Definition trigamma.hpp:134