/*********** ComputeL v1.3, Apr 2006, (c) Tim Dokchitser ***********/
/************ (computing special values of L-functions) ************/
/******* see preprint http://arXiv.org/abs/math.NT/0207280, ********/
/******* appeared in Exper. Math. 13 (2004), no. 2, 137-150 ********/
/*** Questions/comments welcome! -> t.dokchitser@dpmms.cam.ac.uk ***/

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ 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)<n,x=concat(x," "));x

\\ Concatenate up to 4 strings
concatstr(s1="",s2="",s3="",s4="")=concat(Str(s1),concat(Str(s2),concat(Str(s3),Str(s4))))

\\ Print a ``small error'', e.g. 0.00000013 as "1E-7"
{
errprint(x)=if(type(x)=="t_COMPLEX",x=abs(x));
   if(x==0,concatstr("1E-",default(realprecision,,1)+1),
   concatstr(truncate(x/10^floor(log(abs(x))/log(10))),"E",floor(log(abs(x))/log(10))));
}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ gammaseries(z0,terms)                                                   \\
\\ Taylor series expansion of Gamma(z0+x) around 0, z0 arbitrary complex   \\
\\ - up to O(x^(terms+1))                                                  \\
\\ - uses current real precision                                           \\
\\ See Luke "Mathematical functions and their approximations", section 1.4  \
\\ note a misprint there in the recursion formulas [(z-n) term in c3 below] \
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
gammaseries(z0,terms,
    Avec,Bvec,Qvec,n,z,err,res,c0,c1,c2,c3,sinser,reflect,digits,srprec,negint)=
  srprec=default(seriesprecision,,1);
  if (z0==real(round(z0)),z0=real(round(z0)));    \\ you don't want to know
  negint=type(z0)=="t_INT" && z0<=0;              \\ z0 is a pole
  default(seriesprecision,terms+1+negint,1);
  if (terms==0 && !negint,res=gamma(z0)+O(x),     \\ faster to use
  if (terms==1 && !imag(z0) && !negint,           \\   built-in functions
      res=gamma(z0)*(1+psi(z0)*x+O(x^2)),         \\   in special cases
  if (z0==0, res=gamma(1+x)/x,
  if (z0==1, res=gamma(1+x),
  if (z0==2, res=gamma(1+x)*(1+x),
  \\ otherwise use Luke's rational approximations for psi(x)
  digits=default(realprecision,,1);   \\ save working precision
  default(realprecision,digits+3,1);  \\   and work with 3 digits more
  reflect=real(z0)<0.5;               \\ left of 1/2 use reflection formula
  if (reflect,z0=1-z0);
  z=subst(Ser(precision(1.*z0,digits+3)+X),X,x);
    \\ work with z0+x as a variable gives power series in X as an answer
  Avec=[1,(z+6)/2,(z^2+82*z+96)/6,(z^3+387*z^2+2906*z+1920)/12];
  Bvec=[1,4,8*z+28,14*z^2+204*z+310];
  Qvec=[0,0,0,Avec[4]/Bvec[4]];
  n=4;
  until(err<0.1^(digits+1.5),         \\ Luke's recursions for psi(x)
    c1=(2*n-1)*(3*(n-1)*z+7*n^2-9*n-6);
    c2=-(2*n-3)*(z-n-1)*(3*(n-1)*z-7*n^2+19*n-4);
    c3=(2*n-1)*(n-3)*(z-n)*(z-n-1)*(z+n-4);
    c0=(2*n-3)*(n+1);
    Avec=concat(Avec,[(c1*Avec[n]+c2*Avec[n-1]+c3*Avec[n-2])/c0]);
    Bvec=concat(Bvec,[(c1*Bvec[n]+c2*Bvec[n-1]+c3*Bvec[n-2])/c0]);
    Qvec=concat(Qvec,Avec[n+1]/Bvec[n+1]);
    err=vecmax(abs(Vec(Qvec[n+1]-Qvec[n])));
    n++;
  );
  res=gamma(z0)*exp(intformal( psi(1)+2*(z-1)/z*Qvec[n] )); \\ psi->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,digits,1);
  )))));
  default(seriesprecision,srprec,1);
  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,
    digits,GSD)=
  digits=default(realprecision,,1);
  if (lastFGSs!=ss || lastFGSterms!=extraterms,
    GSD=sum(j=1,numpoles,(abs((ss+poles[j])/2-round(real((ss+poles[j])/2)))<10^(2-digits)) * 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,
     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,,1);
  default(seriesprecision,2*d+2,1);
  deltapol=subst(Vec( (sinh(x)/x)^tvar ),tvar,x);
  default(seriesprecision,OldSeriesPrecision,1);

  \\ 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,
    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,
    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,
    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,,1) +
     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(res>err,t1=t,t2=t),if(res<err,t1=t2/2,t2*=2));
  );
  ceil(t2)
}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ initLdata(cfstr, cutoff=1.2 [,cfdualstr])                             \\
\\                                                                        \\
\\ - to be called before the L-values computations.                       \\
\\ - gammaV, conductor and weight must be initialized by now.             \\
\\   also coefgrow(), MaxImaginaryPart and MaxAsympCoeffs are used here   \\
\\ - CFSTR must be a string which evaluates to a function of k            \\
\\   which gives k-th coefficient of the L-series, e.g. "(-1)^k"          \\
\\ - 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.                                                   \\
\\ - if cutoff<0, force the number of coefficients to be -cutoff          \\
\\ - CFDUALSTR must evaluate (like cfstr) to the k-th coefficient of      \\
\\   the dual L-function if it is different from L(s),                    \\
\\   for instance initLdata("a(k)",,"conj(a(k))")   (see e.g. ex-chgen)   \\
\\ - uses current real precision to determine the desired precision       \\
\\   for L-values                                                         \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
initLdata(vstr,cutoff=1.2,vdualstr="",
    len,d,pordtmp,recF,recG,terms)=

  if (type(gammaV)!="t_VEC",
    error("Gamma-factor gammaV has to be defined before calling initLdata()"));
  if (type(conductor)!="t_INT" && type(conductor)!="t_REAL",
    error("conductor has to be defined before calling initLdata()"));

  len=if(cutoff<0,-cutoff,cflength(cutoff));
  cfvec=vector(len,k,eval(vstr));
  if (vdualstr=="",cfdualvec=cfvec,cfdualvec=vector(len,k,eval(vdualstr)));
  if (cutoff<0,len=cflength());

  lastFGs  = -1E99;         \\ globals to track what was calculated last
  lastFGSs = -1E99;         \\ to avoid re-calculating same values each time
  lastGs   = [-1E99,-1E99];
  lastLSs  = [-1E99,-1E99];

  d       = length(gammaV);                   \\ d = number of gamma factors

  \\ Calculate the necessary amount of extra digits

  answerdigits = default(realprecision,,1);
  vA=sqrt(conductor/Pi^d);
  lfundigits = answerdigits +
     max(ceil(-log(abs(fullgamma(0.7*weight+MaxImaginaryPart*I)))/log(10)),0);
  termdigits   = lfundigits + floor(weight-1);
  taylordigits = 2*termdigits;
  asympdigits  = termdigits;

  \\ Exponential factor defined to maximal precision

  default(realprecision,taylordigits,1);
  vA     = sqrt(precision(conductor,taylordigits)/Pi^d); \\ exp. factor
  lastt  = len/vA;

  pordtmp = vector(d,k,1);                    \\ locate poles and their orders
  for(j=1,d,for(k=1,d,if(j!=k,\
    if(type(diff=gammaV[j]-gammaV[k])=="t_INT" & (diff%2==0) & (diff<=0),\
      pordtmp[j]+=pordtmp[k];pordtmp[k]=0))));
  poles      = [];
  PoleOrders = [];
  for(j=1,d,if(pordtmp[j]!=0,\
    poles=concat(poles,gammaV[j]);PoleOrders=concat(PoleOrders,pordtmp[j])));
  numpoles   = length(poles);

  \\ Initialize the asymptotic coefficients at infinity

  default(realprecision,asympdigits,1);

  recFG=RecursionsAtInfinity(gammaV);
  recF=recFG[1];
  recG=recFG[2];
  kill(recFG);

  \\ Maximum number of terms in the asymptotic expansion

  ncoeff=MaxAsympCoeffs;

  \\ Asymptotic behaviour at infinity for phi and G

  expdifff    = (sum(k=1,d,gammaV[k])+1)/d-1;
  asympconstf = 2*prod(k=1,d,gamma(k/d));
  expdiffg    = (sum(k=1,d,gammaV[k])+1)/d-1-2/d;
  asympconstg = prod(k=1,d,gamma(k/d));

  \\ Coefficients for the asymptotic expansion of phi(t) and G(t)

  Fvec=vector(d+ncoeff,X,0);
  Fvec[d]=1.;
  for(y=1,ncoeff,Fvec[d+y]=1.*sum(j=1,d-1,subst(recF[j],n,y)*Fvec[d+y-j]));

  Gvec=vector(d+ncoeff,X,0);
  Gvec[d]=1.;
  for(y=1,ncoeff,Gvec[d+y]=1.*sum(j=1,d,subst(recG[j],n,y)*Gvec[d+y-j]));

  \\ Convert the Fvec (Taylor asymptotic) coefficients into fcf (contfrac coeffs)

  fcf=SeriesToContFrac(vector(ncoeff+1,k,Fvec[d+k-1]));
  fncf=length(fcf)\2;     \\ at most ncoeff+1, less if terminates

  \\ Taylor series coefficients of phi(t) around t=infinity

  if (lastt<35,termstep=1,termstep=floor(lastt^(1/3)));
  phiinfterms=vector(round(lastt/termstep)+1,k,-1);

  terms=fncf;
  PhiCaseBound=0;
  for (k=1,length(phiinfterms),
    t1=(k-1)*termstep;
    while ((k>1)&&(terms>1)&&
      (abs(phiinf(t1,terms-1)-phiinf(t1,terms)))<10^(-termdigits-1),terms-=1);
    if (sum(j=1,terms,fcf[2*j])<ncoeff,phiinfterms[k]=terms,PhiCaseBound=k*termstep);
  );

  \\ Recursions for phi(t) and G(t,s) at the origin

  default(realprecision,taylordigits,1);

  \\ Initial values of the gamma factors for recursions

  InitV = vector(numpoles,j,prod(k=1,d,
          subst(gammaseries((-poles[j]+gammaV[k])/2,PoleOrders[j]-1),x,X/2)));

  \\ Taylor series coefficients of phi(t) around t=0 -> 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,1);
}



\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ phi(t) - inverse Mellin transform of the product of Gamma-factors      \\
\\          computed either with Taylor in 0 or from asymptotics          \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

phi(t)=if(t<PhiCaseBound,phi0(t),phiinf(t,phiinfterms[min(1+floor(abs(t)/termstep),length(phiinfterms))]));

{
RecursephiV()=     \\ compute one more term for the recursions at the origin
  phiVnn++;
  phiV=concat(phiV,[matrix(numpoles,vecmax(PoleOrders),j,k,polcoeff(phiVser[j],-k))]);
  for (j=1,numpoles,for(k=1,length(gammaV),
    phiVser[j]/=(X/2-poles[j]/2-phiVnn+gammaV[k]/2)));
}

{
phi0(t,                              \\ phi(t) using series expansion at t=0
    t2,LogTTerm,TPower,res=0,nn=0,totalold)=
  default(realprecision,taylordigits,1);
  t        = precision(t,taylordigits);
  t2       = t^2;
  LogTTerm = vector(vecmax(PoleOrders),k,(-log(t))^(k-1)/(k-1)!)~;
  TPower   = 1.0*vector(numpoles,j,t^poles[j]);
  until (abs(res-totalold)<10^-(termdigits+1)&&(nn>3),
    totalold=res;
    nn++;
    if(nn>phiVnn,RecursephiV());
    res+=TPower*phiV[nn]*LogTTerm;
    TPower*=t2;
  );
  default(realprecision,termdigits,1);
  res
}

{                          \\ phi(t) using asymptotic expansion at infinity
phiinf(t,ncf=fncf,
    res,d,td2) =
  default(realprecision,asympdigits,1);
  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,1);
  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,
    nn)=
  if(lastGs!=[ss,der] || type(Ginfterms)!="t_VEC",
    initGinf(ss,der);lastGs=[ss,der]);
  if(t<GCaseBound,G0(t,ss,der),
    nn=min(1+floor(abs(t)/termstep),length(Ginfterms));
    Ginf(t,ss,der,Ginfterms[nn]));
}

LogInt(i,j,logt,der)=\
  if (abs(i)<10^(2-termdigits),0,sum(k=0,j-1,binomial((k-j),der)*der!*(-i*logt)^k/k!)/i^(j+der));

{
MakeLogSum(ss,der,
    nn,V,logsum)=
  if(length(LogSum)==phiVnn,                      \\ more phiV's necessary
    for (j=1,floor(phiVnn/10)+1,RecursephiV()));
  for (nn=length(LogSum)+1,phiVnn,                \\ generate logsums
    V=phiV[nn];
    logsum=vector(numpoles,j,sum(k=1,PoleOrders[j],V[j,k]*LogInt(poles[j]+2*(nn-1)+ss,k,lt,der)))~;
    LogSum=concat(LogSum,[logsum]);
  );
  lastLSs=[ss,der];
}

{
G0(t,ss,der,                 \\ G(t,s,der) computed using Taylor series at 0
    t2,LT,TPower,res,nn,term,gmser,gmcf,dgts)=
  default(realprecision,taylordigits,1);
  ss     = precision(ss,taylordigits);
  if ([ss,der]!=lastLSs,LogSum=[]);
  t      = precision(t,taylordigits);
  t2     = t^2;       \\ time
  LT     = log(t);    \\ = money
  TPower = vector(numpoles,j,t^poles[j]);
  res    = 0;
  nn     = 0;
  term   = 1;
  until ((nn>3) && abs(term)<10^-(termdigits+1),
    nn++;
    if(nn>length(LogSum),MakeLogSum(ss,der));
    term=TPower*subst(LogSum[nn],lt,LT);
    res+=term;
    TPower*=t2;
  );
  gmser=fullgammaseries(ss,der)/t^(S+ss);
  gmcf=polcoeff(gmser,der,S)*der!;
  res=(gmcf-res);
  default(realprecision,termdigits,1);
  res
}


{
Ginf(t,ss,der,ncf=-1,     \\ G(t,s,der) computed using asymptotic expansion
    res,d,tt) =           \\ at infinity and associated continued fraction
  default(realprecision,asympdigits,1);
  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,1);
  res;
}


{
initGinf(ss,der,          \\ pre-compute asymptotic expansions for a given s
    d,gvec,gncf,terms,t1)=
  default(realprecision,asympdigits,1);
  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])<ncoeff,Ginfterms[k]=terms,GCaseBound=k*termstep);
  );
  default(realprecision,termdigits,1);
}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ ltheta(t) = sum a(k)*phi(k*t/vA)                                        \\
\\ satisfies ltheta(1/t)=sgn*t^weight*ldualtheta(t) + residue contribution \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

ltheta(t) = t=precision(t/vA,taylordigits);\
  sum(k=1,min(floor(lastt/t+1),length(cfvec)),cfvec[k]*if(cfvec[k],phi(k*t),0));
ldualtheta(t) = t=precision(t/vA,taylordigits);\
  sum(k=1,min(floor(lastt/t+1),length(cfvec)),cfdualvec[k]*if(cfdualvec[k],phi(k*t),0));

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ checkfeq(t=1.2) = verify the functional equation by evaluating LHS-RHS  \\
\\                   for func.eq for ltheta(t), should return approx. 0    \\
\\ - also determines residues if Lresidues is set to automatic             \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
checkfeq(t=6/5,             \\ determine residues if they are not yet set
    nlp,lpx,lpv,lpm,res)=   \\ .. and check the functional equation
  if (Lresidues==automatic && Lpoles!=[],
    nlp=length(Lpoles);
    lpx=vector(nlp,k,1.15+if(k==1,0,k-1)/10);
    lpv=vector(nlp,k,tq=lpx[k];sgn*tq^weight*ldualtheta(tq)-ltheta(1/tq));
    lpm=matrix(nlp,nlp,k,j,tq=lpx[k];ss=Lpoles[j];tq^ss-sgn*tq^(weight-ss));
    Lresidues=matsolve(lpm,lpv~)~;
    \\for (k=1,nlp,print("Residue at ",Lpoles[k]," = ",Lresidues[k]));
  );
  res=ltheta(t)-sgn*t^(-weight)*ldualtheta(1/t)+
    sum(k=1,length(Lpoles),Lresidues[k]*(t^-Lpoles[k]-sgn*t^(-weight+Lpoles[k])));
  default(realprecision,answerdigits,1);
  res;
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ L(ss,cutoff,k) = k-th derivative of L(s) at s=ss                       \\
\\                  cutoff = 1 by default (cutoff point), >=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,
    FGSeries,LGSeries,res)=
  default(realprecision,lfundigits,1);
  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,1);
  res;
}


{
Lstar(ss,cutoff=1,der=0,                \\ Lstar(s) = L(s) * Gamma factor
    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,1);
  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;
}