Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.
|
|
|
|
incomplete-Gamma-based 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 machine-readable Derive expressions with corrections of some
misprints and misconceptions should facilitate their implementation in
computer-algebra 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^(s-1)*SUM(ld^n/(n!*(n+s-1))*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*(s-1)*GAMMA(s/2+1))+1/GAMMA(s/2)*~
(SUM(incomplete_gamma(s/2,pi*n^2)/n^s,n,1,m)+pi^(s-1/2)*SUM(inco~
mplete_gamma((1-s)/2,pi*u^2)/u^(1-s),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^(s-1)*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)/(s-1)*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,-m-1,m))+pi^(s-1/2)*SUM(1/u^(1-~
s)*(incomplete_gamma((1-s)/2,pi*u^2)/GAMMA(s/2)*COS(2*pi*a*u)+in~
complete_gamma(1-s/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)^(s-1)*(#i^(1-s)*SUM(inc~
omplete_gamma(1-s,ld*(n+1/2+LN(-z)/(2*pi*#i)))/(n+1/2+LN(-z)/(2*~
pi*#i))^(1-s),n,0,m1)+#i^(s-1)*SUM(incomplete_gamma(1-s,ld*(n+1/~
2-LN(-z)/(2*pi*#i)))/(n+1/2-LN(-z)/(2*pi*#i))^(1-s),n,0,m1)-pi/l~
d^s*SUM((-#i*ld)^n/n!*sinc(pi/2*(n-s))*MY_BERNOULLI_POLY(n,1/2-L~
N(-z)/(2*pi*#i)),n,0,m2))
Here, the omission of a prefactor (2 pi)^(s-1) 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/(s-1),-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^(s-1/2)/2*SUM(1/((u+LN(z)/(2*pi*#i))^2)^((1-s)/2)*(inco~
mplete_gamma((1-s)/2,pi*(u+LN(z)/(2*pi*#i))^2)/GAMMA(s/2)+#i*(in~
complete_gamma(1-s/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/((s-1)*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,-m-1,m))+~
pi^(s-1/2)/2*SUM(EXP(-(2*pi*#i*u+LN(z))*a)/((u+LN(z)/(2*pi*#i))^~
2)^((1-s)/2)*(incomplete_gamma((1-s)/2,pi*(u+LN(z)/(2*pi*#i))^2)~
/GAMMA(s/2)+#i*(incomplete_gamma(1-s/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([s-z+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-(s-i_)*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 m-th convergent: the pole loci form
chains that sweep across the z-plane 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) =
= (s-1)*gamma_series4(s-1,z,m) + z^(s-1)*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_,-s-n_,-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 computer-algebra systems compute Gamma(10^-100,-10),
say, to ten decimal digits? Do they try to evaluate continued fractions
only away from denominator roots?
Martin.
|
|
|
|
