Stan Math Library  2.20.0
reverse mode automatic differentiation
trigamma.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_TRIGAMMA_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_TRIGAMMA_HPP
3 
5 // Reference:
6 // BE Schneider,
7 // Algorithm AS 121:
8 // Trigamma Function,
9 // Applied Statistics,
10 // Volume 27, Number 1, pages 97-99, 1978.
11 
15 #include <cmath>
16 
17 namespace stan {
18 namespace math {
19 
32 template <typename T>
33 inline T trigamma_impl(const T& x) {
34  using std::floor;
35  using std::sin;
36 
37  double small = 0.0001;
38  double large = 5.0;
39  T value;
40  T y;
41  T z;
42 
43  // bernoulli numbers
44  double b2 = inv(6.0);
45  double b4 = -inv(30.0);
46  double b6 = inv(42.0);
47  double b8 = -inv(30.0);
48 
49  // negative integers and zero return postiive infinity
50  // see http://mathworld.wolfram.com/PolygammaFunction.html
51  if ((x <= 0.0) && (floor(x) == x)) {
52  value = positive_infinity();
53  return value;
54  }
55 
56  // negative non-integers: use the reflection formula
57  // see http://mathworld.wolfram.com/PolygammaFunction.html
58  if ((x <= 0) && (floor(x) != x)) {
59  value = -trigamma_impl(-x + 1.0) + square(pi() / sin(-pi() * x));
60  return value;
61  }
62 
63  // small value approximation if x <= small.
64  if (x <= small)
65  return inv_square(x);
66 
67  // use recurrence relation until x >= large
68  // see http://mathworld.wolfram.com/PolygammaFunction.html
69  z = x;
70  value = 0.0;
71  while (z < large) {
72  value += inv_square(z);
73  z += 1.0;
74  }
75 
76  // asymptotic expansion as a Laurent series if x >= large
77  // see http://en.wikipedia.org/wiki/Trigamma_function
78  y = inv_square(z);
79  value += 0.5 * y + (1.0 + y * (b2 + y * (b4 + y * (b6 + y * b8)))) / z;
80 
81  return value;
82 }
83 
115 inline double trigamma(double u) { return trigamma_impl(u); }
116 
124 inline double trigamma(int u) { return trigamma(static_cast<double>(u)); }
125 
126 } // namespace math
127 } // namespace stan
128 
129 #endif
T trigamma_impl(const T &x)
Return the trigamma function applied to the argument.
Definition: trigamma.hpp:33
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:11
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:20
fvar< T > inv_square(const fvar< T > &x)
Definition: inv_square.hpp:12
double positive_infinity()
Return positive infinity.
Definition: constants.hpp:108
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:12
double pi()
Return the value of pi.
Definition: constants.hpp:80
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

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