/************* ComputeL v1.3.6, 2001-2018, (c) Tim Dokchitser ************/
/**************** computing special values of L-functions ****************/
/* arXiv.org/abs/math.NT/0207280, Exper. Math. 13 (2004), no. 2, 137-150 */
/****** Questions/comments welcome! -> tim.dokchitser@bristol.ac.uk ******/
\\ ACKNOWLEDGEMENTS: I'd like to thank Mark Watkins, Steve Donnelly,
\\ William Stein, Anton Mellit, Almasa Odzak, Karim Belabas, Myoungil Kim,
\\ Chris King, F. Patrick Rabarison, Neil Dummigan, Maciej Radziejewski,
\\ François Brunault and Alex Best, for examples, bug fixes and suggestions
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Distributed under the terms of the GNU General Public License (GPL)
\\ This code is distributed in the hope that it will be useful,
\\ but WITHOUT ANY WARRANTY; without even the implied warranty of
\\ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
\\ GNU General Public License for more details.
\\ The full text of the GPL is available at:
\\ http://www.gnu.org/licenses/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ USAGE: Take an L-function L(s) = sum of a(n)/n^s over complex numbers
\\ e.g. Riemann zeta-function, Dedekind zeta-function,
\\ Dirichlet L-function of a character, L-function
\\ of a curve over a number field, L-function of a modular form,
\\ any ``motivic'' L-function, Shintani's zeta-function etc.
\\ assuming L(s) satisfies a functional equation of a standard form,
\\ this package computes L(s) or its k-th derivative for some k
\\ for a given complex s to required precision
\\ - a short usage guide is provided below
\\ - or (better) just look at the example files ex-*
\\ they are hopefully self-explanatory
\\
\\ ASSUMED: L^*(s) = Gamma-factor * L(s) satisfies functional equation
\\ L^*(s) = sgn * L^*(weight-s),
\\ [ more generally L^*(s) = sgn * Ldual^*(weight-s) ]
\\ Gamma-factor = A^s * product of Gamma((s+gammaV[k])/2)
\\ where A = sqrt(conductor/Pi^d)
\\
\\ gammaV = list of Gamma-factor parameters,
\\ e.g. [0] for Riemann zeta, [0,1] for ell.curves
\\ conductor = exponential factor (real>0, usually integer),
\\ e.g. 1 for Riemann zeta and modular forms under SL_2(Z)
\\ e.g. |discriminant| for number fields
\\ e.g. conductor for H^1 of curves/Q
\\ weight = real > 0 (usually integer, =1 by default)
\\ e.g. 1 for Riemann zeta, 2 for H^1 of curves/Q
\\ sgn = complex number (=1 by default)
\\
\\ 1. Read the package (\rcomputel)
\\ 2. Set the required working precision (say \p28)
\\
\\ 3. DEFINE gammaV, conductor, weight, sgn,
\\ Lpoles = vector of points where L^*(s) has (simple) poles
\\ Only poles with Re(s)>weight/2 are to be included
\\ Lresidues = vector of residues of L^*(s) in those poles
\\ or set Lresidues = automatic (default value; see ex-nf)
\\ if necessary, re-define coefgrow(), MaxImaginaryPart (see below)
\\
\\ [4.] CALL cflength() determine how many coefficients a(n) are necessary
\\ [optional] to perform L-function computations
\\
\\ 5. CALL initLdata(cfstr) where cfstr (e.g. "(-1)^k") is a string which
\\ evaluates to k-th coefficient a(k) in L-series, e.g.
\\ N = cflength(); \\ say returns N=10
\\ Avec = [1,-1,0,1,-1,0,1,-1,0,1,-1,0]; \\ must be at least 10 long
\\ initLdata("Avec[k]");
\\ If Ldual(s)<>L(s), in other words, if the functional equation involves
\\ another L-function, its coefficients are passed as a 3rd parameter,
\\ initLdata("Avec[k]",,"conj(Avec[k])"); see ex-chgen as an example
\\
\\ [7.] CALL checkfeq() verify how well numerically the functional
\\ [optional] equation is satisfied
\\ also determines the residues if Lpoles!=[]
\\ and Lresidues=automatic
\\ More specifically: for T>1 (default 1.2), checkfeq(T) should ideally
\\ return 0 (with current precision, e.g. 3.2341E-29 for \p28 is good)
\\ * if what checkfeq() returns does not look like 0 at all,
\\ probably functional equation is wrong
\\ (i.e. some of the parameters gammaV, conductor etc., or the coeffs)
\\ * if checkfeq(T) is to be used, more coefficients have to be
\\ generated (approximately T times more), e.g. call
\\ cflength(1.3), initLdata("a(k)",1.3), checkfeq(1.3)
\\ * T=1 always (!) returns 0, so T has to be away from 1
\\ * default value T=1.2 seems to give a reasonable balance
\\ * if you don't have to verify the functional equation or the L-values,
\\ call cflength(1) and initLdata("a(k)",1),
\\ you need slightly less coefficients then
\\
\\ 8. CALL L(s0) to determine the value of L-function L(s) in s=s0
\\ CALL L(s0,c) with c>1 to get the same value with a different cutoff
\\ point (c close to 1); should return the same answer,
\\ good to check if everything works with right precision
\\ (if it doesn't, email me!)
\\ needs generally more coefficients for larger ex
\\ if L(s0,ex)-L(s0) is large, either the functional eq.
\\ is wrong or loss of precision (should get a warning)
\\ CALL L(s0,,k) to determine k-th derivative of L(s) in s=s0
\\ see ex-bsw for example
\\ CALL Lseries(s,,k) to get first k terms of Taylor series expansion
\\ L(s)+L'(s)S+L''(s)*S^2/2!+...
\\ faster than k calls to L(s)
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Default values for the L-function parameters \\
\\ All may be (and conductor and gammaV must be) re-defined \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ MUST be re-defined, gives error if unchanged
conductor = automatic;
gammaV = automatic;
\\ MAY be re-defined
weight = 1; \\ by default L(s)<->L(1-s)
sgn = 1; \\ with sign=1 in functional equation
Lpoles = []; \\ and L*(s) has no poles
Lresidues = automatic; \\ if this is not changed to [r1,r2,...] by hand,
\\ checkfeq() tries to determine residues automatically
\\ see ex-nf for instance
{
coefgrow(n) = if(length(Lpoles), \\ default bound for coeffs. a(n)
1.5*n^(vecmax(real(Lpoles))-1), \\ you may redefine coefgrow() by hand
sqrt(4*n)^(weight-1)); \\ if your a(n) have different growth
} \\ see ex-delta for example
\\ - For s with large imaginary part there is a lot of cancellation when
\\ computing L(s), so a precision loss occurs. You then get a warning message
\\ - If you want to compute L(s), say, for s=1/2+100*I,
\\ set MaxImaginaryPart=100 before calling initLdata()
\\ - global variable PrecisionLoss holds the number of digits lost in
\\ the last calculation (indepedently of the MaxImaginaryPart setting)
MaxImaginaryPart = 0; \\ re-define this if you want to compute L(s)
\\ for large imaginary s (see ex-zeta2 for example)
MaxAsympCoeffs = 40; \\ At most this number of terms is generated
\\ in asymptotic series for phi(t) and G(s,t)
\\ default value of 40 seems to work generally well
/******************* IMPLEMENTATION OF THE PACKAGE ************************/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Some helfpul functions \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Extraction operations on vectors
vecleft(v,n) = vecextract(v,concat("1..",Str(n)));
vecright(v,n) = vecextract(v,concat(Str(length(v)-n+1),concat("..",Str(length(v)))));
\\ Tabulate a string to n characters, e.g. StrTab(3,2)="3 ";
StrTab(x,n) = x=Str(x);while(length(x)gamma
if (reflect, \\ reflect if necessary
sinser=Vec(sin(Pi*z));
if (negint,sinser[1]=0); \\ taking slight care at integers<0
res=subst(Pi/res/Ser(sinser),x,-x);
);
default(realprecision,digts);
)))));
default(seriesprecision,srprec);
res;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ fullgamma(ss) - the full gamma factor (at s=ss) \\
\\ vA^s*Gamma((s+gammaV[1])/2)*...*Gamma((s+gammaV[d])/2) \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
fullgamma(ss) =
if(ss!=lastFGs,lastFGs=ss;
lastFGval=prod(j=1,length(gammaV),gamma((ss+gammaV[j])/2),vA^ss)
);
lastFGval;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ fullgammaseries(ss,extraterms) - Laurent series for the gamma factor \\
\\ without the exponential factor, i.e. \\
\\ Gamma((s+gammaV[1])/2)*...*Gamma((s+gammaV[d])/2) \\
\\ around s=ss with a given number of extra terms. The series variable is S.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
fullgammaseries(ss,extraterms)=
local(digts,GSD);
digts=default(realprecision);
if (lastFGSs!=ss || lastFGSterms!=extraterms,
GSD=sum(j=1,numpoles,(abs((ss+poles[j])/2-round(real((ss+poles[j])/2)))<10^(2-digts)) * PoleOrders[j] )+extraterms;
lastFGSs=ss;
lastFGSterms=extraterms;
lastFGSval=subst(prod(j=1,length(gammaV),gammaseries((ss+gammaV[j])/2,GSD)),x,S/2);
);
lastFGSval;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ RecursionsAtInfinity(gammaV) \\
\\ Recursions for the asymptotic expansion coefficients \\
\\ for phi(x) and G(s,x) at infinity. \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
RecursionsAtInfinity(gammaV)=
local(d,p,j,k,symvec,modsymvec,deltapol,recF,recG);
\\ d = number of Gamma-factors in question
\\ gammaV[k] = Gamma-factors
\\ symvec = vector of elementary symmetric functions
\\ 1, gammaV[1]+...+gammaV[d], ... , gammaV[1]*...*gammaV[d], 0
\\ modsymvec = symmetric expressions used in the formula
d = length(gammaV);
symvec = concat(Vec(prod(k=1,d,(x+gammaV[k]))),[0]);
modsymvec = vector(d+2,k,0);
for (j=0,d,
for (k=0,j,
modsymvec[j+1]+=(-symvec[2])^k*d^(j-1-k)*binomial(k+d-j,k)*symvec[j-k+1];
));
\\ Delta polynomials
OldSeriesPrecision = default(seriesprecision);
default(seriesprecision,2*d+2);
deltapol=subst(Vec( (sinh(x)/x)^tvar ),tvar,x);
default(seriesprecision,OldSeriesPrecision);
\\ recursion coefficients for phi at infinity
recF=vector(d+1,p,
-1/2^p/d^(p-1)/n*sum(m=0,p,modsymvec[m+1]*prod(j=m,p-1,d-j)*
sum(k=0,floor((p-m)/2),(2*n-p+1)^(p-m-2*k)/(p-m-2*k)!*subst(deltapol[2*k+1],x,d-p))));
\\ recursion coefficients for G at infinity
recG=vector(d,p,recF[p+1]-(symvec[2]+d*(s-1)-2*(n-p)-1)/2/d*recF[p]);
[vector(d-1,p,recF[p+1]),recG] \\ return them both
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SeriesToContFrac(vec) \\
\\ Convert a power series vec[1]+vec[2]*x+vec[3]*x^2+... \\
\\ into a continued fraction expansion. \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
SeriesToContFrac(vec)=
local(res=[],ind);
vec=1.*vec;
while (1,
res=concat(res,[vec[1]]);
ind=0;
until(ind==length(vec) || abs(vec[ind+1])>10^-asympdigits,ind++;vec[ind]=0);
if(ind>=length(vec),break);
res=concat(res,[ind]);
vec=Vec(x^ind/Ser(vec));
);
res;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ EvaluateContFrac(cfvec,terms,t) \\
\\ Evaluate a continued fraction at x=t, using a given number of terms \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
EvaluateContFrac(cfvec,terms,t)=
local(res);
if (terms<0 || terms>length(cfvec)\2,terms=length(cfvec)\2);
res=cfvec[2*terms+1];
while(terms>0,res=if(res,cfvec[2*terms-1]+t^cfvec[2*terms]/res,10^asympdigits);terms--);
res;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ cflength( cutoff=1.2 ) \\
\\ \\
\\ number of coefficients necessary to work with L-series with \\
\\ current Gamma-factor gammaV, conductor, weight and working precision \\
\\ - CUTOFF specifies largest t used as a cutoff point in checkfeq(t) \\
\\ and L(...,t,...). Default is 1.2. Set it to 1 if checkfeq() \\
\\ is not to be used. \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
cflength(cutoff=1.2)=
local(vA,d,expdifff,asympconstf,err,t,t1=0,t2=2,tt,res,lfundigits);
vA = sqrt(conductor/Pi^length(gammaV));
d = length(gammaV);
lfundigits = default(realprecision) +
max(ceil(-log(abs(fullgamma(0.7*weight+MaxImaginaryPart*I)))/log(10)),0);
expdifff = (sum(k=1,d,gammaV[k])+1)/d-1;
asympconstf = 2*prod(k=1,d,gamma(k/d));
err = 10^(-lfundigits-0.5);
until (t2-t1<=1,
t = if(t1,(t1+t2)\2,t2);
tt = t/cutoff/vA;
res = coefgrow(t) * asympconstf*exp(-d*tt^(2/d))*tt^expdifff;
if (t1,if(abs(res)>err,t1=t,t2=t),if(abs(res)1)&&(terms>1)&&
(abs(phiinf(t1,terms-1)-phiinf(t1,terms)))<10^(-termdigits-1),terms-=1);
if (sum(j=1,terms,fcf[2*j]) phiVser
phiV = [];
phiVnn = 0;
phiVser = InitV;
until((phiVnn>3)&&(vecmax(abs(phiV[phiVnn]))*((PhiCaseBound+1)*vA)^(2*phiVnn)<10^(-termdigits-1)),
RecursephiV());
default(realprecision,answerdigits);
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ phi(t) - inverse Mellin transform of the product of Gamma-factors \\
\\ computed either with Taylor in 0 or from asymptotics \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
phi(t)=if(t3),
totalold=res;
nn++;
if(nn>phiVnn,RecursephiV());
res+=TPower*phiV[nn]*LogTTerm;
TPower*=t2;
);
default(realprecision,termdigits);
res;
}
{ \\ phi(t) using asymptotic expansion at infinity
phiinf(t,ncf=fncf)=
local(res,d,td2);
default(realprecision,asympdigits);
t=precision(t,asympdigits);
d=length(gammaV);
td2=t^(-2/d);
res=EvaluateContFrac(fcf,ncf,td2);
res=res*asympconstf*exp(-d/td2)*t^expdifff;
default(realprecision,termdigits);
res;
}
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ G(t,s) - incomplete Mellin transform of phi(t) divided by x^s \\
\\ computed either with Taylor in 0 or from asymptotics \\
\\ G(t,s,k) - its k-th derivative (default 0) \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{
G(t,ss,der=0)=
local(nn);
if(lastGs!=[ss,der] || type(Ginfterms)!="t_VEC",
initGinf(ss,der);lastGs=[ss,der]);
if(t3) && abs(term)<10^-(termdigits+1),
nn++;
if(nn>length(LogSum),MakeLogSum(ss,der));
term=TPower*subst(LogSum[nn],lt,LT);
res+=term;
TPower*=t2;
);
default(realprecision,lfundigits);
gmser=fullgammaseries(ss,der)/t^(S+ss);
gmcf=polcoeff(gmser,der,S)*der!;
res=(gmcf-res);
default(realprecision,termdigits);
res;
}
{
Ginf(t,ss,der,ncf=-1)= \\ G(t,s,der) computed using asymptotic expansion
local(res,d,tt); \\ at infinity and associated continued fraction
default(realprecision,asympdigits);
ss=precision(ss,asympdigits);
t=precision(t,asympdigits);
if (ncf==-1,ncf=gncf);
d=length(gammaV);
tt=t^(-2/d);
res=EvaluateContFrac(gcf,ncf,tt);
res=asympconstg*exp(-d/tt)*t^expdiffg*tt^der*res;
default(realprecision,termdigits);
res;
}
{
initGinf(ss,der)= \\ pre-compute asymptotic expansions for a given s
local(d,gvec,gncf,terms,t1);
default(realprecision,asympdigits);
ss=precision(ss,asympdigits);
d=length(gammaV);
gvec=Gvec;
for (k=1,der,gvec=deriv(gvec,s);gvec=concat(vecright(gvec,length(gvec)-1),1));
gcf=SeriesToContFrac(vector(ncoeff+1,k,subst(gvec[d+k-1],s,ss)));
gncf=length(gcf)\2;
Ginfterms=vector(round(lastt/termstep)+1,k,-1);
terms=gncf;
GCaseBound=0;
for (k=1,length(Ginfterms),
t1=(k-1)*termstep;
while ((k>1)&&(terms>1)&&
(abs(Ginf(t1,ss,der,terms-1)-Ginf(t1,ss,der,terms)))<10^(-termdigits-2),terms-=1);
if (sum(j=1,terms,gcf[2*j])=1 in general \\
\\ must be equal to 1 if k>0 \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
L(ss,cutoff=1,der=0)=polcoeff(Lseries(ss,cutoff,der),der)*der!
{ \\ Lseries(s,,der) = L(s)+L'(s)S+L''(s)*S^2/2!+... ,first der terms
Lseries(ss,cutoff=1,der=0)=
local(FGSeries,LGSeries,res);
default(realprecision,lfundigits);
FGSeries = fullgammaseries(ss,der)*vA^(ss+S);
if (length(Lpoles) && (vecmin(abs(vector(length(Lpoles),k,Lpoles[k]-ss)))<10^(-answerdigits) ||
vecmin(abs(vector(length(Lpoles),k,weight-Lpoles[k]-ss)))<10^(-answerdigits)),
error("L*(s) has a pole at s=",ss));
\\ if (vecmin(abs(vector(length(Lpoles),k,Lpoles[k]-ss)))<10^(-answerdigits) ||
\\ vecmin(abs(vector(length(Lpoles),k,weight-Lpoles[k]-ss)))<10^(-answerdigits),
\\ error("L*(s) has a pole at s=",ss));
PrecisionLoss = ceil(-log(vecmax(abs(Vec(FGSeries))))/log(10))
-(lfundigits-answerdigits);
if (PrecisionLoss>1,
print("Warning: Loss of ",PrecisionLoss," digits due to cancellation"));
LSSeries = sum(k=0,der,Lstar(ss,cutoff,k)*S^k/k!)+O(S^(der+1));
res=LSSeries/FGSeries;
default(realprecision,answerdigits);
res;
}
{
Lstar(ss,cutoff=1,der=0)= \\ Lstar(s) = L(s) * Gamma factor
local(res,ncf1,ncf2);
if (der & (cutoff!=1),error("L(s,cutoff,k>0) is only implemented for cutoff=1"));
ss = precision(ss,taylordigits);
cutoff = precision(cutoff,taylordigits);
ncf1 = min(round(lastt*vA*cutoff),length(cfvec));
ncf2 = min(round(lastt*vA/cutoff),length(cfvec));
default(realprecision,termdigits);
res=(-sum(k=1,length(Lpoles),(-1)^der*der!*Lresidues[k]/(ss-Lpoles[k])^(der+1)*cutoff^(-Lpoles[k]))
-sgn*sum(k=1,length(Lpoles),Lresidues[k]*der!/(weight-Lpoles[k]-ss)^(der+1)*cutoff^(-weight+Lpoles[k]))
+sgn*sum(k=1,ncf1,if(cfdualvec[k],cfdualvec[k]*(-1)^der*G(k/vA/cutoff,weight-ss,der),0)/cutoff^weight)
+sum(k=1,ncf2,if(cfvec[k],cfvec[k]*G(k*cutoff/vA,ss,der),0))
)*cutoff^ss;
res;
}