/* version 1.2: prodeulerrat / sumeulerrat: used insufficient precision + let * pin = 2 by default (KB) */ MYDEBUG = 0; /* hurwitz zeta function zeta(x,s). Valid for all reasonable values of x and s, with s not equal to 1. */ { zetahurwitz(x,s)= local(a,res,tes,in,sig,t,m,pr,lim,nn,in2,s1,s2); if (MYDEBUG,gettime); sig=real(x); if (sig>1.5, m=floor(sig-0.5); return (zetahurwitz(x-m,s)-sum(i=1,m,(x-i)^(-s))) ); if (sig<=0, m=ceil(-sig+0.5); return (zetahurwitz(x+m,s)+sum(i=0,m-1,(x+i)^(-s))) ); pr=precision(1.); sig=real(s); t=imag(s); default(realprecision,9); res=s-1.;if(abs(res)<0.1,res=-1,res=log(res)); lim=(pr*log(10)-real((s-.5)*res)+(1.*sig)*log(2.*Pi))/2; lim=max(2,ceil(max(lim,abs(s*1.)/2))); nn=ceil(sqrt((lim+sig/2-.25)^2+(t*1.)^2/4)/Pi); if (MYDEBUG,print("lim, nn: ",[lim,nn])); default(realprecision,pr+5); a=(x+nn+0.)^(-s); res=sum(n=0,nn-1,(x+n)^(-s),a/2); if (MYDEBUG,print("sum from 0 to N-1: ",gettime)); in=x+nn; in2=1./(in*in); s1=2*s-1; s2=s*(s-1); tes=bernreal(2*lim); if (MYDEBUG,print("Bernoullis: ",gettime)); forstep (k=2*lim-2,2,-2, tes=bernreal(k)+in2*(k*k+s1*k+s2)*tes/((k+1)*(k+2)) ); tes=in*(1+in2*s2*tes/2); if (MYDEBUG,print("Bernoulli sum: ",gettime)); res+=tes*a/(s-1); default(realprecision,pr); precision(res,pr); } /* L(chi,s), where chi is a vector of complex values, assumed to be the values from 1 to m of a periodic function. In addition, with zero sum, such as a nontrivial character, if s=1. This is slow if m is large, other implementations should be used. Only given as an illustration of zetahurwitz. */ { L(chi,s)= local(m); m=length(chi); if (s==1, -sum(r=1,m,chi[r]*psi(r/m))/m , sum(r=1,m,chi[r]*zetahurwitz(r/m,s))/m^s ); } /* L(chi_D,s) where chi_D is the quadratic character kronecker(D,n) */ Lquad(D,s)=local(chi);chi=vector(abs(D),i,kronecker(D,i));L(chi,s); /* partial fraction decomposition of rational function F. The answer is a vector of 4-component vectors [al,v,rr,sal], where al ranges through the distinct poles of F, and only one among 2 conjugates if F is real, v is the order of the pole, rr=2 if F is real and al nonreal, and sal is the power series giving the polar decomposition around al, so that if sal=sum(k=0,v-1,ck*x^k)+O(x^v) then around al we have F=sum(k=0,v-1,ck/(x-al)^(v-k))+O(1). */ { ratdec(F)= local(vx,D,N,d,fl,ropro,ro,vfl,al,ct,rr,Dal,sal); vx=variable(F); D=denominator(F); N=numerator(F); d=poldegree(D); fl=(imag(F)==0); ropro=polroots(D); ro=[]; vfl=vector(d,i,1); for (j=1,d, if (vfl[j], al=ropro[j]; ct=1; vfl[j]=0; for (i=j+1,d, if (vfl[i], if (ropro[i]==al, ct++;vfl[i]=0, if (fl && ropro[i]==conj(al), vfl[i]=0) ) ) ); rr=1; if (imag(al)==0,al=real(al),if (fl,rr=2)); Dal=D\((vx-al)^ct); sal=subst(N/Dal,vx,vx+al)+O(vx^ct); ro=concat(ro,[[al,ct,rr,sal]]); ) ); ro; } /* Computes the integral from N to infinity of rational function F. N must be substantially larger than the square of the moduli of the poles of F. */ { intinfrat(F,N,flerr=0)= local(s,vx,pr,lim,G,ro,r); s=0.; if (poldegree(F)>=-1,error("infinite integral in intinfrat")); vx=variable(F); G=subst(F,vx,1/vx); pr=precision(1.); default(realprecision,18); ro=polroots(denominator(G)); r=1;for(k=1,length(ro),r=min(r,norml2(ro[k]))); r=1/r; if (N=1 (default 5) is only used for speeding up the computation. */ { sumrat(F,in=0,la=5)= local(vx,pr,lim,nn,res,sal); if (MYDEBUG,gettime); vx=variable(F); if (in,F=subst(F,vx,vx+in)); if (poldegree(F)>=-1,error("infinite sum in sumrat")); pr=precision(1.); default(realprecision,9); la=max(la+0.,1); lim=ceil((pr*log(10)+2)/(2*(1+log(la)))); nn=ceil((lim+.25)/Pi*la); default(realprecision,pr+5); res=intinfrat(F,nn,1); if (type(res)=="t_VEC", nn=max(nn,ceil(2*res[1])); default(realprecision,9); lim=ceil(nn*Pi/la+.75); default(realprecision,pr+5); res=intinfrat(F,nn); ); if (MYDEBUG,print("lim, nn: ",[lim,nn])); if (MYDEBUG,print("integral: ",gettime)); res+=sum(ijkl=0,nn-1,subst(F,vx,ijkl),subst(F,vx,nn)/2.); if (MYDEBUG,print("sum from 0 to N-1: ",gettime)); sal=1.*subst(F,vx,vx+nn)+O(vx^(2*lim)); if (MYDEBUG,print("power series: ",gettime)); bernreal(2*lim); if (MYDEBUG,print("Bernoulli numbers: ",gettime)); for (k=1,lim, res-=bernreal(2*k)/(2*k)*polcoeff(sal,2*k-1) ); if (MYDEBUG,print("Bernoulli sum: ",gettime)); default(realprecision,pr); precision(res,pr); } /* F rational fraction. Computes the sum from in to infinity of (-1)^n*F(n). By default in=0. Optional parameter la is here set to 3, but only controls the speed. */ { sumratalt(F,in=0,la=3)= local(vx,res); if (type(in)!="t_INT",error("not an integer starting value in sumratalt")); vx=variable(F); if (in,F=subst(F,vx,vx+in)); res=sumrat(2*subst(F,vx,2*vx)-F,0,la); if (poldegree(F)==-1, res-=pollead(numerator(F))/pollead(denominator(F))*log(2) ); if (in%2,res=-res); res; } /* F rational function. Computes the integral from N to infty of log(F). Again N must be sufficiently large. */ { intinflograt(F,N,flerr=0)= local(s,vx,pr,lim,G,ro,r); s=0.; if (poldegree(F-1)>=-1,error("infinite integral in intinflograt")); vx=variable(F); G=subst(F,vx,1/vx); G=G'/G; pr=precision(1.); default(realprecision,18); ro=polroots(denominator(G)); r=1;for(k=1,length(ro),r=min(r,norml2(ro[k]))); r=1/r; if (N=-1,error("infinite product in prodrat")); pr=precision(1.); default(realprecision,9); la=max(la+0.,1); lim=ceil((pr*log(10)+2)/(2*(1+log(la)))); nn=ceil((lim+.25)/Pi*la); default(realprecision,pr+5); tes=intinflograt(F,nn,1); if (type(tes)=="t_VEC", nn=max(nn,ceil(2*tes[1])); default(realprecision,9); lim=ceil(nn*Pi/la+.75); default(realprecision,pr+5); tes=intinflograt(F,nn); ); if (MYDEBUG,print("lim, nn : ",[lim,nn])); if (MYDEBUG,print("integral: ",gettime)); res=prod(ijkl=0,nn-1,subst(F,vx,ijkl),sqrt(subst(F,vx,nn))); if (MYDEBUG,print("product from 0 to N-1: ",gettime)); sal=1.*subst(F'/F,vx,vx+nn)+O(vx^(2*lim-1)); if (MYDEBUG,print("power series: ",gettime)); bernreal(2*lim); if (MYDEBUG,print("Bernoulli numbers: ",gettime)); for (k=1,lim, tes-=bernreal(2*k)*polcoeff(sal,2*k-2)/((2*k)*(2*k-1)) ); if (MYDEBUG,print("Bernoulli sum: ",gettime)); res*=exp(tes); default(realprecision,pr); precision(res,pr); } logzetaa(N,A) = log(zeta(N) * prodeuler(p=2,A,1-1/p^N)); eval_sal(sal, A, pr, lim) = { local(v, m); \\ for efficiency, realprecision had better be low at this point v = abs( Vec(sal * 1.) ); m = 0; for (i = 1, length(v), if (v[i] && v[i] < m, m = v[i])); if (!m, m = 1); default(realprecision, pr + max(0, floor(log(vecmax(v) / m)/log(10)))); sal *= 1.; sum(N = 2, lim, logzetaa(N,A) * sumdiv(N,k,moebius(k)*polcoeff(sal,N/k)/k)); } restore_pr(res, pr) = { default(realprecision, pr); if (precision(res) > pr, precision(res,pr), res); } /* F rational function. Computes the sum over primes from p=pin to infinity of F(p). We must have F(x)=O(1/x^2) as x tends to infty. */ sumeulerrat(F, pin = 2, A = 30) = { local(vx, pr, mro, lim, sal, v, res, m); vx = variable(F); pr = precision(1.); if (poldegree(F) > -2, error("infinite sum in sumeulerrat")); default(realprecision, 9); mro = max(1, vecmax(abs(polroots(denominator(F))))); A = max(pin, ceil(max(A, 3*mro))); lim = 1 + ceil(pr*log(10) / log(A/mro)); sal = subst(F,vx,1/x) + O(x^(lim+1)); res = eval_sal(sal, A, pr, lim); forprime (p = pin, A, res += subst(F,vx,p)); restore_pr(res, pr); } /* F rational function. Computes the product over primes from p=pin to infinity of F(p). We must have F(x)=1+O(1/x^2) as x tends to infty. */ prodeulerrat(F, pin = 2, A = 30) = { local(vx, pr, mro, lim, sal, v, res); vx = variable(F); pr = precision(1.); if (poldegree(F-1) > -2, error("infinite product in prodrat")); default(realprecision, 9); mro = max(1, vecmax(abs(polroots(denominator(F))))); mro = max(mro, vecmax(abs(polroots(numerator(F))))); A = max(pin, ceil(max(A, 3*mro))); lim = 1 + ceil(pr*log(10) / log(A/mro)); sal = log(subst(F,vx,1/x) + O(x^(lim+1))); res = exp( eval_sal(sal, A, pr, lim) ); forprime (p = pin, A, res *= subst(F,vx,p)); restore_pr(res, pr); } /* it is easy if desired in all the above functions involving a rational function F to change so that they apply to F(x^s) instead of F(x). We have not done so. */ /* Examples: \p67 zetahurwitz(1/Pi,exp(-1)) --> 0.1118519821089416198607069686097165736039264753798338741426286156391 chi=[1,-2,1]; L(chi,2) --> 0.6236422637286539664766427555541301381413848710823286510879640179935 L(chi,1) --> 0.3575935377830540795994165103598150138182878451349955619175953008752 Lquad(-7,2) --> 1.151925470544491047101692397320549964797821404686566914083968636166 Lquad(-7,1) --> 1.187410411723725948784625297949363029992334686165035757515202385858 Lquad(5,Pi) --> 0.8688308816459794683875714932200039777264912896134149291409207123151 Lquad(5,1) --> 0.4304089409640040388894332329506054254245706825402896547570061039925 intinfrat(1/(1+x^2),2) --> 0.4636476090008061162142562314612144020285370542861202638109330887201 sumrat(1/(x^3-x-1),2) --> 0.2855375779049113467340954287778781353463574412385555894529422823165 prodrat(1+1/(x^3-x-1),2) --> 1.305705162839975919409358170936919736921608394751995738269456496730 sumeulerrat(1/x^2) --> 0.4522474200410654985065433648322479341732313432398924217364189303511 sumeulerrat(1/(x^3-x-1)) --> 0.2566880688731808665475747518276722430877368417623783159983139346562 prodeulerrat(1+2/x^2) --> 2.190870085553255796350193794718822371567119299935772109133015722465 */