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