libstdc++-v3/include/tr1/riemann_zeta.tcc

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