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

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