47 #ifndef _TR1_GAMMA_TCC
48 #define _TR1_GAMMA_TCC 1
69 template <
typename _Tp>
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)
94 return -_Tp(1) / _Tp(2);
106 if ((__n / 2) % 2 == 0)
108 for (
unsigned int __k = 1; __k <= __n; ++__k)
113 for (
unsigned int __i = 1; __i < 1000; ++__i)
115 _Tp __term =
std::pow(_Tp(__i), -_Tp(__n));
121 return __fact * __sum;
131 template<
typename _Tp>
135 return __bernoulli_series<_Tp>(__n);
147 template<
typename _Tp>
155 const _Tp __xx = __x * __x;
156 _Tp __help = _Tp(1) / __x;
157 for (
unsigned int __i = 1; __i < 20; ++__i )
159 const _Tp __2i = _Tp(2 * __i);
160 __help /= __2i * (__2i - _Tp(1)) * __xx;
161 __lg += __bernoulli<_Tp>(2 * __i) * __help;
175 template<
typename _Tp>
179 const _Tp __xm1 = __x - _Tp(1);
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)
193 static const _Tp __LOGROOT2PI
194 = _Tp(0.9189385332046727417803297364056176L);
196 _Tp __sum = __lanczos_cheb_7[0];
197 for(
unsigned int __k = 1; __k < 9; ++__k)
198 __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
200 const _Tp __term1 = (__xm1 + _Tp(0.5L))
203 const _Tp __term2 = __LOGROOT2PI +
std::log(__sum);
204 const _Tp __result = __term1 + (__term2 - _Tp(7));
219 template<
typename _Tp>
229 if (__sin_fact == _Tp(0))
230 std::__throw_domain_error(__N(
"Argument is nonpositive integer "
246 template<
typename _Tp>
256 if (__sin_fact > _Tp(0))
258 else if (__sin_fact < _Tp(0))
277 template<
typename _Tp>
282 static const _Tp __max_bincoeff
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));
308 template<
typename _Tp>
310 __bincoef(
const unsigned int __n,
const unsigned int __k)
313 static const _Tp __max_bincoeff
317 const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
318 if (__log_coeff > __max_bincoeff)
331 template<
typename _Tp>
352 template<
typename _Tp>
357 const unsigned int __max_iter = 100000;
358 for (
unsigned int __k = 1; __k < __max_iter; ++__k)
360 const _Tp __term = __x / (__k * (__k + __x));
382 template<
typename _Tp>
386 _Tp __sum =
std::log(__x) - _Tp(0.5L) / __x;
387 const _Tp __xx = __x * __x;
389 const unsigned int __max_iter = 100;
390 for (
unsigned int __k = 1; __k < __max_iter; ++__k)
392 const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
413 template<
typename _Tp>
417 const int __n =
static_cast<int>(__x + 0.5L);
419 if (__n <= 0 &&
std::abs(__x - _Tp(__n)) < __eps)
421 else if (__x < _Tp(0))
424 return __psi(_Tp(1) - __x)
427 else if (__x > _Tp(100))
442 template<
typename _Tp>
444 __psi(
const unsigned int __n,
const _Tp __x)
447 std::__throw_domain_error(__N(
"Argument out of range "
454 #if _GLIBCXX_USE_C99_MATH_TR1
455 const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
459 _Tp __result =
std::exp(__ln_nfact) * __hzeta;
461 __result = -__result;
470 #endif // _TR1_GAMMA_TCC
_Tp __psi_asymp(const _Tp __x)
Return the digamma function for large argument. The digamma or function is defined by ...
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
_Tp __log_gamma_bernoulli(const _Tp __x)
Return by asymptotic expansion with Bernoulli number coefficients. This is like Sterling's approxima...
A structure for numeric constants.
_Tp __bernoulli_series(unsigned int __n)
This returns Bernoulli numbers from a table or by summation for larger values.
static _Tp __gamma_e()
Constant Euler's constant .
complex< _Tp > cos(const complex< _Tp > &)
Return complex cosine of z.
_Tp __psi_series(const _Tp __x)
Return the digamma function by series expansion. The digamma or function is defined by ...
_Tp __hurwitz_zeta(const _Tp __a, const _Tp __s)
Return the Hurwitz zeta function for all s != 1 and x > -1.
_Tp __bincoef(const unsigned int __n, const unsigned int __k)
Return the binomial coefficient. The binomial coefficient is given by: .
Properties of fundamental types.
_Tp __log_bincoef(const unsigned int __n, const unsigned int __k)
Return the logarithm of the binomial coefficient. The binomial coefficient is given by: ...
_Tp __psi(const _Tp __x)
Return the digamma function. The digamma or function is defined by For negative argument the reflec...
_Tp __log_gamma_lanczos(const _Tp __x)
Return by the Lanczos method. This method dominates all others on the positive axis I think...
static _Tp __lnpi()
Constant .
_Tp __log_gamma(const _Tp __x)
Return . This will return values even for . To recover the sign of for any argument use __log_gamma_...
complex< _Tp > sin(const complex< _Tp > &)
Return complex sine of z.
complex< _Tp > pow(const complex< _Tp > &, const _Tp &)
Return x to the y'th power.
_Size __lg(_Size __n)
This is a helper function for the sort routines. Precondition: __n > 0.
_Tp __bernoulli(const int __n)
This returns Bernoulli number .
static _Tp __pi()
Constant .
static _Tp __euler()
Constant Euler-Mascheroni .
_Tp __log_gamma_sign(const _Tp __x)
Return the sign of . At nonpositive integers zero is returned.
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
_Tp __gamma(const _Tp __x)
Return .
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.