libstdc++
gamma.tcc
1 // Special functions -*- C++ -*-
2 
3 // Copyright (C) 2006, 2007, 2008, 2009
4 // Free Software Foundation, Inc.
5 //
6 // This file is part of the GNU ISO C++ Library. This library is free
7 // software; you can redistribute it and/or modify it under the
8 // terms of the GNU General Public License as published by the
9 // Free Software Foundation; either version 3, or (at your option)
10 // any later version.
11 //
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // Under Section 7 of GPL version 3, you are granted additional
18 // permissions described in the GCC Runtime Library Exception, version
19 // 3.1, as published by the Free Software Foundation.
20 
21 // You should have received a copy of the GNU General Public License and
22 // a copy of the GCC Runtime Library Exception along with this program;
23 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 // <http://www.gnu.org/licenses/>.
25 
26 /** @file tr1/gamma.tcc
27  * This is an internal header file, included by other library headers.
28  * You should not attempt to use it directly.
29  */
30 
31 //
32 // ISO C++ 14882 TR1: 5.2 Special functions
33 //
34 
35 // Written by Edward Smith-Rowland based on:
36 // (1) Handbook of Mathematical Functions,
37 // ed. Milton Abramowitz and Irene A. Stegun,
38 // Dover Publications,
39 // Section 6, pp. 253-266
40 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
41 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
42 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
43 // 2nd ed, pp. 213-216
44 // (4) Gamma, Exploring Euler's Constant, Julian Havil,
45 // Princeton, 2003.
46 
47 #ifndef _TR1_GAMMA_TCC
48 #define _TR1_GAMMA_TCC 1
49 
50 #include "special_function_util.h"
51 
52 namespace std
53 {
54 namespace tr1
55 {
56  // Implementation-space details.
57  namespace __detail
58  {
59 
60  /**
61  * @brief This returns Bernoulli numbers from a table or by summation
62  * for larger values.
63  *
64  * Recursion is unstable.
65  *
66  * @param __n the order n of the Bernoulli number.
67  * @return The Bernoulli number of order n.
68  */
69  template <typename _Tp>
70  _Tp __bernoulli_series(unsigned int __n)
71  {
72 
73  static const _Tp __num[28] = {
74  _Tp(1UL), -_Tp(1UL) / _Tp(2UL),
75  _Tp(1UL) / _Tp(6UL), _Tp(0UL),
76  -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
77  _Tp(1UL) / _Tp(42UL), _Tp(0UL),
78  -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
79  _Tp(5UL) / _Tp(66UL), _Tp(0UL),
80  -_Tp(691UL) / _Tp(2730UL), _Tp(0UL),
81  _Tp(7UL) / _Tp(6UL), _Tp(0UL),
82  -_Tp(3617UL) / _Tp(510UL), _Tp(0UL),
83  _Tp(43867UL) / _Tp(798UL), _Tp(0UL),
84  -_Tp(174611) / _Tp(330UL), _Tp(0UL),
85  _Tp(854513UL) / _Tp(138UL), _Tp(0UL),
86  -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),
87  _Tp(8553103UL) / _Tp(6UL), _Tp(0UL)
88  };
89 
90  if (__n == 0)
91  return _Tp(1);
92 
93  if (__n == 1)
94  return -_Tp(1) / _Tp(2);
95 
96  // Take care of the rest of the odd ones.
97  if (__n % 2 == 1)
98  return _Tp(0);
99 
100  // Take care of some small evens that are painful for the series.
101  if (__n < 28)
102  return __num[__n];
103 
104 
105  _Tp __fact = _Tp(1);
106  if ((__n / 2) % 2 == 0)
107  __fact *= _Tp(-1);
108  for (unsigned int __k = 1; __k <= __n; ++__k)
109  __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
110  __fact *= _Tp(2);
111 
112  _Tp __sum = _Tp(0);
113  for (unsigned int __i = 1; __i < 1000; ++__i)
114  {
115  _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
116  if (__term < std::numeric_limits<_Tp>::epsilon())
117  break;
118  __sum += __term;
119  }
120 
121  return __fact * __sum;
122  }
123 
124 
125  /**
126  * @brief This returns Bernoulli number \f$B_n\f$.
127  *
128  * @param __n the order n of the Bernoulli number.
129  * @return The Bernoulli number of order n.
130  */
131  template<typename _Tp>
132  inline _Tp
133  __bernoulli(const int __n)
134  {
135  return __bernoulli_series<_Tp>(__n);
136  }
137 
138 
139  /**
140  * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion
141  * with Bernoulli number coefficients. This is like
142  * Sterling's approximation.
143  *
144  * @param __x The argument of the log of the gamma function.
145  * @return The logarithm of the gamma function.
146  */
147  template<typename _Tp>
148  _Tp
149  __log_gamma_bernoulli(const _Tp __x)
150  {
151  _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
152  + _Tp(0.5L) * std::log(_Tp(2)
153  * __numeric_constants<_Tp>::__pi());
154 
155  const _Tp __xx = __x * __x;
156  _Tp __help = _Tp(1) / __x;
157  for ( unsigned int __i = 1; __i < 20; ++__i )
158  {
159  const _Tp __2i = _Tp(2 * __i);
160  __help /= __2i * (__2i - _Tp(1)) * __xx;
161  __lg += __bernoulli<_Tp>(2 * __i) * __help;
162  }
163 
164  return __lg;
165  }
166 
167 
168  /**
169  * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.
170  * This method dominates all others on the positive axis I think.
171  *
172  * @param __x The argument of the log of the gamma function.
173  * @return The logarithm of the gamma function.
174  */
175  template<typename _Tp>
176  _Tp
177  __log_gamma_lanczos(const _Tp __x)
178  {
179  const _Tp __xm1 = __x - _Tp(1);
180 
181  static const _Tp __lanczos_cheb_7[9] = {
182  _Tp( 0.99999999999980993227684700473478L),
183  _Tp( 676.520368121885098567009190444019L),
184  _Tp(-1259.13921672240287047156078755283L),
185  _Tp( 771.3234287776530788486528258894L),
186  _Tp(-176.61502916214059906584551354L),
187  _Tp( 12.507343278686904814458936853L),
188  _Tp(-0.13857109526572011689554707L),
189  _Tp( 9.984369578019570859563e-6L),
190  _Tp( 1.50563273514931155834e-7L)
191  };
192 
193  static const _Tp __LOGROOT2PI
194  = _Tp(0.9189385332046727417803297364056176L);
195 
196  _Tp __sum = __lanczos_cheb_7[0];
197  for(unsigned int __k = 1; __k < 9; ++__k)
198  __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
199 
200  const _Tp __term1 = (__xm1 + _Tp(0.5L))
201  * std::log((__xm1 + _Tp(7.5L))
202  / __numeric_constants<_Tp>::__euler());
203  const _Tp __term2 = __LOGROOT2PI + std::log(__sum);
204  const _Tp __result = __term1 + (__term2 - _Tp(7));
205 
206  return __result;
207  }
208 
209 
210  /**
211  * @brief Return \f$ log(|\Gamma(x)|) \f$.
212  * This will return values even for \f$ x < 0 \f$.
213  * To recover the sign of \f$ \Gamma(x) \f$ for
214  * any argument use @a __log_gamma_sign.
215  *
216  * @param __x The argument of the log of the gamma function.
217  * @return The logarithm of the gamma function.
218  */
219  template<typename _Tp>
220  _Tp
221  __log_gamma(const _Tp __x)
222  {
223  if (__x > _Tp(0.5L))
224  return __log_gamma_lanczos(__x);
225  else
226  {
227  const _Tp __sin_fact
228  = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));
229  if (__sin_fact == _Tp(0))
230  std::__throw_domain_error(__N("Argument is nonpositive integer "
231  "in __log_gamma"));
232  return __numeric_constants<_Tp>::__lnpi()
233  - std::log(__sin_fact)
234  - __log_gamma_lanczos(_Tp(1) - __x);
235  }
236  }
237 
238 
239  /**
240  * @brief Return the sign of \f$ \Gamma(x) \f$.
241  * At nonpositive integers zero is returned.
242  *
243  * @param __x The argument of the gamma function.
244  * @return The sign of the gamma function.
245  */
246  template<typename _Tp>
247  _Tp
248  __log_gamma_sign(const _Tp __x)
249  {
250  if (__x > _Tp(0))
251  return _Tp(1);
252  else
253  {
254  const _Tp __sin_fact
255  = std::sin(__numeric_constants<_Tp>::__pi() * __x);
256  if (__sin_fact > _Tp(0))
257  return (1);
258  else if (__sin_fact < _Tp(0))
259  return -_Tp(1);
260  else
261  return _Tp(0);
262  }
263  }
264 
265 
266  /**
267  * @brief Return the logarithm of the binomial coefficient.
268  * The binomial coefficient is given by:
269  * @f[
270  * \left( \right) = \frac{n!}{(n-k)! k!}
271  * @f]
272  *
273  * @param __n The first argument of the binomial coefficient.
274  * @param __k The second argument of the binomial coefficient.
275  * @return The binomial coefficient.
276  */
277  template<typename _Tp>
278  _Tp
279  __log_bincoef(const unsigned int __n, const unsigned int __k)
280  {
281  // Max e exponent before overflow.
282  static const _Tp __max_bincoeff
283  = std::numeric_limits<_Tp>::max_exponent10
284  * std::log(_Tp(10)) - _Tp(1);
285 #if _GLIBCXX_USE_C99_MATH_TR1
286  _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n))
287  - std::tr1::lgamma(_Tp(1 + __k))
288  - std::tr1::lgamma(_Tp(1 + __n - __k));
289 #else
290  _Tp __coeff = __log_gamma(_Tp(1 + __n))
291  - __log_gamma(_Tp(1 + __k))
292  - __log_gamma(_Tp(1 + __n - __k));
293 #endif
294  }
295 
296 
297  /**
298  * @brief Return the binomial coefficient.
299  * The binomial coefficient is given by:
300  * @f[
301  * \left( \right) = \frac{n!}{(n-k)! k!}
302  * @f]
303  *
304  * @param __n The first argument of the binomial coefficient.
305  * @param __k The second argument of the binomial coefficient.
306  * @return The binomial coefficient.
307  */
308  template<typename _Tp>
309  _Tp
310  __bincoef(const unsigned int __n, const unsigned int __k)
311  {
312  // Max e exponent before overflow.
313  static const _Tp __max_bincoeff
314  = std::numeric_limits<_Tp>::max_exponent10
315  * std::log(_Tp(10)) - _Tp(1);
316 
317  const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
318  if (__log_coeff > __max_bincoeff)
319  return std::numeric_limits<_Tp>::quiet_NaN();
320  else
321  return std::exp(__log_coeff);
322  }
323 
324 
325  /**
326  * @brief Return \f$ \Gamma(x) \f$.
327  *
328  * @param __x The argument of the gamma function.
329  * @return The gamma function.
330  */
331  template<typename _Tp>
332  inline _Tp
333  __gamma(const _Tp __x)
334  {
335  return std::exp(__log_gamma(__x));
336  }
337 
338 
339  /**
340  * @brief Return the digamma function by series expansion.
341  * The digamma or @f$ \psi(x) @f$ function is defined by
342  * @f[
343  * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
344  * @f]
345  *
346  * The series is given by:
347  * @f[
348  * \psi(x) = -\gamma_E - \frac{1}{x}
349  * \sum_{k=1}^{\infty} \frac{x}{k(x + k)}
350  * @f]
351  */
352  template<typename _Tp>
353  _Tp
354  __psi_series(const _Tp __x)
355  {
356  _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
357  const unsigned int __max_iter = 100000;
358  for (unsigned int __k = 1; __k < __max_iter; ++__k)
359  {
360  const _Tp __term = __x / (__k * (__k + __x));
361  __sum += __term;
362  if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
363  break;
364  }
365  return __sum;
366  }
367 
368 
369  /**
370  * @brief Return the digamma function for large argument.
371  * The digamma or @f$ \psi(x) @f$ function is defined by
372  * @f[
373  * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
374  * @f]
375  *
376  * The asymptotic series is given by:
377  * @f[
378  * \psi(x) = \ln(x) - \frac{1}{2x}
379  * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}
380  * @f]
381  */
382  template<typename _Tp>
383  _Tp
384  __psi_asymp(const _Tp __x)
385  {
386  _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
387  const _Tp __xx = __x * __x;
388  _Tp __xp = __xx;
389  const unsigned int __max_iter = 100;
390  for (unsigned int __k = 1; __k < __max_iter; ++__k)
391  {
392  const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
393  __sum -= __term;
394  if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
395  break;
396  __xp *= __xx;
397  }
398  return __sum;
399  }
400 
401 
402  /**
403  * @brief Return the digamma function.
404  * The digamma or @f$ \psi(x) @f$ function is defined by
405  * @f[
406  * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
407  * @f]
408  * For negative argument the reflection formula is used:
409  * @f[
410  * \psi(x) = \psi(1-x) - \pi \cot(\pi x)
411  * @f]
412  */
413  template<typename _Tp>
414  _Tp
415  __psi(const _Tp __x)
416  {
417  const int __n = static_cast<int>(__x + 0.5L);
418  const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
419  if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)
420  return std::numeric_limits<_Tp>::quiet_NaN();
421  else if (__x < _Tp(0))
422  {
423  const _Tp __pi = __numeric_constants<_Tp>::__pi();
424  return __psi(_Tp(1) - __x)
425  - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
426  }
427  else if (__x > _Tp(100))
428  return __psi_asymp(__x);
429  else
430  return __psi_series(__x);
431  }
432 
433 
434  /**
435  * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$.
436  *
437  * The polygamma function is related to the Hurwitz zeta function:
438  * @f[
439  * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)
440  * @f]
441  */
442  template<typename _Tp>
443  _Tp
444  __psi(const unsigned int __n, const _Tp __x)
445  {
446  if (__x <= _Tp(0))
447  std::__throw_domain_error(__N("Argument out of range "
448  "in __psi"));
449  else if (__n == 0)
450  return __psi(__x);
451  else
452  {
453  const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);
454 #if _GLIBCXX_USE_C99_MATH_TR1
455  const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
456 #else
457  const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
458 #endif
459  _Tp __result = std::exp(__ln_nfact) * __hzeta;
460  if (__n % 2 == 1)
461  __result = -__result;
462  return __result;
463  }
464  }
465 
466  } // namespace std::tr1::__detail
467 }
468 }
469 
470 #endif // _TR1_GAMMA_TCC
471