Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



incompleteGammabased zeta, Li, and Phi numerics
Posted:
Nov 19, 2012 4:55 PM


I have looked into R.E. Crandall's formulae for the numerical evaluation of the Riemann zeta function, the polylogarithm, and the Lerch transcendent, which appeared earlier this year:
<http://www.perfscipress.com/papers/UniversalTOC25.pdf>
These functions are here evaluated for complex orders and arguments via rapidly convergent series in the upper incomplete gamma function. The following machinereadable Derive expressions with corrections of some misprints and misconceptions should facilitate their implementation in computeralgebra systems. The auxiliary functions S(r) := r/(r^2)^(1/2) and sinc(x) := sin(x)/x are used.
1. "Bernoulli series" formula for the Riemann zeta function (eq. 15):
rec_zeta(s,m1:=12,m2:=30,ld:=pi):=1/GAMMA(s)*(SUM(incomplete_ga~ mma(s,ld*n)/n^s,n,1,m1)+ld^(s1)*SUM(ld^n/(n!*(n+s1))*BERNOULLI~ (n),n,0,m2))
The summation limits are preset for about 10 decimal digits accuracy. This expansion beomes impractical for large Im(s) since the terms grow exponentially with Im(s) whereas zeta(s) grows only logarithmically. Note that the Bernoulli numbers BERNOULLI(n) vanish for odd n > 1.
2. "Riemann splitting" formula for the Riemann zeta function (eq. 16, 17):
rec_zeta(s,m:=3):=pi^(s/2)/(2*(s1)*GAMMA(s/2+1))+1/GAMMA(s/2)*~ (SUM(incomplete_gamma(s/2,pi*n^2)/n^s,n,1,m)+pi^(s1/2)*SUM(inco~ mplete_gamma((1s)/2,pi*u^2)/u^(1s),u,1,m))
Here too, the summation limits are preset for about 10 digits accuracy, and the expansion beomes impractical for large Im(s). Both remarks apply to the other functions too, and will not be repeated.
3. "Bernoulli series" formula for the Hurwitz zeta function (eq. 29):
rec_zeta(s,a,m1:=12,m2:=30,ld:=pi):=1/GAMMA(s)*(SUM(incomplete_~ gamma(s,ld*(n+a))/(n+a)^s,n,0,m1)+ld^(s1)*SUM((ld)^n/(n!*(n+s~ 1))*MY_BERNOULLI_POLY(n,a),n,0,m2))
This expansion holds for all a /= 0, 1, 2, ... (My redefinition of BERNOULLI_POLY(n,x) removes a minor problem with the implementation of the Bernoulli polynomials in Derive.)
4. "Riemann splitting" formula for the Hurwitz zeta function (eq. 30):
rec_zeta(s,a,m:=3):=pi^(s/2)/(s1)*IF(a=1,1/(2*GAMMA(s/2+1)),1/~ GAMMA(s/2))+1/2*SUM(1/((n+a)^2)^(s/2)*(incomplete_gamma(s/2,pi*(~ n+a)^2)/GAMMA(s/2)+incomplete_gamma((s+1)/2,pi*(n+a)^2)/GAMMA((s~ +1)/2)*S(n+a)),n,SELECT(n+a/=0,n,m1,m))+pi^(s1/2)*SUM(1/u^(1~ s)*(incomplete_gamma((1s)/2,pi*u^2)/GAMMA(s/2)*COS(2*pi*a*u)+in~ complete_gamma(1s/2,pi*u^2)/GAMMA((s+1)/2)*SIN(2*pi*a*u)),u,1,m)
Here, pi^s has been corrected to pi^(s/2). For 0 < Re(a) <= 1 if Im(a) <= 0, and for 0 <= Re(a) < 1 if Im(a) > 0, this formula returns the customary Hurwitz zeta function. Application of the functional equation zeta(s,a) = zeta(s,a+1) + a^(s) then allows to compute zeta(s,a) for any other value of a. The Riemann zeta function is recovered for a=1: zeta(s,1) = zeta(s).
5. "Bernoulli series" formula for the polylogarithm (eq. 34):
rec_li(s,z,m1:=12,m2:=30,ld:=pi):=(2*pi)^(s1)*(#i^(1s)*SUM(inc~ omplete_gamma(1s,ld*(n+1/2+LN(z)/(2*pi*#i)))/(n+1/2+LN(z)/(2*~ pi*#i))^(1s),n,0,m1)+#i^(s1)*SUM(incomplete_gamma(1s,ld*(n+1/~ 2LN(z)/(2*pi*#i)))/(n+1/2LN(z)/(2*pi*#i))^(1s),n,0,m1)pi/l~ d^s*SUM((#i*ld)^n/n!*sinc(pi/2*(ns))*MY_BERNOULLI_POLY(n,1/2L~ N(z)/(2*pi*#i)),n,0,m2))
Here, the omission of a prefactor (2 pi)^(s1) has been corrected, and the polylogarithm argument has been converted from exp(x) to z. This expansion holds for any z /= 1, for z=1 one has to substitute Li(s,1) = zeta(s). Note that Li(s,z) has a pole at z=1 if Re(s) <= 1.
6. "Riemann splitting" formula for the polylogarithm (eq. 37):
rec_li(s,z,m:=3):=pi^(s/2)/(2*GAMMA(s/2+1))*IF(z=1,1/(s1),1)+S~ UM(1/n^s*(incomplete_gamma(s/2,pi*n^2)/GAMMA(s/2)*COSH(n*LN(z))+~ incomplete_gamma((s+1)/2,pi*n^2)/GAMMA((s+1)/2)*SINH(n*LN(z))),n~ ,1,m)+pi^(s1/2)/2*SUM(1/((u+LN(z)/(2*pi*#i))^2)^((1s)/2)*(inco~ mplete_gamma((1s)/2,pi*(u+LN(z)/(2*pi*#i))^2)/GAMMA(s/2)+#i*(in~ complete_gamma(1s/2,pi*(u+LN(z)/(2*pi*#i))^2)/GAMMA((s+1)/2))*S~ (u+LN(z)/(2*pi*#i))),u,SELECT(u+LN(z)/(2*pi*#i)/=0,u,m,m))
Here, the factors z^n and exp( 2 pi i a u) have been corrected to 1. This expansion hold for all complex z, although it returns Li(s,1) = zeta(s) for Re(s) <= 1 too. Note that 2 cosh(n ln(z)) = z^n + z^(n) and 2 sinh(n ln(z)) = z^n  z^(n).
7. "Riemann splitting" formula for the Lerch transcendent (eq. 14):
rec_phi(z,s,a,m:=3):=pi^(s/2)*(IF(z=1,1/((s1)*GAMMA(s/2)),0)~ IF(a=1,1/(2*z*GAMMA(s/2+1)),0))+1/2*SUM(z^n/((n+a)^2)^(s/2)*(inc~ omplete_gamma(s/2,pi*(n+a)^2)/GAMMA(s/2)+incomplete_gamma((s+1)/~ 2,pi*(n+a)^2)/GAMMA((s+1)/2)*S(n+a)),n,SELECT(n+a/=0,n,m1,m))+~ pi^(s1/2)/2*SUM(EXP((2*pi*#i*u+LN(z))*a)/((u+LN(z)/(2*pi*#i))^~ 2)^((1s)/2)*(incomplete_gamma((1s)/2,pi*(u+LN(z)/(2*pi*#i))^2)~ /GAMMA(s/2)+#i*(incomplete_gamma(1s/2,pi*(u+LN(z)/(2*pi*#i))^2)~ /GAMMA((s+1)/2))*S(u+LN(z)/(2*pi*#i))),u,SELECT(u+LN(z)/(2*pi*#i~ )/=0,u,m,m))
For 0 < Re(a) <= 1 if Im(a) <= 0, and for 0 <= Re(a) < 1 if Im(a) > 0, this formula returns the customary Lerch transcendent (not some variant). Application of the functional equation Phi(z,s,a) = z Phi(z,s,a+1) + a^(s) then allows to compute Phi(z,s,a) for any other value of a. The polylogarithm is recovered for a=1: z Phi(z,s,1) = Li(s,z). The expansion returns Phi(1,s,a) = zeta(s,a) for Re(s) <= 1 too, whereas Phi(z,s,a) has a pole at z=1 if Re(s) <= 1. Note that exp( (2 pi i u + ln(z)) a) = exp( 2 pi i u a) / z^a. The functional equation (7) for the Lerch transcendent is not entirely correct, the corrected equation was posted earlier this year in the thread "Inversion Lerch Phi".

Because the Derive kernel doesn't implement the upper incomplete gamma function, I had to define it in order to check the above expansions:
gamma_fraction1(s,z,m:=20):=GAMMA(s)z^s*EXP(z)/(ITERATE([sz+i~ _*z/(q_+i_),i_1],[q_,i_],[s,m],m)) SUB 1
gamma_fraction2(s,z,m:=20):=z^s*EXP(z)/(ITERATE([z(si_)*q_/(q~ _+i_),i_1],[q_,i_],[z,m],m)) SUB 1
incomplete_gamma(s,z,m:=30):=IF([23*RE(z),6.5*IM(z)]^2<m^2 OR (R~ E(z)<0 AND 6.5*ABS(IM(z))<m),gamma_fraction1(s,z,m),gamma_fracti~ on2(s,z,m))
This delivers 10 decimal digits for arbitrary complex z, but only as long as s stays in the complex vicinity of s = 1/2. So this simple code is in need of refinement: Nothing is done about the cancellation of diverging terms near nonpositive integer s, and the fixed evaluation depth of the continued fractions produces inaccurate results if z happens to be close to a pole of the mth convergent: the pole loci form chains that sweep across the zplane as s and m vary, whereas my switching between the continued fractions is ignorant of s.
The series expansion of Gamma(s,z) in terms of generalized Laguerre polynomials (eq. 25) simply reproduces the convergents of Legendre's continued fraction (eq. 24):
gamma_fraction2(s,z,m) = = (s1)*gamma_series4(s1,z,m) + z^(s1)*exp(z)
where
gamma_series4(s,z,m:=20,ll):=PROG(ll:=VECTOR(MY_GENERALIZED_LAGU~ ERRE(n_,s,z),n_,0,m),z^s*EXP(z)*SUM(PRODUCT(p_s,p_,1,n_1)/(~ n_!*ll SUB n_*ll SUB (n_+1)),n_,1,m))
(My redefinition of GENERALIZED_LAGUERRE(n,a,x) removes a minor problem with the implementation of the generalized Laguerre polynomials in Derive.) Thus, the series can be used to analyze the convergence behavior of the continued fraction, but it makes no sense at all to compute Gamma(s,z) in this way.
Convergents of the complementary continued fraction obey a very similar relation:
gamma_fraction1(s,z,m) = 1/s*(gamma_series3(s+1,z,m)  z^s*exp(z))
where
gamma_series3(s,z,m:=20,ll):=PROG(ll:=VECTOR(MY_GENERALIZED_LAGU~ ERRE(n_,sn_,z),n_,0,m),GAMMA(s)+z^s*EXP(z)*SUM((z)^(n_1)/(~ n_!*ll SUB n_*ll SUB (n_+1)),n_,1,m))
How do the free computeralgebra systems compute Gamma(10^100,10), say, to ten decimal digits? Do they try to evaluate continued fractions only away from denominator roots?
Martin.



