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

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