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