l4/pkg/libstdc++-v3/contrib/libstdc++-v3-4.3.3/include/tr1/riemann_zeta.tcc

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