```Date: Nov 19, 2012 4:55 PM
Author: clicliclic@freenet.de
Subject: incomplete-Gamma-based zeta, Li, and Phi numerics

I have looked into R.E. Crandall's formulae for the numerical evaluationof the Riemann zeta function, the polylogarithm, and the Lerchtranscendent, which appeared earlier this year:  <http://www.perfscipress.com/papers/UniversalTOC25.pdf>These functions are here evaluated for complex orders and arguments viarapidly convergent series in the upper incomplete gamma function. Thefollowing machine-readable Derive expressions with corrections of somemisprints and misconceptions should facilitate their implementation incomputer-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 growexponentially 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 applyto 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 ofBERNOULLI_POLY(n,x) removes a minor problem with the implementation ofthe 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 thecustomary Hurwitz zeta function. Application of the functional equationzeta(s,a) = zeta(s,a+1) + a^(-s) then allows to compute zeta(s,a) forany 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, andthe polylogarithm argument has been converted from exp(-x) to z. Thisexpansion 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) and2 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 somevariant). Application of the functional equation Phi(z,s,a) = zPhi(z,s,a+1) + a^(-s) then allows to compute Phi(z,s,a) for any othervalue 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) <= 1too, 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 correctedequation was posted earlier this year in the thread "Inversion LerchPhi".-Because the Derive kernel doesn't implement the upper incomplete gammafunction, 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 1gamma_fraction2(s,z,m:=20):=z^s*EXP(-z)/(ITERATE([z-(s-i_)*q_/(q~_+i_),i_-1],[q_,i_],[z,m],m)) SUB 1incomplete_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 aslong as s stays in the complex vicinity of s = 1/2. So this simple codeis in need of refinement: Nothing is done about the cancellation ofdiverging terms near nonpositive integer s, and the fixed evaluationdepth of the continued fractions produces inaccurate results if zhappens to be close to a pole of the m-th convergent: the pole loci formchains that sweep across the z-plane as s and m vary, whereas myswitching between the continued fractions is ignorant of s.The series expansion of Gamma(s,z) in terms of generalized Laguerrepolynomials (eq. 25) simply reproduces the convergents of Legendre'scontinued fraction (eq. 24):  gamma_fraction2(s,z,m) =  = (s-1)*gamma_series4(s-1,z,m) + z^(s-1)*exp(-z)wheregamma_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 problemwith the implementation of the generalized Laguerre polynomials inDerive.) Thus, the series can be used to analyze the convergencebehavior of the continued fraction, but it makes no sense at all tocompute Gamma(s,z) in this way.Convergents of the complementary continued fraction obey a very similarrelation:  gamma_fraction1(s,z,m) = 1/s*(gamma_series3(s+1,z,m) - z^s*exp(-z))wheregamma_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 fractionsonly away from denominator roots?Martin.
```