| 1 | // Special functions -*- C++ -*-
|
|---|
| 2 |
|
|---|
| 3 | // Copyright (C) 2006-2021 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 |
|
|---|
| 51 | namespace 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)
|
|---|
| 58 | namespace 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 |
|
|---|
