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