STLdoc
STLdocumentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Macros | Functions
gamma.tcc File Reference
#include "special_function_util.h"

Macros

#define _GLIBCXX_TR1_GAMMA_TCC   1
 

Functions

namespace std _GLIBCXX_VISIBILITY (default)
 

Detailed Description

This is an internal header file, included by other library headers. Do not attempt to use it directly. {tr1/cmath}

Macro Definition Documentation

#define _GLIBCXX_TR1_GAMMA_TCC   1

Function Documentation

namespace std _GLIBCXX_VISIBILITY ( default  )

This returns Bernoulli numbers from a table or by summation for larger values.

Recursion is unstable.

Parameters
__nthe order n of the Bernoulli number.
Returns
The Bernoulli number of order n.

This returns Bernoulli number $B_n$.

Parameters
__nthe order n of the Bernoulli number.
Returns
The Bernoulli number of order n.

Return $log(\Gamma(x))$ by asymptotic expansion with Bernoulli number coefficients. This is like Sterling's approximation.

Parameters
__xThe argument of the log of the gamma function.
Returns
The logarithm of the gamma function.

Return $log(\Gamma(x))$ by the Lanczos method. This method dominates all others on the positive axis I think.

Parameters
__xThe argument of the log of the gamma function.
Returns
The logarithm of the gamma function.

Return $ log(|\Gamma(x)|) $. This will return values even for $ x < 0 $. To recover the sign of $ \Gamma(x) $ for any argument use __log_gamma_sign.

Parameters
__xThe argument of the log of the gamma function.
Returns
The logarithm of the gamma function.

Return the sign of $ \Gamma(x) $. At nonpositive integers zero is returned.

Parameters
__xThe argument of the gamma function.
Returns
The sign of the gamma function.

Return the logarithm of the binomial coefficient. The binomial coefficient is given by:

\[ \left( \right) = \frac{n!}{(n-k)! k!} \]

Parameters
__nThe first argument of the binomial coefficient.
__kThe second argument of the binomial coefficient.
Returns
The binomial coefficient.

Return the binomial coefficient. The binomial coefficient is given by:

\[ \left( \right) = \frac{n!}{(n-k)! k!} \]

Parameters
__nThe first argument of the binomial coefficient.
__kThe second argument of the binomial coefficient.
Returns
The binomial coefficient.

Return $ \Gamma(x) $.

Parameters
__xThe argument of the gamma function.
Returns
The gamma function.

Return the digamma function by series expansion. The digamma or $ \psi(x) $ function is defined by

\[ \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} \]

The series is given by:

\[ \psi(x) = -\gamma_E - \frac{1}{x} \sum_{k=1}^{\infty} \frac{x}{k(x + k)} \]

Return the digamma function for large argument. The digamma or $ \psi(x) $ function is defined by

\[ \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} \]

The asymptotic series is given by:

\[ \psi(x) = \ln(x) - \frac{1}{2x} - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}} \]

Return the digamma function. The digamma or $ \psi(x) $ function is defined by

\[ \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} \]

For negative argument the reflection formula is used:

\[ \psi(x) = \psi(1-x) - \pi \cot(\pi x) \]

Return the polygamma function $ \psi^{(n)}(x) $.

The polygamma function is related to the Hurwitz zeta function:

\[ \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x) \]

52 {
53 namespace tr1
54 {
55  // Implementation-space details.
56  namespace __detail
57  {
58  _GLIBCXX_BEGIN_NAMESPACE_VERSION
59 
69  template <typename _Tp>
70  _Tp
71  __bernoulli_series(unsigned int __n)
72  {
73 
74  static const _Tp __num[28] = {
75  _Tp(1UL), -_Tp(1UL) / _Tp(2UL),
76  _Tp(1UL) / _Tp(6UL), _Tp(0UL),
77  -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
78  _Tp(1UL) / _Tp(42UL), _Tp(0UL),
79  -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
80  _Tp(5UL) / _Tp(66UL), _Tp(0UL),
81  -_Tp(691UL) / _Tp(2730UL), _Tp(0UL),
82  _Tp(7UL) / _Tp(6UL), _Tp(0UL),
83  -_Tp(3617UL) / _Tp(510UL), _Tp(0UL),
84  _Tp(43867UL) / _Tp(798UL), _Tp(0UL),
85  -_Tp(174611) / _Tp(330UL), _Tp(0UL),
86  _Tp(854513UL) / _Tp(138UL), _Tp(0UL),
87  -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),
88  _Tp(8553103UL) / _Tp(6UL), _Tp(0UL)
89  };
90 
91  if (__n == 0)
92  return _Tp(1);
93 
94  if (__n == 1)
95  return -_Tp(1) / _Tp(2);
96 
97  // Take care of the rest of the odd ones.
98  if (__n % 2 == 1)
99  return _Tp(0);
100 
101  // Take care of some small evens that are painful for the series.
102  if (__n < 28)
103  return __num[__n];
104 
105 
106  _Tp __fact = _Tp(1);
107  if ((__n / 2) % 2 == 0)
108  __fact *= _Tp(-1);
109  for (unsigned int __k = 1; __k <= __n; ++__k)
110  __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
111  __fact *= _Tp(2);
112 
113  _Tp __sum = _Tp(0);
114  for (unsigned int __i = 1; __i < 1000; ++__i)
115  {
116  _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
117  if (__term < std::numeric_limits<_Tp>::epsilon())
118  break;
119  __sum += __term;
120  }
121 
122  return __fact * __sum;
123  }
124 
125 
132  template<typename _Tp>
133  inline _Tp
134  __bernoulli(int __n)
135  { return __bernoulli_series<_Tp>(__n); }
136 
137 
146  template<typename _Tp>
147  _Tp
148  __log_gamma_bernoulli(_Tp __x)
149  {
150  _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
151  + _Tp(0.5L) * std::log(_Tp(2)
152  * __numeric_constants<_Tp>::__pi());
153 
154  const _Tp __xx = __x * __x;
155  _Tp __help = _Tp(1) / __x;
156  for ( unsigned int __i = 1; __i < 20; ++__i )
157  {
158  const _Tp __2i = _Tp(2 * __i);
159  __help /= __2i * (__2i - _Tp(1)) * __xx;
160  __lg += __bernoulli<_Tp>(2 * __i) * __help;
161  }
162 
163  return __lg;
164  }
165 
166 
174  template<typename _Tp>
175  _Tp
176  __log_gamma_lanczos(_Tp __x)
177  {
178  const _Tp __xm1 = __x - _Tp(1);
179 
180  static const _Tp __lanczos_cheb_7[9] = {
181  _Tp( 0.99999999999980993227684700473478L),
182  _Tp( 676.520368121885098567009190444019L),
183  _Tp(-1259.13921672240287047156078755283L),
184  _Tp( 771.3234287776530788486528258894L),
185  _Tp(-176.61502916214059906584551354L),
186  _Tp( 12.507343278686904814458936853L),
187  _Tp(-0.13857109526572011689554707L),
188  _Tp( 9.984369578019570859563e-6L),
189  _Tp( 1.50563273514931155834e-7L)
190  };
191 
192  static const _Tp __LOGROOT2PI
193  = _Tp(0.9189385332046727417803297364056176L);
194 
195  _Tp __sum = __lanczos_cheb_7[0];
196  for(unsigned int __k = 1; __k < 9; ++__k)
197  __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
198 
199  const _Tp __term1 = (__xm1 + _Tp(0.5L))
200  * std::log((__xm1 + _Tp(7.5L))
201  / __numeric_constants<_Tp>::__euler());
202  const _Tp __term2 = __LOGROOT2PI + std::log(__sum);
203  const _Tp __result = __term1 + (__term2 - _Tp(7));
204 
205  return __result;
206  }
207 
208 
218  template<typename _Tp>
219  _Tp
220  __log_gamma(_Tp __x)
221  {
222  if (__x > _Tp(0.5L))
223  return __log_gamma_lanczos(__x);
224  else
225  {
226  const _Tp __sin_fact
227  = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));
228  if (__sin_fact == _Tp(0))
229  std::__throw_domain_error(__N("Argument is nonpositive integer "
230  "in __log_gamma"));
231  return __numeric_constants<_Tp>::__lnpi()
232  - std::log(__sin_fact)
233  - __log_gamma_lanczos(_Tp(1) - __x);
234  }
235  }
236 
237 
245  template<typename _Tp>
246  _Tp
247  __log_gamma_sign(_Tp __x)
248  {
249  if (__x > _Tp(0))
250  return _Tp(1);
251  else
252  {
253  const _Tp __sin_fact
254  = std::sin(__numeric_constants<_Tp>::__pi() * __x);
255  if (__sin_fact > _Tp(0))
256  return (1);
257  else if (__sin_fact < _Tp(0))
258  return -_Tp(1);
259  else
260  return _Tp(0);
261  }
262  }
263 
264 
276  template<typename _Tp>
277  _Tp
278  __log_bincoef(unsigned int __n, unsigned int __k)
279  {
280  // Max e exponent before overflow.
281  static const _Tp __max_bincoeff
282  = std::numeric_limits<_Tp>::max_exponent10
283  * std::log(_Tp(10)) - _Tp(1);
284 #if _GLIBCXX_USE_C99_MATH_TR1
285  _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n))
286  - std::tr1::lgamma(_Tp(1 + __k))
287  - std::tr1::lgamma(_Tp(1 + __n - __k));
288 #else
289  _Tp __coeff = __log_gamma(_Tp(1 + __n))
290  - __log_gamma(_Tp(1 + __k))
291  - __log_gamma(_Tp(1 + __n - __k));
292 #endif
293  }
294 
295 
307  template<typename _Tp>
308  _Tp
309  __bincoef(unsigned int __n, unsigned int __k)
310  {
311  // Max e exponent before overflow.
312  static const _Tp __max_bincoeff
313  = std::numeric_limits<_Tp>::max_exponent10
314  * std::log(_Tp(10)) - _Tp(1);
315 
316  const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
317  if (__log_coeff > __max_bincoeff)
318  return std::numeric_limits<_Tp>::quiet_NaN();
319  else
320  return std::exp(__log_coeff);
321  }
322 
323 
330  template<typename _Tp>
331  inline _Tp
332  __gamma(_Tp __x)
333  { return std::exp(__log_gamma(__x)); }
334 
335 
349  template<typename _Tp>
350  _Tp
351  __psi_series(_Tp __x)
352  {
353  _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
354  const unsigned int __max_iter = 100000;
355  for (unsigned int __k = 1; __k < __max_iter; ++__k)
356  {
357  const _Tp __term = __x / (__k * (__k + __x));
358  __sum += __term;
359  if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
360  break;
361  }
362  return __sum;
363  }
364 
365 
379  template<typename _Tp>
380  _Tp
381  __psi_asymp(_Tp __x)
382  {
383  _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
384  const _Tp __xx = __x * __x;
385  _Tp __xp = __xx;
386  const unsigned int __max_iter = 100;
387  for (unsigned int __k = 1; __k < __max_iter; ++__k)
388  {
389  const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
390  __sum -= __term;
391  if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
392  break;
393  __xp *= __xx;
394  }
395  return __sum;
396  }
397 
398 
410  template<typename _Tp>
411  _Tp
412  __psi(_Tp __x)
413  {
414  const int __n = static_cast<int>(__x + 0.5L);
415  const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
416  if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)
417  return std::numeric_limits<_Tp>::quiet_NaN();
418  else if (__x < _Tp(0))
419  {
420  const _Tp __pi = __numeric_constants<_Tp>::__pi();
421  return __psi(_Tp(1) - __x)
422  - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
423  }
424  else if (__x > _Tp(100))
425  return __psi_asymp(__x);
426  else
427  return __psi_series(__x);
428  }
429 
430 
439  template<typename _Tp>
440  _Tp
441  __psi(unsigned int __n, _Tp __x)
442  {
443  if (__x <= _Tp(0))
444  std::__throw_domain_error(__N("Argument out of range "
445  "in __psi"));
446  else if (__n == 0)
447  return __psi(__x);
448  else
449  {
450  const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);
451 #if _GLIBCXX_USE_C99_MATH_TR1
452  const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
453 #else
454  const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
455 #endif
456  _Tp __result = std::exp(__ln_nfact) * __hzeta;
457  if (__n % 2 == 1)
458  __result = -__result;
459  return __result;
460  }
461  }
462 
463  _GLIBCXX_END_NAMESPACE_VERSION
464  } // namespace std::tr1::__detail
465 }
466 }
Generated on Sat Jul 22 2017 14:42:36 for STLdoc by   doxygen 1.8.6