1*e4b17023SJohn Marino // Special functions -*- C++ -*- 2*e4b17023SJohn Marino 3*e4b17023SJohn Marino // Copyright (C) 2006, 2007, 2008, 2009, 2010 4*e4b17023SJohn Marino // Free Software Foundation, Inc. 5*e4b17023SJohn Marino // 6*e4b17023SJohn Marino // This file is part of the GNU ISO C++ Library. This library is free 7*e4b17023SJohn Marino // software; you can redistribute it and/or modify it under the 8*e4b17023SJohn Marino // terms of the GNU General Public License as published by the 9*e4b17023SJohn Marino // Free Software Foundation; either version 3, or (at your option) 10*e4b17023SJohn Marino // any later version. 11*e4b17023SJohn Marino // 12*e4b17023SJohn Marino // This library is distributed in the hope that it will be useful, 13*e4b17023SJohn Marino // but WITHOUT ANY WARRANTY; without even the implied warranty of 14*e4b17023SJohn Marino // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15*e4b17023SJohn Marino // GNU General Public License for more details. 16*e4b17023SJohn Marino // 17*e4b17023SJohn Marino // Under Section 7 of GPL version 3, you are granted additional 18*e4b17023SJohn Marino // permissions described in the GCC Runtime Library Exception, version 19*e4b17023SJohn Marino // 3.1, as published by the Free Software Foundation. 20*e4b17023SJohn Marino 21*e4b17023SJohn Marino // You should have received a copy of the GNU General Public License and 22*e4b17023SJohn Marino // a copy of the GCC Runtime Library Exception along with this program; 23*e4b17023SJohn Marino // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24*e4b17023SJohn Marino // <http://www.gnu.org/licenses/>. 25*e4b17023SJohn Marino 26*e4b17023SJohn Marino /** @file tr1/riemann_zeta.tcc 27*e4b17023SJohn Marino * This is an internal header file, included by other library headers. 28*e4b17023SJohn Marino * Do not attempt to use it directly. @headername{tr1/cmath} 29*e4b17023SJohn Marino */ 30*e4b17023SJohn Marino 31*e4b17023SJohn Marino // 32*e4b17023SJohn Marino // ISO C++ 14882 TR1: 5.2 Special functions 33*e4b17023SJohn Marino // 34*e4b17023SJohn Marino 35*e4b17023SJohn Marino // Written by Edward Smith-Rowland based on: 36*e4b17023SJohn Marino // (1) Handbook of Mathematical Functions, 37*e4b17023SJohn Marino // Ed. by Milton Abramowitz and Irene A. Stegun, 38*e4b17023SJohn Marino // Dover Publications, New-York, Section 5, pp. 807-808. 39*e4b17023SJohn Marino // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 40*e4b17023SJohn Marino // (3) Gamma, Exploring Euler's Constant, Julian Havil, 41*e4b17023SJohn Marino // Princeton, 2003. 42*e4b17023SJohn Marino 43*e4b17023SJohn Marino #ifndef _GLIBCXX_TR1_RIEMANN_ZETA_TCC 44*e4b17023SJohn Marino #define _GLIBCXX_TR1_RIEMANN_ZETA_TCC 1 45*e4b17023SJohn Marino 46*e4b17023SJohn Marino #include "special_function_util.h" 47*e4b17023SJohn Marino 48*e4b17023SJohn Marino namespace std _GLIBCXX_VISIBILITY(default) 49*e4b17023SJohn Marino { 50*e4b17023SJohn Marino namespace tr1 51*e4b17023SJohn Marino { 52*e4b17023SJohn Marino // [5.2] Special functions 53*e4b17023SJohn Marino 54*e4b17023SJohn Marino // Implementation-space details. 55*e4b17023SJohn Marino namespace __detail 56*e4b17023SJohn Marino { 57*e4b17023SJohn Marino _GLIBCXX_BEGIN_NAMESPACE_VERSION 58*e4b17023SJohn Marino 59*e4b17023SJohn Marino /** 60*e4b17023SJohn Marino * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$ 61*e4b17023SJohn Marino * by summation for s > 1. 62*e4b17023SJohn Marino * 63*e4b17023SJohn Marino * The Riemann zeta function is defined by: 64*e4b17023SJohn Marino * \f[ 65*e4b17023SJohn Marino * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 66*e4b17023SJohn Marino * \f] 67*e4b17023SJohn Marino * For s < 1 use the reflection formula: 68*e4b17023SJohn Marino * \f[ 69*e4b17023SJohn Marino * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) 70*e4b17023SJohn Marino * \f] 71*e4b17023SJohn Marino */ 72*e4b17023SJohn Marino template<typename _Tp> 73*e4b17023SJohn Marino _Tp __riemann_zeta_sum(const _Tp __s)74*e4b17023SJohn Marino __riemann_zeta_sum(const _Tp __s) 75*e4b17023SJohn Marino { 76*e4b17023SJohn Marino // A user shouldn't get to this. 77*e4b17023SJohn Marino if (__s < _Tp(1)) 78*e4b17023SJohn Marino std::__throw_domain_error(__N("Bad argument in zeta sum.")); 79*e4b17023SJohn Marino 80*e4b17023SJohn Marino const unsigned int max_iter = 10000; 81*e4b17023SJohn Marino _Tp __zeta = _Tp(0); 82*e4b17023SJohn Marino for (unsigned int __k = 1; __k < max_iter; ++__k) 83*e4b17023SJohn Marino { 84*e4b17023SJohn Marino _Tp __term = std::pow(static_cast<_Tp>(__k), -__s); 85*e4b17023SJohn Marino if (__term < std::numeric_limits<_Tp>::epsilon()) 86*e4b17023SJohn Marino { 87*e4b17023SJohn Marino break; 88*e4b17023SJohn Marino } 89*e4b17023SJohn Marino __zeta += __term; 90*e4b17023SJohn Marino } 91*e4b17023SJohn Marino 92*e4b17023SJohn Marino return __zeta; 93*e4b17023SJohn Marino } 94*e4b17023SJohn Marino 95*e4b17023SJohn Marino 96*e4b17023SJohn Marino /** 97*e4b17023SJohn Marino * @brief Evaluate the Riemann zeta function @f$ \zeta(s) @f$ 98*e4b17023SJohn Marino * by an alternate series for s > 0. 99*e4b17023SJohn Marino * 100*e4b17023SJohn Marino * The Riemann zeta function is defined by: 101*e4b17023SJohn Marino * \f[ 102*e4b17023SJohn Marino * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 103*e4b17023SJohn Marino * \f] 104*e4b17023SJohn Marino * For s < 1 use the reflection formula: 105*e4b17023SJohn Marino * \f[ 106*e4b17023SJohn Marino * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) 107*e4b17023SJohn Marino * \f] 108*e4b17023SJohn Marino */ 109*e4b17023SJohn Marino template<typename _Tp> 110*e4b17023SJohn Marino _Tp __riemann_zeta_alt(const _Tp __s)111*e4b17023SJohn Marino __riemann_zeta_alt(const _Tp __s) 112*e4b17023SJohn Marino { 113*e4b17023SJohn Marino _Tp __sgn = _Tp(1); 114*e4b17023SJohn Marino _Tp __zeta = _Tp(0); 115*e4b17023SJohn Marino for (unsigned int __i = 1; __i < 10000000; ++__i) 116*e4b17023SJohn Marino { 117*e4b17023SJohn Marino _Tp __term = __sgn / std::pow(__i, __s); 118*e4b17023SJohn Marino if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) 119*e4b17023SJohn Marino break; 120*e4b17023SJohn Marino __zeta += __term; 121*e4b17023SJohn Marino __sgn *= _Tp(-1); 122*e4b17023SJohn Marino } 123*e4b17023SJohn Marino __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s); 124*e4b17023SJohn Marino 125*e4b17023SJohn Marino return __zeta; 126*e4b17023SJohn Marino } 127*e4b17023SJohn Marino 128*e4b17023SJohn Marino 129*e4b17023SJohn Marino /** 130*e4b17023SJohn Marino * @brief Evaluate the Riemann zeta function by series for all s != 1. 131*e4b17023SJohn Marino * Convergence is great until largish negative numbers. 132*e4b17023SJohn Marino * Then the convergence of the > 0 sum gets better. 133*e4b17023SJohn Marino * 134*e4b17023SJohn Marino * The series is: 135*e4b17023SJohn Marino * \f[ 136*e4b17023SJohn Marino * \zeta(s) = \frac{1}{1-2^{1-s}} 137*e4b17023SJohn Marino * \sum_{n=0}^{\infty} \frac{1}{2^{n+1}} 138*e4b17023SJohn Marino * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (k+1)^{-s} 139*e4b17023SJohn Marino * \f] 140*e4b17023SJohn Marino * Havil 2003, p. 206. 141*e4b17023SJohn Marino * 142*e4b17023SJohn Marino * The Riemann zeta function is defined by: 143*e4b17023SJohn Marino * \f[ 144*e4b17023SJohn Marino * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 145*e4b17023SJohn Marino * \f] 146*e4b17023SJohn Marino * For s < 1 use the reflection formula: 147*e4b17023SJohn Marino * \f[ 148*e4b17023SJohn Marino * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) 149*e4b17023SJohn Marino * \f] 150*e4b17023SJohn Marino */ 151*e4b17023SJohn Marino template<typename _Tp> 152*e4b17023SJohn Marino _Tp __riemann_zeta_glob(const _Tp __s)153*e4b17023SJohn Marino __riemann_zeta_glob(const _Tp __s) 154*e4b17023SJohn Marino { 155*e4b17023SJohn Marino _Tp __zeta = _Tp(0); 156*e4b17023SJohn Marino 157*e4b17023SJohn Marino const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 158*e4b17023SJohn Marino // Max e exponent before overflow. 159*e4b17023SJohn Marino const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10 160*e4b17023SJohn Marino * std::log(_Tp(10)) - _Tp(1); 161*e4b17023SJohn Marino 162*e4b17023SJohn Marino // This series works until the binomial coefficient blows up 163*e4b17023SJohn Marino // so use reflection. 164*e4b17023SJohn Marino if (__s < _Tp(0)) 165*e4b17023SJohn Marino { 166*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 167*e4b17023SJohn Marino if (std::tr1::fmod(__s,_Tp(2)) == _Tp(0)) 168*e4b17023SJohn Marino return _Tp(0); 169*e4b17023SJohn Marino else 170*e4b17023SJohn Marino #endif 171*e4b17023SJohn Marino { 172*e4b17023SJohn Marino _Tp __zeta = __riemann_zeta_glob(_Tp(1) - __s); 173*e4b17023SJohn Marino __zeta *= std::pow(_Tp(2) 174*e4b17023SJohn Marino * __numeric_constants<_Tp>::__pi(), __s) 175*e4b17023SJohn Marino * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) 176*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 177*e4b17023SJohn Marino * std::exp(std::tr1::lgamma(_Tp(1) - __s)) 178*e4b17023SJohn Marino #else 179*e4b17023SJohn Marino * std::exp(__log_gamma(_Tp(1) - __s)) 180*e4b17023SJohn Marino #endif 181*e4b17023SJohn Marino / __numeric_constants<_Tp>::__pi(); 182*e4b17023SJohn Marino return __zeta; 183*e4b17023SJohn Marino } 184*e4b17023SJohn Marino } 185*e4b17023SJohn Marino 186*e4b17023SJohn Marino _Tp __num = _Tp(0.5L); 187*e4b17023SJohn Marino const unsigned int __maxit = 10000; 188*e4b17023SJohn Marino for (unsigned int __i = 0; __i < __maxit; ++__i) 189*e4b17023SJohn Marino { 190*e4b17023SJohn Marino bool __punt = false; 191*e4b17023SJohn Marino _Tp __sgn = _Tp(1); 192*e4b17023SJohn Marino _Tp __term = _Tp(0); 193*e4b17023SJohn Marino for (unsigned int __j = 0; __j <= __i; ++__j) 194*e4b17023SJohn Marino { 195*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 196*e4b17023SJohn Marino _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i)) 197*e4b17023SJohn Marino - std::tr1::lgamma(_Tp(1 + __j)) 198*e4b17023SJohn Marino - std::tr1::lgamma(_Tp(1 + __i - __j)); 199*e4b17023SJohn Marino #else 200*e4b17023SJohn Marino _Tp __bincoeff = __log_gamma(_Tp(1 + __i)) 201*e4b17023SJohn Marino - __log_gamma(_Tp(1 + __j)) 202*e4b17023SJohn Marino - __log_gamma(_Tp(1 + __i - __j)); 203*e4b17023SJohn Marino #endif 204*e4b17023SJohn Marino if (__bincoeff > __max_bincoeff) 205*e4b17023SJohn Marino { 206*e4b17023SJohn Marino // This only gets hit for x << 0. 207*e4b17023SJohn Marino __punt = true; 208*e4b17023SJohn Marino break; 209*e4b17023SJohn Marino } 210*e4b17023SJohn Marino __bincoeff = std::exp(__bincoeff); 211*e4b17023SJohn Marino __term += __sgn * __bincoeff * std::pow(_Tp(1 + __j), -__s); 212*e4b17023SJohn Marino __sgn *= _Tp(-1); 213*e4b17023SJohn Marino } 214*e4b17023SJohn Marino if (__punt) 215*e4b17023SJohn Marino break; 216*e4b17023SJohn Marino __term *= __num; 217*e4b17023SJohn Marino __zeta += __term; 218*e4b17023SJohn Marino if (std::abs(__term/__zeta) < __eps) 219*e4b17023SJohn Marino break; 220*e4b17023SJohn Marino __num *= _Tp(0.5L); 221*e4b17023SJohn Marino } 222*e4b17023SJohn Marino 223*e4b17023SJohn Marino __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s); 224*e4b17023SJohn Marino 225*e4b17023SJohn Marino return __zeta; 226*e4b17023SJohn Marino } 227*e4b17023SJohn Marino 228*e4b17023SJohn Marino 229*e4b17023SJohn Marino /** 230*e4b17023SJohn Marino * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$ 231*e4b17023SJohn Marino * using the product over prime factors. 232*e4b17023SJohn Marino * \f[ 233*e4b17023SJohn Marino * \zeta(s) = \Pi_{i=1}^\infty \frac{1}{1 - p_i^{-s}} 234*e4b17023SJohn Marino * \f] 235*e4b17023SJohn Marino * where @f$ {p_i} @f$ are the prime numbers. 236*e4b17023SJohn Marino * 237*e4b17023SJohn Marino * The Riemann zeta function is defined by: 238*e4b17023SJohn Marino * \f[ 239*e4b17023SJohn Marino * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 240*e4b17023SJohn Marino * \f] 241*e4b17023SJohn Marino * For s < 1 use the reflection formula: 242*e4b17023SJohn Marino * \f[ 243*e4b17023SJohn Marino * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) 244*e4b17023SJohn Marino * \f] 245*e4b17023SJohn Marino */ 246*e4b17023SJohn Marino template<typename _Tp> 247*e4b17023SJohn Marino _Tp __riemann_zeta_product(const _Tp __s)248*e4b17023SJohn Marino __riemann_zeta_product(const _Tp __s) 249*e4b17023SJohn Marino { 250*e4b17023SJohn Marino static const _Tp __prime[] = { 251*e4b17023SJohn Marino _Tp(2), _Tp(3), _Tp(5), _Tp(7), _Tp(11), _Tp(13), _Tp(17), _Tp(19), 252*e4b17023SJohn Marino _Tp(23), _Tp(29), _Tp(31), _Tp(37), _Tp(41), _Tp(43), _Tp(47), 253*e4b17023SJohn Marino _Tp(53), _Tp(59), _Tp(61), _Tp(67), _Tp(71), _Tp(73), _Tp(79), 254*e4b17023SJohn Marino _Tp(83), _Tp(89), _Tp(97), _Tp(101), _Tp(103), _Tp(107), _Tp(109) 255*e4b17023SJohn Marino }; 256*e4b17023SJohn Marino static const unsigned int __num_primes = sizeof(__prime) / sizeof(_Tp); 257*e4b17023SJohn Marino 258*e4b17023SJohn Marino _Tp __zeta = _Tp(1); 259*e4b17023SJohn Marino for (unsigned int __i = 0; __i < __num_primes; ++__i) 260*e4b17023SJohn Marino { 261*e4b17023SJohn Marino const _Tp __fact = _Tp(1) - std::pow(__prime[__i], -__s); 262*e4b17023SJohn Marino __zeta *= __fact; 263*e4b17023SJohn Marino if (_Tp(1) - __fact < std::numeric_limits<_Tp>::epsilon()) 264*e4b17023SJohn Marino break; 265*e4b17023SJohn Marino } 266*e4b17023SJohn Marino 267*e4b17023SJohn Marino __zeta = _Tp(1) / __zeta; 268*e4b17023SJohn Marino 269*e4b17023SJohn Marino return __zeta; 270*e4b17023SJohn Marino } 271*e4b17023SJohn Marino 272*e4b17023SJohn Marino 273*e4b17023SJohn Marino /** 274*e4b17023SJohn Marino * @brief Return the Riemann zeta function @f$ \zeta(s) @f$. 275*e4b17023SJohn Marino * 276*e4b17023SJohn Marino * The Riemann zeta function is defined by: 277*e4b17023SJohn Marino * \f[ 278*e4b17023SJohn Marino * \zeta(s) = \sum_{k=1}^{\infty} k^{-s} for s > 1 279*e4b17023SJohn Marino * \frac{(2\pi)^s}{pi} sin(\frac{\pi s}{2}) 280*e4b17023SJohn Marino * \Gamma (1 - s) \zeta (1 - s) for s < 1 281*e4b17023SJohn Marino * \f] 282*e4b17023SJohn Marino * For s < 1 use the reflection formula: 283*e4b17023SJohn Marino * \f[ 284*e4b17023SJohn Marino * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) 285*e4b17023SJohn Marino * \f] 286*e4b17023SJohn Marino */ 287*e4b17023SJohn Marino template<typename _Tp> 288*e4b17023SJohn Marino _Tp __riemann_zeta(const _Tp __s)289*e4b17023SJohn Marino __riemann_zeta(const _Tp __s) 290*e4b17023SJohn Marino { 291*e4b17023SJohn Marino if (__isnan(__s)) 292*e4b17023SJohn Marino return std::numeric_limits<_Tp>::quiet_NaN(); 293*e4b17023SJohn Marino else if (__s == _Tp(1)) 294*e4b17023SJohn Marino return std::numeric_limits<_Tp>::infinity(); 295*e4b17023SJohn Marino else if (__s < -_Tp(19)) 296*e4b17023SJohn Marino { 297*e4b17023SJohn Marino _Tp __zeta = __riemann_zeta_product(_Tp(1) - __s); 298*e4b17023SJohn Marino __zeta *= std::pow(_Tp(2) * __numeric_constants<_Tp>::__pi(), __s) 299*e4b17023SJohn Marino * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) 300*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 301*e4b17023SJohn Marino * std::exp(std::tr1::lgamma(_Tp(1) - __s)) 302*e4b17023SJohn Marino #else 303*e4b17023SJohn Marino * std::exp(__log_gamma(_Tp(1) - __s)) 304*e4b17023SJohn Marino #endif 305*e4b17023SJohn Marino / __numeric_constants<_Tp>::__pi(); 306*e4b17023SJohn Marino return __zeta; 307*e4b17023SJohn Marino } 308*e4b17023SJohn Marino else if (__s < _Tp(20)) 309*e4b17023SJohn Marino { 310*e4b17023SJohn Marino // Global double sum or McLaurin? 311*e4b17023SJohn Marino bool __glob = true; 312*e4b17023SJohn Marino if (__glob) 313*e4b17023SJohn Marino return __riemann_zeta_glob(__s); 314*e4b17023SJohn Marino else 315*e4b17023SJohn Marino { 316*e4b17023SJohn Marino if (__s > _Tp(1)) 317*e4b17023SJohn Marino return __riemann_zeta_sum(__s); 318*e4b17023SJohn Marino else 319*e4b17023SJohn Marino { 320*e4b17023SJohn Marino _Tp __zeta = std::pow(_Tp(2) 321*e4b17023SJohn Marino * __numeric_constants<_Tp>::__pi(), __s) 322*e4b17023SJohn Marino * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) 323*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 324*e4b17023SJohn Marino * std::tr1::tgamma(_Tp(1) - __s) 325*e4b17023SJohn Marino #else 326*e4b17023SJohn Marino * std::exp(__log_gamma(_Tp(1) - __s)) 327*e4b17023SJohn Marino #endif 328*e4b17023SJohn Marino * __riemann_zeta_sum(_Tp(1) - __s); 329*e4b17023SJohn Marino return __zeta; 330*e4b17023SJohn Marino } 331*e4b17023SJohn Marino } 332*e4b17023SJohn Marino } 333*e4b17023SJohn Marino else 334*e4b17023SJohn Marino return __riemann_zeta_product(__s); 335*e4b17023SJohn Marino } 336*e4b17023SJohn Marino 337*e4b17023SJohn Marino 338*e4b17023SJohn Marino /** 339*e4b17023SJohn Marino * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$ 340*e4b17023SJohn Marino * for all s != 1 and x > -1. 341*e4b17023SJohn Marino * 342*e4b17023SJohn Marino * The Hurwitz zeta function is defined by: 343*e4b17023SJohn Marino * @f[ 344*e4b17023SJohn Marino * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s} 345*e4b17023SJohn Marino * @f] 346*e4b17023SJohn Marino * The Riemann zeta function is a special case: 347*e4b17023SJohn Marino * @f[ 348*e4b17023SJohn Marino * \zeta(s) = \zeta(1,s) 349*e4b17023SJohn Marino * @f] 350*e4b17023SJohn Marino * 351*e4b17023SJohn Marino * This functions uses the double sum that converges for s != 1 352*e4b17023SJohn Marino * and x > -1: 353*e4b17023SJohn Marino * @f[ 354*e4b17023SJohn Marino * \zeta(x,s) = \frac{1}{s-1} 355*e4b17023SJohn Marino * \sum_{n=0}^{\infty} \frac{1}{n + 1} 356*e4b17023SJohn Marino * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (x+k)^{-s} 357*e4b17023SJohn Marino * @f] 358*e4b17023SJohn Marino */ 359*e4b17023SJohn Marino template<typename _Tp> 360*e4b17023SJohn Marino _Tp __hurwitz_zeta_glob(const _Tp __a,const _Tp __s)361*e4b17023SJohn Marino __hurwitz_zeta_glob(const _Tp __a, const _Tp __s) 362*e4b17023SJohn Marino { 363*e4b17023SJohn Marino _Tp __zeta = _Tp(0); 364*e4b17023SJohn Marino 365*e4b17023SJohn Marino const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 366*e4b17023SJohn Marino // Max e exponent before overflow. 367*e4b17023SJohn Marino const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10 368*e4b17023SJohn Marino * std::log(_Tp(10)) - _Tp(1); 369*e4b17023SJohn Marino 370*e4b17023SJohn Marino const unsigned int __maxit = 10000; 371*e4b17023SJohn Marino for (unsigned int __i = 0; __i < __maxit; ++__i) 372*e4b17023SJohn Marino { 373*e4b17023SJohn Marino bool __punt = false; 374*e4b17023SJohn Marino _Tp __sgn = _Tp(1); 375*e4b17023SJohn Marino _Tp __term = _Tp(0); 376*e4b17023SJohn Marino for (unsigned int __j = 0; __j <= __i; ++__j) 377*e4b17023SJohn Marino { 378*e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 379*e4b17023SJohn Marino _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i)) 380*e4b17023SJohn Marino - std::tr1::lgamma(_Tp(1 + __j)) 381*e4b17023SJohn Marino - std::tr1::lgamma(_Tp(1 + __i - __j)); 382*e4b17023SJohn Marino #else 383*e4b17023SJohn Marino _Tp __bincoeff = __log_gamma(_Tp(1 + __i)) 384*e4b17023SJohn Marino - __log_gamma(_Tp(1 + __j)) 385*e4b17023SJohn Marino - __log_gamma(_Tp(1 + __i - __j)); 386*e4b17023SJohn Marino #endif 387*e4b17023SJohn Marino if (__bincoeff > __max_bincoeff) 388*e4b17023SJohn Marino { 389*e4b17023SJohn Marino // This only gets hit for x << 0. 390*e4b17023SJohn Marino __punt = true; 391*e4b17023SJohn Marino break; 392*e4b17023SJohn Marino } 393*e4b17023SJohn Marino __bincoeff = std::exp(__bincoeff); 394*e4b17023SJohn Marino __term += __sgn * __bincoeff * std::pow(_Tp(__a + __j), -__s); 395*e4b17023SJohn Marino __sgn *= _Tp(-1); 396*e4b17023SJohn Marino } 397*e4b17023SJohn Marino if (__punt) 398*e4b17023SJohn Marino break; 399*e4b17023SJohn Marino __term /= _Tp(__i + 1); 400*e4b17023SJohn Marino if (std::abs(__term / __zeta) < __eps) 401*e4b17023SJohn Marino break; 402*e4b17023SJohn Marino __zeta += __term; 403*e4b17023SJohn Marino } 404*e4b17023SJohn Marino 405*e4b17023SJohn Marino __zeta /= __s - _Tp(1); 406*e4b17023SJohn Marino 407*e4b17023SJohn Marino return __zeta; 408*e4b17023SJohn Marino } 409*e4b17023SJohn Marino 410*e4b17023SJohn Marino 411*e4b17023SJohn Marino /** 412*e4b17023SJohn Marino * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$ 413*e4b17023SJohn Marino * for all s != 1 and x > -1. 414*e4b17023SJohn Marino * 415*e4b17023SJohn Marino * The Hurwitz zeta function is defined by: 416*e4b17023SJohn Marino * @f[ 417*e4b17023SJohn Marino * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s} 418*e4b17023SJohn Marino * @f] 419*e4b17023SJohn Marino * The Riemann zeta function is a special case: 420*e4b17023SJohn Marino * @f[ 421*e4b17023SJohn Marino * \zeta(s) = \zeta(1,s) 422*e4b17023SJohn Marino * @f] 423*e4b17023SJohn Marino */ 424*e4b17023SJohn Marino template<typename _Tp> 425*e4b17023SJohn Marino inline _Tp __hurwitz_zeta(const _Tp __a,const _Tp __s)426*e4b17023SJohn Marino __hurwitz_zeta(const _Tp __a, const _Tp __s) 427*e4b17023SJohn Marino { 428*e4b17023SJohn Marino return __hurwitz_zeta_glob(__a, __s); 429*e4b17023SJohn Marino } 430*e4b17023SJohn Marino 431*e4b17023SJohn Marino _GLIBCXX_END_NAMESPACE_VERSION 432*e4b17023SJohn Marino } // namespace std::tr1::__detail 433*e4b17023SJohn Marino } 434*e4b17023SJohn Marino } 435*e4b17023SJohn Marino 436*e4b17023SJohn Marino #endif // _GLIBCXX_TR1_RIEMANN_ZETA_TCC 437