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

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