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