/*-*- compile-command: "/usr/bin/gcc -c -o computel5.gp.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I"/usr/local/include" computel5.gp.c && /usr/bin/gcc -o computel5.gp.so -shared -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -Wl,-shared computel5.gp.o -lc -lm -L"/usr/local/lib" -lpari"; -*-*/
#include <pari/pari.h>
/*
GP;install("init_computel5","v","init_computel5","./computel5.gp.so");
GP;install("initglobs","","initglobs","./computel5.gp.so");
GP;install("getconductor","","getconductor","./computel5.gp.so");
GP;install("getgammaV","","getgammaV","./computel5.gp.so");
GP;install("getweight","","getweight","./computel5.gp.so");
GP;install("getsgn","","getsgn","./computel5.gp.so");
GP;install("getLpoles","","getLpoles","./computel5.gp.so");
GP;install("getLresidues","","getLresidues","./computel5.gp.so");
GP;install("getMaxImaginaryPart","","getMaxImaginaryPart","./computel5.gp.so");
GP;install("getMaxAsympCoeffs","","getMaxAsympCoeffs","./computel5.gp.so");
GP;install("setglobs","D0,G,","setglobs","./computel5.gp.so");
GP;install("initall","D0,G,","initall","./computel5.gp.so");
GP;install("setconductor","D0,G,","setconductor","./computel5.gp.so");
GP;install("setgammaV","D0,G,","setgammaV","./computel5.gp.so");
GP;install("setweight","D0,G,","setweight","./computel5.gp.so");
GP;install("setsgn","D0,G,","setsgn","./computel5.gp.so");
GP;install("setLpoles","D0,G,","setLpoles","./computel5.gp.so");
GP;install("setLresidues","D0,G,","setLresidues","./computel5.gp.so");
GP;install("setMaxImaginaryPart","D0,G,","setMaxImaginaryPart","./computel5.gp.so");
GP;install("setMaxAsympCoeffs","D0,G,","setMaxAsympCoeffs","./computel5.gp.so");
GP;install("setmycoefgrow","D0,G,","setmycoefgrow","./computel5.gp.so");
GP;install("coefgrow","D0,G,p","coefgrow","./computel5.gp.so");
GP;install("vecleft","D0,G,D0,G,","vecleft","./computel5.gp.so");
GP;install("vecright","D0,G,D0,G,","vecright","./computel5.gp.so");
GP;install("StrTab","D0,G,D0,G,","StrTab","./computel5.gp.so");
GP;install("concatstr","DGDGDGDG","concatstr","./computel5.gp.so");
GP;install("errprint","D0,G,p","errprint","./computel5.gp.so");
GP;install("gammaseries","D0,G,D0,G,p","gammaseries","./computel5.gp.so");
GP;install("fullgamma","D0,G,p","fullgamma","./computel5.gp.so");
GP;install("fullgammaseries","D0,G,D0,G,p","fullgammaseries","./computel5.gp.so");
GP;install("RecursionsAtInfinity","p","RecursionsAtInfinity","./computel5.gp.so");
GP;install("SeriesToContFrac","D0,G,p","SeriesToContFrac","./computel5.gp.so");
GP;install("EvaluateContFrac","D0,G,D0,G,D0,G,p","EvaluateContFrac","./computel5.gp.so");
GP;install("cflength","DGp","cflength","./computel5.gp.so");
GP;install("initLdata","lD0,G,DGDGp","initLdata","./computel5.gp.so");
GP;install("phi","D0,G,p","phi","./computel5.gp.so");
GP;install("RecursephiV","v","RecursephiV","./computel5.gp.so");
GP;install("phi0","D0,G,p","phi0","./computel5.gp.so");
GP;install("phiinf","D0,G,DGp","phiinf","./computel5.gp.so");
GP;install("G","D0,G,D0,G,DGp","G","./computel5.gp.so");
GP;install("LogInt","D0,G,D0,G,D0,G,D0,G,p","LogInt","./computel5.gp.so");
GP;install("MakeLogSum","D0,G,D0,G,p","MakeLogSum","./computel5.gp.so");
GP;install("G0","D0,G,D0,G,D0,G,p","G0","./computel5.gp.so");
GP;install("Ginf","D0,G,D0,G,D0,G,DGp","Ginf","./computel5.gp.so");
GP;install("initGinf","lD0,G,D0,G,p","initGinf","./computel5.gp.so");
GP;install("ltheta","D0,G,p","ltheta","./computel5.gp.so");
GP;install("ldualtheta","D0,G,p","ldualtheta","./computel5.gp.so");
GP;install("checkfeq","DGp","checkfeq","./computel5.gp.so");
GP;install("L","D0,G,DGDGp","L","./computel5.gp.so");
GP;install("Lseries","D0,G,DGDGp","Lseries","./computel5.gp.so");
GP;install("Lstar","D0,G,DGDGp","Lstar","./computel5.gp.so");
*/
void init_computel5(void);
GEN initglobs(void);
GEN getconductor(void), getgammaV(void), getweight(void), getsgn(void), getLpoles(void), getLresidues(void), getMaxImaginaryPart(void), getMaxAsympCoeffs(void);
GEN setglobs(GEN v);
GEN initall(GEN v);
GEN setconductor(GEN c), setgammaV(GEN c), setweight(GEN c), setsgn(GEN c), setLpoles(GEN c), setLresidues(GEN c), setMaxImaginaryPart(GEN c), setMaxAsympCoeffs(GEN c), setmycoefgrow(GEN c);
GEN coefgrow(GEN n, long prec);
GEN vecleft(GEN v, GEN n);
GEN vecright(GEN v, GEN n);
GEN StrTab(GEN x, GEN n);
GEN concatstr(GEN s1, GEN s2, GEN s3, GEN s4);
GEN errprint(GEN x, long prec);
GEN gammaseries(GEN z0, GEN terms, long prec);
GEN fullgamma(GEN ss, long prec);
GEN fullgammaseries(GEN ss, GEN extraterms, long prec);
GEN RecursionsAtInfinity(long prec);
GEN SeriesToContFrac(GEN vec, long prec);
GEN EvaluateContFrac(GEN cfvec, GEN terms, GEN t, long prec);
GEN cflength(GEN cutoff, long prec);
long initLdata(GEN vstr, GEN cutoff, GEN vdualstr, long prec);
GEN phi(GEN t, long prec);
void RecursephiV(void);
GEN phi0(GEN t, long prec);
GEN phiinf(GEN t, GEN ncf, long prec);
GEN G(GEN t, GEN ss, GEN der, long prec);
GEN LogInt(GEN i, GEN j, GEN logt, GEN der, long prec);
GEN MakeLogSum(GEN ss, GEN der, long prec);
GEN G0(GEN t, GEN ss, GEN der, long prec);
GEN Ginf(GEN t, GEN ss, GEN der, GEN ncf, long prec);
long initGinf(GEN ss, GEN der, long prec);
GEN ltheta(GEN t, long prec);
GEN ldualtheta(GEN t, long prec);
GEN checkfeq(GEN t, long prec);
GEN L(GEN ss, GEN cutoff, GEN der, long prec);
GEN Lseries(GEN ss, GEN cutoff, GEN der, long prec);
GEN Lstar(GEN ss, GEN cutoff, GEN der, long prec);
/*End of prototype*/

static GEN conductor;
static GEN gammaV;
static GEN weight;
static GEN sgn;
static GEN Lpoles;
static GEN Lresidues;
static GEN MaxImaginaryPart;
static GEN MaxAsympCoeffs;
static GEN lastFGs;
static GEN lastFGSs;
static GEN lastGs;
static GEN lastLSs;
static GEN lastFGval;
static GEN lastFGSval;
static GEN lastFGSterms;
static GEN cfvec;
static GEN cfdualvec;
static GEN answerdigits;
static GEN vA;
static GEN lfundigits;
static GEN termdigits;
static GEN taylordigits;
static GEN asympdigits;
static GEN lastt;
static GEN poles;
static GEN PoleOrders;
static GEN numpoles;
static GEN ncoeff;
static GEN expdifff;
static GEN asympconstf;
static GEN expdiffg;
static GEN asympconstg;
static GEN Fvec;
static GEN Gvec;
static GEN fcf;
static GEN fncf;
static GEN termstep;
static GEN phiinfterms;
static GEN PhiCaseBound;
static GEN InitV;
static GEN phiV;
static GEN phiVnn;
static GEN phiVser;
static GEN Ginfterms;
static GEN GCaseBound;
static GEN LogSum;
static GEN gncf;
static GEN gcf_T;
static GEN mycoefgrow;
/*End of global vars*/

void
init_computel5(void)	  /* void */
{
  /*********** 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! -> tim.dokchitser@durham.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 */
  
  MaxAsympCoeffs = pol_x(fetch_user_var("MaxAsympCoeffs"));
  MaxImaginaryPart = pol_x(fetch_user_var("MaxImaginaryPart"));
  Lresidues = pol_x(fetch_user_var("Lresidues"));
  Lpoles = pol_x(fetch_user_var("Lpoles"));
  sgn = pol_x(fetch_user_var("sgn"));
  weight = pol_x(fetch_user_var("weight"));
  gammaV = pol_x(fetch_user_var("gammaV"));
  conductor = pol_x(fetch_user_var("conductor"));
  lastLSs = pol_x(fetch_user_var("lastLSs"));
  lastGs = pol_x(fetch_user_var("lastGs"));
  lastFGSs = pol_x(fetch_user_var("lastFGSs"));
  lastFGs = pol_x(fetch_user_var("lastFGs"));
  lastFGSterms = pol_x(fetch_user_var("lastFGSterms"));
  lastFGSval = pol_x(fetch_user_var("lastFGSval"));
  lastFGval = pol_x(fetch_user_var("lastFGval"));
  phiVser = pol_x(fetch_user_var("phiVser"));
  phiVnn = pol_x(fetch_user_var("phiVnn"));
  phiV = pol_x(fetch_user_var("phiV"));
  InitV = pol_x(fetch_user_var("InitV"));
  PhiCaseBound = pol_x(fetch_user_var("PhiCaseBound"));
  phiinfterms = pol_x(fetch_user_var("phiinfterms"));
  termstep = pol_x(fetch_user_var("termstep"));
  fncf = pol_x(fetch_user_var("fncf"));
  fcf = pol_x(fetch_user_var("fcf"));
  Gvec = pol_x(fetch_user_var("Gvec"));
  Fvec = pol_x(fetch_user_var("Fvec"));
  asympconstg = pol_x(fetch_user_var("asympconstg"));
  expdiffg = pol_x(fetch_user_var("expdiffg"));
  asympconstf = pol_x(fetch_user_var("asympconstf"));
  expdifff = pol_x(fetch_user_var("expdifff"));
  ncoeff = pol_x(fetch_user_var("ncoeff"));
  numpoles = pol_x(fetch_user_var("numpoles"));
  PoleOrders = pol_x(fetch_user_var("PoleOrders"));
  poles = pol_x(fetch_user_var("poles"));
  lastt = pol_x(fetch_user_var("lastt"));
  asympdigits = pol_x(fetch_user_var("asympdigits"));
  taylordigits = pol_x(fetch_user_var("taylordigits"));
  termdigits = pol_x(fetch_user_var("termdigits"));
  lfundigits = pol_x(fetch_user_var("lfundigits"));
  vA = pol_x(fetch_user_var("vA"));
  answerdigits = pol_x(fetch_user_var("answerdigits"));
  cfdualvec = pol_x(fetch_user_var("cfdualvec"));
  cfvec = pol_x(fetch_user_var("cfvec"));
  GCaseBound = pol_x(fetch_user_var("GCaseBound"));
  Ginfterms = pol_x(fetch_user_var("Ginfterms"));
  LogSum = pol_x(fetch_user_var("LogSum"));
  gcf_T = pol_x(fetch_user_var("gcf_T"));
  gncf = pol_x(fetch_user_var("gncf"));
  mycoefgrow = pol_x(fetch_user_var("mycoefgrow"));
  return;
}

GEN
initglobs(void)
{
  GEN automatic = pol_x(fetch_user_var("automatic"));
  conductor = gclone(automatic);
  gammaV = gclone(automatic);
  weight = gen_1;
  sgn = gen_1;
  Lpoles = gclone(cgetg(1, t_VEC));
  Lresidues = gclone(automatic);
  MaxImaginaryPart = gen_0;
  MaxAsympCoeffs = gclone(stoi(40));
  mycoefgrow = gen_0;
  return gen_0;
}

GEN
setglobs(GEN v)
{
  conductor = gclone(gel(v, 1));
  gammaV = gclone(gel(v, 2));
  if (glength(v) >= 3) weight = gclone(gel(v, 3));
  if (glength(v) >= 4) sgn = gclone(gel(v, 4));
  if (glength(v) >= 5) Lpoles = gclone(gel(v, 5));
  if (glength(v) >= 6) Lresidues = gclone(gel(v, 6));
  if (glength(v) >= 7) MaxImaginaryPart = gclone(gel(v, 7));
  if (glength(v) >= 8) MaxAsympCoeffs = gclone(gel(v, 8));
  if (glength(v) >= 9) mycoefgrow = gclone(gel(v, 9));
  return gen_0;
}

GEN
initall(GEN v)
{
  init_computel5();
  initglobs();
  setglobs(v);
  return gen_0;
}

GEN
setconductor(GEN c)
{
  return (conductor = gclone(c));
}

GEN
getconductor(void)
{
  return conductor;
}

GEN
setgammaV(GEN c)
{
    return (gammaV = gclone(c));
}

GEN
getgammaV(void)
{
  return gammaV;
}

GEN
setweight(GEN c)
{
  return (weight = gclone(c));
}

GEN
getweight(void)
{
  return weight;
}

GEN
setsgn(GEN c)
{
  return (sgn = gclone(c));
}

GEN
getsgn(void)
{
  return sgn;
}

GEN
setLpoles(GEN c)
{
  return (Lpoles = gclone(c));
}

GEN
getLpoles(void)
{
  return Lpoles;
}

GEN
setLresidues(GEN c)
{
  return (Lresidues = gclone(c));
}

GEN
getLresidues(void)
{
  return Lresidues;
}

GEN
setMaxImaginaryPart(GEN c)
{
  return (MaxImaginaryPart = gclone(c));
}

GEN
getMaxImaginaryPart(void)
{
  return MaxImaginaryPart;
}

GEN
setMaxAsympCoeffs(GEN c)
{
  return(MaxAsympCoeffs = gclone(c));
}

GEN
getMaxAsympCoeffs(void)
{
  return MaxAsympCoeffs;
}

GEN
setmycoefgrow(GEN c)
{
  return(mycoefgrow = gclone(c));
}

/* by default L(s)<->L(1-s) */
/*            with sign=1 in functional equation */
/*            and L*(s) has no poles */
/* if this is not changed to [r1,r2,...] by hand, */
/* checkfeq() tries to determine residues automatically */
/* see ex-nf for instance */

GEN
coefgrow(GEN n, long prec)
{
  GEN p1;
  if (typ(mycoefgrow)!=t_INT) return closure_callgen1(mycoefgrow,n);
  if (glength(Lpoles))
    p1 = gmul(/* default bound for coeffs. a(n) */
    strtor("1.5", prec), gpow(n, gsubgs(vecmax(greal(Lpoles)), 1), prec));
  else
    p1 = gpow(/* you may redefine coefgrow() by hand */
    gsqrt(gmulsg(4, n), prec), gsubgs(weight, 1), prec);
  return p1;
}

/* 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) */

/* re-define this if you want to compute L(s) */
/* for large imaginary s (see ex-zeta2 for example) */

/* 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 */
GEN
vecleft(GEN v, GEN n)	  /* vec */
{
  return extract0(v, concat(strtoGENstr("1.."), GENtoGENstr(n)), NULL);
}

GEN
vecright(GEN v, GEN n)	  /* vec */
{
  return extract0(v, concat(GENtoGENstr(gaddgs(gsubsg(glength(v), n), 1)), concat(strtoGENstr(".."), GENtoGENstr(stoi(glength(v))))), NULL);
}

/* Tabulate a string to n characters, e.g. StrTab(3,2)="3 "; */
GEN
StrTab(GEN x, GEN n)
{
  x = GENtoGENstr(x);
  while (gcmpsg(glength(x), n) < 0)
    x = concat(x, strtoGENstr(" "));
  return x;
}

/* Concatenate up to 4 strings */
GEN
concatstr(GEN s1, GEN s2, GEN s3, GEN s4)	  /* genstr */
{
  if (!s1) s1 = strtoGENstr("");
  if (!s2) s2 = strtoGENstr("");
  if (!s3) s3 = strtoGENstr("");
  if (!s4) s4 = strtoGENstr("");
  if (typ(s1)!=t_STR) s1=GENtoGENstr(s1);
  if (typ(s2)!=t_STR) s2=GENtoGENstr(s2);
  if (typ(s3)!=t_STR) s3=GENtoGENstr(s3);
  if (typ(s4)!=t_STR) s4=GENtoGENstr(s4);
  return concat(s1, concat(s2, concat(s3, s4)));
}

/* Print a ``small error'', e.g. 0.00000013 as "1E-7" */
GEN
errprint(GEN x, long prec)
{
  GEN p1;
  if (typ(x) == t_COMPLEX) x = gabs(x, prec);
  if (gequal0(x))
    p1 = concatstr(strtoGENstr("1E-"), stoi(getrealprecision() + 1), NULL, NULL);
  else
    p1 = concatstr(gtrunc(gdiv(x, gpow(stoi(10), gfloor(gdiv(glog(gabs(x, prec), prec), glog(stoi(10), prec))), prec))), strtoGENstr("E"), gfloor(gdiv(glog(gabs(x, prec), prec), glog(stoi(10), prec))), NULL);
  return p1;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* 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] \ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
gammaseries(GEN z0, GEN terms, long prec)
{
  pari_sp ltop = avma;
  GEN Avec = gen_0, Bvec = gen_0, Qvec = gen_0, n = gen_0, z = gen_0, err = gen_0, res = gen_0, c0 = gen_0, c1 = gen_0, c2 = gen_0, c3 = gen_0, sinser = gen_0, reflect = gen_0, digits = gen_0, srprec = gen_0, negint = gen_0, x = pol_x(fetch_user_var("x")), X = pol_x(fetch_user_var("X"));
  GEN p1 = gen_0, p2 = gen_0, p3 = gen_0;	  /* vec */
  srprec = stoi(precdl);
  if (gequal(z0, greal(ground(z0)))) z0 = greal(ground(z0));
  /* you don't want to know */
  negint = stoi((typ(z0) == t_INT) && (gcmpgs(z0, 0) <= 0));
  /* z0 is a pole */
  precdl = gtos(gadd(gaddgs(terms, 1), negint));
  if (gequal0(terms) && (gequal0(negint)))
    res = gadd(ggamma(z0, prec), ggrando(x, 1));
  else
  {
    if ((gequal1(terms) && (gequal0(gimag(z0)))) && (gequal0(negint)))
      /*   built-in functions */
      res = gmul(ggamma(z0, prec), gadd(gaddsg(1, gmul(gpsi(z0, prec), x)), ggrando(x, 2)));
    else
    {
      if (gequal0(z0))
        res = gdiv(ggamma(gaddsg(1, x), prec), x);
      else
      {
        if (gequal1(z0))
          res = ggamma(gaddsg(1, x), prec);
        else
        {
          if (gequalgs(z0, 2))
            res = gmul(ggamma(gaddsg(1, x), prec), gaddsg(1, x));
          else
          {
            /* otherwise use Luke's rational approximations for psi(x) */
            digits = stoi(getrealprecision());
            /* save working precision */
            setrealprecision(gtos(gaddgs(digits, 3)), &prec);
            /*   and work with 3 digits more */
            reflect = stoi(gcmp(greal(z0), strtor("0.5", prec)) < 0);
            /* left of 1/2 use reflection formula */
            if (!gequal0(reflect))
              z0 = gsubsg(1, z0);
            z = gsubst(gtoser(gadd(precision0(gmul(real_1(prec), z0), gtos(gaddgs(digits, 3))), X), -1, precdl), gvar(X), x);
            p1 = cgetg(5, t_VEC);
            gel(p1, 1) = gen_1;
            gel(p1, 2) = gdivgs(gaddgs(z, 6), 2);
            gel(p1, 3) = gdivgs(gaddgs(gadd(gsqr(z), gmulsg(82, z)), 96), 6);
            gel(p1, 4) = gdivgs(gaddgs(gadd(gadd(gpowgs(z, 3), gmulsg(387, gsqr(z))), gmulsg(2906, z)), 1920), 12);
            /* work with z0+x as a variable gives power series in X as an answer */
            Avec = p1;
            p2 = cgetg(5, t_VEC);
            gel(p2, 1) = gen_1;
            gel(p2, 2) = stoi(4);
            gel(p2, 3) = gaddgs(gmulsg(8, z), 28);
            gel(p2, 4) = gaddgs(gadd(gmulsg(14, gsqr(z)), gmulsg(204, z)), 310);
            Bvec = p2;
            p3 = cgetg(5, t_VEC);
            gel(p3, 1) = gen_0;
            gel(p3, 2) = gen_0;
            gel(p3, 3) = gen_0;
            gel(p3, 4) = gdiv(gel(Avec, 4), gel(Bvec, 4));
            Qvec = p3;
            n = stoi(4);
            {
              pari_sp btop = avma, st_lim = stack_lim(btop, 1);
              GEN p4 = gen_0, p5 = gen_0;	  /* vec */
              do
              {
                /* Luke's recursions for psi(x) */
                c1 = gmul(gsubgs(gmulsg(2, n), 1), gsubgs(gsub(gadd(gmul(gmulsg(3, gsubgs(n, 1)), z), gmulsg(7, gsqr(n))), gmulsg(9, n)), 6));
                c2 = gmul(gmul(gneg(gsubgs(gmulsg(2, n), 3)), gsubgs(gsub(z, n), 1)), gsubgs(gadd(gsub(gmul(gmulsg(3, gsubgs(n, 1)), z), gmulsg(7, gsqr(n))), gmulsg(19, n)), 4));
                c3 = gmul(gmul(gmul(gmul(gsubgs(gmulsg(2, n), 1), gsubgs(n, 3)), gsub(z, n)), gsubgs(gsub(z, n), 1)), gsubgs(gadd(z, n), 4));
                c0 = gmul(gsubgs(gmulsg(2, n), 3), gaddgs(n, 1));
                p4 = cgetg(2, t_VEC);
                gel(p4, 1) = gdiv(gadd(gadd(gmul(c1, gel(Avec, gtos(n))), gmul(c2, gel(Avec, gtos(gsubgs(n, 1))))), gmul(c3, gel(Avec, gtos(gsubgs(n, 2))))), c0);
                Avec = concat(Avec, p4);
                p5 = cgetg(2, t_VEC);
                gel(p5, 1) = gdiv(gadd(gadd(gmul(c1, gel(Bvec, gtos(n))), gmul(c2, gel(Bvec, gtos(gsubgs(n, 1))))), gmul(c3, gel(Bvec, gtos(gsubgs(n, 2))))), c0);
                Bvec = concat(Bvec, p5);
                Qvec = concat(Qvec, gdiv(gel(Avec, gtos(gaddgs(n, 1))), gel(Bvec, gtos(gaddgs(n, 1)))));
                err = vecmax(gabs(gtovec(gsub(gel(Qvec, gtos(gaddgs(n, 1))), gel(Qvec, gtos(n)))), prec));
                n = gaddgs(n, 1);
                if (low_stack(st_lim, stack_lim(btop, 1)))
                  gerepileall(btop, 11, &c1, &c2, &c3, &c0, &p4, &Avec, &p5, &Bvec, &Qvec, &err, &n);
              } while(!(gcmp(err, gpow(strtor("0.1", prec), gadd(digits, strtor("1.5", prec)), prec)) < 0));
            }
            res = gmul(ggamma(z0, prec), gexp(integ(gadd(gpsi(gen_1, prec), gmul(gdiv(gmulsg(2, gsubgs(z, 1)), z), gel(Qvec, gtos(n)))), -1), prec));
            /* psi->gamma */
            if (!gequal0(reflect))
            {
              /* reflect if necessary */
              sinser = gtovec(gsin(gmul(mppi(prec), z), prec));
              if (!gequal0(negint))
                gel(sinser, 1) = gen_0;
              /* taking slight care at integers<0 */
              res = gsubst(gdiv(gdiv(mppi(prec), res), gtoser(sinser, -1, precdl)), gvar(x), gneg(x));
            }
            setrealprecision(gtos(digits), &prec);
          }
        }
      }
    }
  }
  precdl = gtos(srprec);
  res = gerepilecopy(ltop, res);
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* fullgamma(ss) - the full gamma factor (at s=ss)                        \\ */
/*   vA^s*Gamma((s+gammaV[1])/2)*...*Gamma((s+gammaV[d])/2)               \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
fullgamma(GEN ss, long prec)
{
  long l1;
  GEN p2 = gen_0;
  if (!gequal(ss, lastFGs))
  {
    lastFGs = gclone(ss);
    l1 = glength(gammaV);
    p2 = gpow(vA, ss, prec);
    {
      pari_sp btop = avma, st_lim = stack_lim(btop, 1);
      long j;
      for (j = 1; j <= l1; ++j)
      {
        p2 = gmul(p2, ggamma(gdivgs(gadd(ss, gel(gammaV, j)), 2), prec));
        if (low_stack(st_lim, stack_lim(btop, 1)))
          p2 = gerepilecopy(btop, p2);
      }
    }
    lastFGval = gclone(p2);
  }
  return 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. */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
fullgammaseries(GEN ss, GEN extraterms, long prec)
{
  GEN digits = gen_0, GSD = gen_0, p1 = gen_0;
  long l2;
  GEN p3 = gen_0, x = pol_x(fetch_user_var("x")), S = pol_x(fetch_user_var("S"));
  digits = stoi(getrealprecision());
  if (!gequal(lastFGSs, ss) || !gequal(lastFGSterms, extraterms))
  {
    {
      long j;
      p1 = gen_0;
      for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
        p1 = gadd(p1, gmulsg(gcmp(gabs(gsub(gdivgs(gadd(ss, gel(poles, j)), 2), ground(greal(gdivgs(gadd(ss, gel(poles, j)), 2)))), prec), gpow(stoi(10), gsubsg(2, digits), prec)) < 0, gel(PoleOrders, j)));
    }
    GSD = gadd(p1, extraterms);
    lastFGSs = gclone(ss);
    lastFGSterms = gclone(extraterms);
    l2 = glength(gammaV);
    {
      pari_sp btop = avma, st_lim = stack_lim(btop, 1);
      long j;
      p3 = gen_1;
      for (j = 1; j <= l2; ++j)
      {
        p3 = gmul(p3, gammaseries(gdivgs(gadd(ss, gel(gammaV, j)), 2), GSD, prec));
        if (low_stack(st_lim, stack_lim(btop, 1)))
          p3 = gerepilecopy(btop, p3);
      }
    }
    lastFGSval = gclone(gsubst(p3, gvar(x), gdivgs(S, 2)));
  }
  return lastFGSval;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/*         RecursionsAtInfinity(gammaV)                                   \\ */
/*   Recursions for the asymptotic expansion coefficients                 \\ */
/*   for phi(x) and G(s,x) at infinity.                                   \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
RecursionsAtInfinity(long prec)	  /* vec */
{
  GEN symvec = gen_0, modsymvec = gen_0, deltapol = gen_0, recF = gen_0, recG = gen_0, OldSeriesPrecision = gen_0, p1 = gen_0, x = pol_x(fetch_user_var("x"));
  GEN p2 = gen_0;	  /* vec */
  GEN p4 = gen_0;	  /* vec */
  GEN tvar = pol_x(fetch_user_var("tvar"));
  GEN p6 = gen_0;	  /* vec */
  GEN n = pol_x(fetch_user_var("n"));
  GEN p7 = gen_0;	  /* vec */
  GEN s = pol_x(fetch_user_var("s"));
  GEN p8 = gen_0;	  /* vec */
  long p9;
  GEN p10 = gen_0;	  /* vec */
  GEN gd;
  long d, j, k, p3, p5;
  d = glength(gammaV); gd = stoi(d);
  p1 = gen_1;
  for (k = 1; k <= d; ++k) p1 = gmul(p1, gadd(x, gel(gammaV, k)));
  p2 = cgetg(2, t_VEC); gel(p2, 1) = gen_0;
  symvec = concat(gtovec(p1), p2);
  p3 = d + 2;
  p4 = cgetg(p3+1, t_VEC);
  for (k = 1; k <= p3; ++k) gel(p4, k) = gen_0;
  modsymvec = p4;
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    for (j = 0; j <= d; ++j)
    {
      long k;
      long l11, l12;
      for (k = 0; k <= j; ++k)
      {
        l11 = j + 1;
	l12 = j - k + 1;
	gel(modsymvec, l11) = gadd(gel(modsymvec, l11), gmul(gmul(gmul(gpowgs(gneg(gel(symvec, 2)), k), gpowgs(gd, j - 1 - k)), binomial(stoi(d + k - j), k)), gel(symvec, l12)));
      }
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 1, &modsymvec);
    }
  }
  /* Delta polynomials */
  OldSeriesPrecision = stoi(precdl);
  precdl = 2*d + 2;
  deltapol = gsubst(gtovec(gpow(gdiv(gsh(x, prec), x), tvar, prec)), gvar(tvar), x);
  precdl = gtos(OldSeriesPrecision);
  p5 = d + 1;
  {
    long p;
    GEN p13 = gen_0;
    p6 = cgetg(p5+1, t_VEC);
    for (p = 1; p <= p5; ++p)
    {
      {
        pari_sp btop = avma, st_lim = stack_lim(btop, 1);
        long m, l14, l15;
        GEN p16 = gen_0, p18 = gen_0;
	long p17;
        p13 = gen_0;
        for (m = 0; m <= p; ++m)
        {
          l14 = m + 1;
          l15 = p - 1;
	  p16 = gen_1;
	  for (j = m; j <= l15; ++j) p16 = gmulgs(p16, d - j);
          p17 = (p - m) >> 1;
          {
            long l19;
            p18 = gen_0;
            for (k = 0; k <= p17; ++k)
            {
              l19 = 2*k + 1;
              p18 = gadd(p18, gmul(gdiv(gpowgs(gaddgs(gsubgs(gmulsg(2, n), p), 1), p - m - 2*k), mpfact(p - m - 2*k)), gsubst(gel(deltapol, l19), gvar(x), stoi(d - p))));
            }
          }
          p13 = gadd(p13, gmul(gmul(gel(modsymvec, l14), p16), p18));
          if (low_stack(st_lim, stack_lim(btop, 1)))
            gerepileall(btop, 3, &p13, &p16, &p18);
        }
      }
      gel(p6, p) = gmul(gdiv(gdiv(gdivsg(-1, gpowgs(gen_2, p)), gpowgs(gd, p - 1)), n), p13);
    }
  }
  /* recursion coefficients for phi at infinity */
  recF = p6;
  {
    long p;
    p7 = cgetg(d+1, t_VEC);
    for (p = 1; p <= d; ++p)
      gel(p7, p) = gsub(gel(recF, p + 1), gmul(gdivgs(gdivgs(gsubgs(gsub(gadd(gel(symvec, 2), gmulsg(d, gsubgs(s, 1))), gmulsg(2, gsubgs(n, p))), 1), 2), d), gel(recF, p)));
  }
  /* recursion coefficients for G at infinity */
  recG = p7;
  p8 = cgetg(3, t_VEC);
  p9 = d - 1;
  {
    long p;
    p10 = cgetg(p9+1, t_VEC);
    for (p = 1; p <= p9; ++p) gel(p10, p) = gcopy(gel(recF, p + 1));
  }
  gel(p8, 1) = p10;
  gel(p8, 2) = gcopy(recG);
  return p8;
}

/* return them both */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/*                 SeriesToContFrac(vec)                                  \\ */
/*    Convert a power series vec[1]+vec[2]*x+vec[3]*x^2+...               \\ */
/*    into a continued fraction expansion.                                \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
SeriesToContFrac(GEN vec, long prec)
{
  pari_sp ltop = avma;
  GEN res = gen_0, ind = gen_0, x = pol_x(fetch_user_var("x"));
  vec = gcopy(vec);
  vec = gmul(real_1(prec), vec);
  res = cgetg(1, t_VEC);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    GEN p1 = gen_0, p2 = gen_0;	  /* vec */
    while (1)
    {
      p1 = cgetg(2, t_VEC);
      gel(p1, 1) = gcopy(gel(vec, 1));
      res = concat(res, p1);
      ind = gen_0;
      {
        pari_sp btop = avma, st_lim = stack_lim(btop, 1);
        do
        {
          ind = gaddgs(ind, 1);
          gel(vec, gtos(ind)) = gen_0;
          if (low_stack(st_lim, stack_lim(btop, 1)))
            gerepileall(btop, 2, &ind, &vec);
        } while(!(gequalgs(ind, glength(vec)) || (gcmp(gabs(gel(vec, gtos(gaddgs(ind, 1))), prec), gpow(stoi(10), gneg(asympdigits), prec)) > 0)));
      }
      if (gcmpgs(ind, glength(vec)) >= 0)
        break;
      p2 = cgetg(2, t_VEC);
      gel(p2, 1) = gcopy(ind);
      res = concat(res, p2);
      vec = gtovec(gdiv(gpow(x, ind, prec), gtoser(vec, -1, precdl)));
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 5, &p1, &res, &ind, &vec, &p2);
    }
  }
  res = gerepilecopy(ltop, res);
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/*                 EvaluateContFrac(cfvec,terms,t)                        \\ */
/* Evaluate a continued fraction at x=t, using a given number of terms    \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
EvaluateContFrac(GEN cfvec, GEN terms, GEN t, long prec)
{
  pari_sp ltop = avma;
  GEN res = gen_0;
  if ((gcmpgs(terms, 0) < 0) || (gcmpgs(terms, glength(cfvec)/2) > 0))
    terms = stoi(glength(cfvec)/2);
  res = gcopy(gel(cfvec, gtos(gaddgs(gmulsg(2, terms), 1))));
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    GEN p1 = gen_0;
    while (gcmpgs(terms, 0) > 0)
    {
      if (!gequal0(res))
        p1 = gadd(gel(cfvec, gtos(gsubgs(gmulsg(2, terms), 1))), gdiv(gpow(t, gel(cfvec, gtos(gmulsg(2, terms))), prec), res));
      else
        p1 = gpow(stoi(10), asympdigits, prec);
      res = p1;
      terms = gsubgs(terms, 1);
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 3, &p1, &res, &terms);
    }
  }
  res = gerepilecopy(ltop, res);
  return 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.                                                   \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

/* jdemeyer: second argument is never used, it gets overriden in the
* first line */
GEN
cflength(GEN cutoff, long prec)
{
  GEN expdifff = gen_0, asympconstf = gen_0, err = gen_0, t = gen_0, t1 = gen_0, t2 = gen_2, tt = gen_0, res = gen_0, lfundigits = gen_0, p1 = gen_0, p2 = gen_0;
  long d, k;
  GEN gd;
  if (!cutoff)
    cutoff = strtor("1.2", prec);
  vA = gsqrt(gdiv(conductor, gpowgs(mppi(prec), glength(gammaV))), prec);
  vA = gclone(vA);
  d = glength(gammaV); gd = stoi(d);
  lfundigits = gaddsg(getrealprecision(), gmaxgs(gceil(gdiv(gneg(glog(gabs(fullgamma(gadd(gmul(strtor("0.7", prec), weight), gmul(MaxImaginaryPart, gen_I())), prec), prec), prec)), glog(stoi(10), prec))), 0));
  lfundigits = gclone(lfundigits);
  p1 = gen_0;
  for (k = 1; k <= d; ++k) p1 = gadd(p1, gel(gammaV, k));
  expdifff = gclone(gsubgs(gdivgs(gaddgs(p1, 1), d), 1));
  p2 = gen_1;
  for (k = 1; k <= d; ++k) p2 = gmul(p2, ggamma(gdivsg(k, gd), prec));
  asympconstf = gclone(gmulsg(2, p2));
  err = gpow(stoi(10), gsub(gneg(lfundigits), strtor("0.5", prec)), prec);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    do
    {
      t = gequal0(t1) ? t2 : gdiventgs(gadd(t1, t2), 2);
      tt = gdiv(gdiv(t, cutoff), vA);
      res = gmul(gmul(gmul(coefgrow(t, prec), asympconstf), gexp(gmulsg(-d, gpow(tt, gdivsg(2, gd), prec)), prec)), gpow(tt, expdifff, prec));
      if (!gequal0(t1))
      {
        if (gcmp(res, err) > 0) t1 = t;
        else t2 = t;
      }
      else
      {
        if (gcmp(res, err) < 0) t1 = gdivgs(t2, 2);
        else t2 = gmulgs(t2, 2);
      }
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 5, &t, &tt, &res, &t1, &t2);
    } while(!(gcmpgs(gsub(t2, t1), 1) <= 0));
  }
  return gceil(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                                                         \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

long
initLdata(GEN vstr, GEN cutoff, GEN vdualstr, long prec)
{
  GEN pordtmp = gen_0, recF = gen_0, recG = gen_0, terms = gen_0, diff = gen_0, p1 = gen_0;
  GEN p2 = gen_0, p3 = gen_0, p4 = gen_0, p5 = gen_0, p6 = gen_0;	  /* vec */
  GEN recFG = pol_x(fetch_user_var("recFG")), p7 = gen_0, p8 = gen_0, p9 = gen_0, p10 = gen_0;
  GEN n = pol_x(fetch_user_var("n"));
  long p11;
  long p15;
  GEN p16 = gen_0;	  /* vec */
  long p17;
  GEN p18 = gen_0;	  /* vec */
  long l19;
  GEN t1 = pol_x(fetch_user_var("t1"));
  GEN p20 = gen_0;	  /* vec */
  GEN x = pol_x(fetch_user_var("x")), X = pol_x(fetch_user_var("X"));
  long l21;
  long d, len, j, k, sncoeff;
  GEN gd;
  if (!cutoff) cutoff = strtor("1.2", prec);
  if (!vdualstr) vdualstr = strtoGENstr("");
  if (typ(gammaV) != t_VEC)
    pari_err(user, "Gamma-factor gammaV has to be defined before calling initLdata()");
  if ((typ(conductor) != t_INT) && (typ(conductor) != t_REAL))
    pari_err(user, "conductor has to be defined before calling initLdata()");
  if (gcmpgs(cutoff, 0) < 0) p1 = gneg(cutoff);
  else p1 = cflength(cutoff, prec);
  len = gtos(p1);
  p2 = cgetg(len+1, t_VEC);
  for (k = 1; k <= len; ++k)
    gel(p2, k) = closure_callgen1(vstr, stoi(k));
  cfvec = gclone(p2);
  if (gequal(vdualstr, strtoGENstr(""))) cfdualvec = cfvec;
  else
  {
    p3 = cgetg(len+1, t_VEC);
    for (k = 1; k <= len; ++k)
      gel(p3, k) = closure_callgen1(vdualstr, stoi(k));
    cfdualvec = gclone(p3);
  }
  if (gcmpgs(cutoff, 0) < 0) len = gtos(cflength(NULL, prec));
  lastFGs = gclone(negr(strtor("1E99", prec)));
  /* globals to track what was calculated last */
  lastFGSs = gclone(negr(strtor("1E99", prec)));
  p4 = cgetg(3, t_VEC);
  gel(p4, 1) = negr(strtor("1E99", prec));
  gel(p4, 2) = negr(strtor("1E99", prec));
  /* to avoid re-calculating same values each time */
  lastGs = gclone(p4);
  p5 = cgetg(3, t_VEC);
  gel(p5, 1) = negr(strtor("1E99", prec));
  gel(p5, 2) = negr(strtor("1E99", prec));
  lastLSs = gclone(p5);
  d = glength(gammaV); gd =stoi(d);
  /* d = number of gamma factors */
  
  /* Calculate the necessary amount of extra digits */
  
  answerdigits = gclone(stoi(getrealprecision()));
  vA = gsqrt(gdiv(conductor, gpowgs(mppi(prec), d)), prec);
  vA = gclone(vA);
  lfundigits = gadd(answerdigits, gmaxgs(gceil(gdiv(gneg(glog(gabs(fullgamma(gadd(gmul(strtor("0.7", prec), weight), gmul(MaxImaginaryPart, gen_I())), prec), prec), prec)), glog(stoi(10), prec))), 0));
  lfundigits = gclone(lfundigits);
  termdigits = gclone(gadd(lfundigits, gfloor(gsubgs(weight, 1))));
  taylordigits = gclone(gmulsg(2, termdigits));
  asympdigits = termdigits;
  /* Exponential factor defined to maximal precision */
  
  setrealprecision(gtos(taylordigits), &prec);
  vA = gsqrt(gdiv(precision0(conductor, gtos(taylordigits)), gpowgs(mppi(prec), d)), prec);
  vA = gclone(vA);
  /* exp. factor */
  lastt = gclone(gdiv(stoi(len), vA));
  p6 = cgetg(d+1, t_VEC);
  for (k = 1; k <= d; ++k) gel(p6, k) = gen_1;
  pordtmp = p6;
  /* locate poles and their orders */
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    for (j = 1; j <= d; ++j)
    {
      {
        for (k = 1; k <= d; ++k)
        {
          if (j != k)
            if (((typ(diff = gsub(gel(gammaV, j), gel(gammaV, k))) == t_INT) && gequal0(gmodgs(diff, 2))) && (gcmpgs(diff, 0) <= 0))
            {
              gel(pordtmp, j) = gadd(gel(pordtmp, j), gel(pordtmp, k));
              gel(pordtmp, k) = gen_0;
            }
        }
      }
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 2, &diff, &pordtmp);
    }
  }
  poles = cgetg(1, t_VEC);
  PoleOrders = cgetg(1, t_VEC);
  for (j = 1; j <= d; ++j)
  {
    if (!gequalgs(gel(pordtmp, j), 0))
    {
      poles = concat(poles, gel(gammaV, j));
      PoleOrders = concat(PoleOrders, gel(pordtmp, j));
    }
  }
  poles = gclone(poles); PoleOrders = gclone(PoleOrders);
  numpoles = gclone(stoi(glength(poles)));
  /* Initialize the asymptotic coefficients at infinity */
  
  setrealprecision(gtos(asympdigits), &prec);
  recFG = RecursionsAtInfinity(prec);
  recF = gcopy(gel(recFG, 1));
  recG = gcopy(gel(recFG, 2));
  kill0("recFG");
  /* Maximum number of terms in the asymptotic expansion */
  
  ncoeff = MaxAsympCoeffs;
  p7 = gen_0;
  for (k = 1; k <= d; ++k) p7 = gadd(p7, gel(gammaV, k));
  /* Asymptotic behaviour at infinity for phi and G */
  
  expdifff = gclone(gsubgs(gdivgs(gaddgs(p7, 1), d), 1));
  p8 = gen_1;
  for (k = 1; k <= d; ++k) p8 = gmul(p8, ggamma(gdivsg(k, gd), prec));
  asympconstf = gclone(gmulsg(2, p8));
  p9 = gen_0;
  for (k = 1; k <= d; ++k) p9 = gadd(p9, gel(gammaV, k));
  expdiffg = gclone(gsub(gsubgs(gdivgs(gaddgs(p9, 1), d), 1), gdivsg(2, gd)));
  p10 = gen_1;
  for (k = 1; k <= d; ++k) p10 = gmul(p10, ggamma(gdivsg(k, gd), prec));
  asympconstg = gclone(p10);
  sncoeff = gtos(ncoeff);
  p11 = d + sncoeff;
  Fvec = cgetg(p11+1, t_VEC);
  for (j = 1; j <= p11; ++j) gel(Fvec, j) = gen_0;
  /* Coefficients for the asymptotic expansion of phi(t) and G(t) */
  gel(Fvec, d) = real_1(prec);
  Gvec = cgetg(p11+1, t_VEC);
  for (j = 1; j <= p11; ++j) gel(Gvec, j) = gen_0;
  gel(Gvec, d) = real_1(prec);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    GEN p23, p25, gy;
    long y;
    for (y = 1; y <= sncoeff; ++y)
    {
      gy=stoi(y);
      {
        long l24;
        p23 = p25 = gen_0;
        for (j = 1; j <= d; ++j)
        {
          l24 = d + y - j;
	  if (j < d)
            p23 = gadd(p23, gmul(gsubst(gel(recF, j), gvar(n), gy), gel(Fvec, l24)));
	  p25 = gadd(p25, gmul(gsubst(gel(recG, j), gvar(n), gy), gel(Gvec, l24)));
        }
      }
      gel(Fvec, d + y) = gmul(real_1(prec), p23);
      gel(Gvec, d + y) = gmul(real_1(prec), p25);
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 2, &Gvec, &Fvec);
    }
  }
  Fvec = gclone(Fvec); Gvec = gclone(Gvec);
  p15 = sncoeff + 1;
  p16 = cgetg(p15+1, t_VEC);
  for (k = 1; k <= p15; ++k) gel(p16, k) = gel(Fvec, d + k - 1);
  /* Convert the Fvec (Taylor asymptotic) coefficients into fcf (contfrac coeffs) */
  
  fcf = gclone(SeriesToContFrac(p16, prec));
  fncf = gclone(stoi(glength(fcf)/2));
  /* at most ncoeff+1, less if terminates */
  
  /* Taylor series coefficients of phi(t) around t=infinity */
  
  if (gcmpgs(lastt, 35) < 0) termstep = gen_1;
  else termstep = gclone(gfloor(gpow(lastt, ginv(stoi(3)), prec)));
  p17 = gtos(ground(gdiv(lastt, termstep))) + 1;
  p18 = cgetg(p17+1, t_VEC);
  for (k = 1; k <= p17; ++k) gel(p18, k) = gen_m1;
  phiinfterms = gclone(p18);
  terms = fncf;
  PhiCaseBound = gen_0;
  l19 = glength(phiinfterms);
  {
    GEN p27 = gen_0;
    for (k = 1; k <= l19; ++k)
    {
      t1 = gmulsg(k - 1, termstep);
      {
        pari_sp btop = avma, st_lim = stack_lim(btop, 1);
        while (((k > 1) && (gcmpgs(terms, 1) > 0)) && (gcmp(gabs(gsub(phiinf(t1, gsubgs(terms, 1), prec), phiinf(t1, terms, prec)), prec), gpow(stoi(10), gsubgs(gneg(termdigits), 1), prec)) < 0))
        {
          terms = gsubgs(terms, 1);
          if (low_stack(st_lim, stack_lim(btop, 1)))
            terms = gerepilecopy(btop, terms);
        }
      }
      p27 = gen_0;
      for (j = 1; gcmpsg(j, terms) <= 0; ++j)
        p27 = gadd(p27, gel(fcf, 2*j));
      if (gcmp(p27, ncoeff) < 0) gel(phiinfterms, k) = gcopy(terms);
      else PhiCaseBound = gclone(gmulsg(k, termstep));
    }
  }
  /* Recursions for phi(t) and G(t,s) at the origin */
  phiinfterms = gclone(phiinfterms);
  setrealprecision(gtos(taylordigits), &prec);
  {
    GEN p29 = gen_0;
    p20 = cgetg(gtos(numpoles)+1, t_VEC);
    for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
    {
      p29 = gen_1;
      for (k = 1; k <= d; ++k)
        p29 = gmul(p29, gsubst(gammaseries(gdivgs(gadd(gneg(gel(poles, j)), gel(gammaV, k)), 2), gsubgs(gel(PoleOrders, j), 1), prec), gvar(x), gdivgs(X, 2)));
      gel(p20, j) = p29;
    }
  }
  /* Initial values of the gamma factors for recursions */
  
  InitV = p20;
  /* Taylor series coefficients of phi(t) around t=0 -> phiVser */
  
  phiV = cgetg(1, t_VEC);
  phiVnn = gen_0;
  phiVser = gcopy(InitV);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    do
    {
      RecursephiV();
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 3, &phiVnn, &phiV, &phiVser);
    } while(!((gcmpgs(phiVnn, 3) > 0) && (gcmp(gmul(vecmax(gabs(gel(phiV, gtos(phiVnn)), prec)), gpow(gmul(gaddgs(PhiCaseBound, 1), vA), gmulsg(2, phiVnn), prec)), gpow(stoi(10), gsubgs(gneg(termdigits), 1), prec)) < 0)));
  }
  l21 = setrealprecision(gtos(answerdigits), &prec);
  phiV = gclone(phiV); phiVnn = gclone(phiVnn); phiVser = gclone(phiVser);
  return l21;
}

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

GEN
phi(GEN t, long prec)
{
  pari_sp ltop = avma;
  GEN p1;
  if (gcmp(t, PhiCaseBound) < 0) p1 = phi0(t, prec);
  else
    p1 = phiinf(t, gel(phiinfterms, gtos(gmings(gaddsg(1, gfloor(gdiv(gabs(t, prec), termstep))), glength(phiinfterms)))), prec);
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

void
RecursephiV(void)	  /* void */
{
  pari_sp ltop = avma;
  GEN p1 = gen_0;	  /* vec */
  GEN p2 = gen_0;
  GEN p3 = gen_0;	  /* vec */
  GEN X = pol_x(fetch_user_var("X"));
  phiVnn = gaddgs(phiVnn, 1);
  p1 = cgetg(2, t_VEC);
  p2 = vecmax(PoleOrders);
  {
    long j, k;
    p3 = cgetg(gtos(p2)+1, t_MAT);
    for (k = 1; gcmpsg(k, p2) <= 0; ++k)
    {
      gel(p3, k) = cgetg(gtos(numpoles)+1, t_COL);
      for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
        gcoeff(p3, j, k) = polcoeff0(gel(phiVser, j), -k, -1);
    }
  }
  gel(p1, 1) = p3;
  phiV = concat(phiV, p1);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    long j;
    long l4;
    for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
    {
      l4 = glength(gammaV);
      {
        pari_sp btop = avma, st_lim = stack_lim(btop, 1);
        long k;
        for (k = 1; k <= l4; ++k)
        {
          gel(phiVser, j) = gdiv(gel(phiVser, j), gadd(gsub(gsub(gdivgs(X, 2), gdivgs(gel(poles, j), 2)), phiVnn), gdivgs(gel(gammaV, k), 2)));
          if (low_stack(st_lim, stack_lim(btop, 1)))
            phiVser = gerepilecopy(btop, phiVser);
        }
      }
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 1, &phiVser);
    }
  }
  gerepileall(ltop, 3, &phiVnn, &phiV, &phiVser);
  return;
}

GEN
phi0(GEN t, long prec)
{
  pari_sp ltop = avma;
  GEN t2 = gen_0, LogTTerm = gen_0, TPower = gen_0, res = gen_0, nn = gen_0, totalold = gen_0, p1 = gen_0;
  GEN p2 = gen_0, p3 = gen_0;	  /* vec */
  setrealprecision(gtos(taylordigits), &prec);
  t = precision0(t, gtos(taylordigits));
  t2 = gsqr(t);
  p1 = vecmax(PoleOrders);
  {
    long k;
    p2 = cgetg(gtos(p1)+1, t_VEC);
    for (k = 1; gcmpsg(k, p1) <= 0; ++k)
      gel(p2, k) = gdiv(gpowgs(gneg(glog(t, prec)), k - 1), mpfact(k - 1));
  }
  LogTTerm = gtrans(p2);
  {
    long j;
    p3 = cgetg(gtos(numpoles)+1, t_VEC);
    for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
      gel(p3, j) = gpow(t, gel(poles, j), prec);
  }
  TPower = gmul(real_1(prec), p3);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    do
    {
      totalold = res;
      nn = gaddgs(nn, 1);
      res = gadd(res, gmul(gmul(TPower, gel(phiV, gtos(nn))), LogTTerm));
      TPower = gmul(TPower, t2);
      if (low_stack(st_lim, stack_lim(btop, 1)))
        gerepileall(btop, 4, &totalold, &nn, &res, &TPower);
    } while(!((gcmp(gabs(gsub(res, totalold), prec), gpow(stoi(10), gneg(gaddgs(termdigits, 1)), prec)) < 0) && (gcmpgs(nn, 3) > 0)));
  }
  setrealprecision(gtos(termdigits), &prec);
  res = gerepilecopy(ltop, res);
  return res;
}

/* phi(t) using asymptotic expansion at infinity */
GEN
phiinf(GEN t, GEN ncf, long prec)
{
  pari_sp ltop = avma;
  GEN res = gen_0, d = gen_0, td2 = gen_0;
  if (!ncf) ncf = gcopy(fncf);
  setrealprecision(gtos(asympdigits), &prec);
  t = precision0(t, gtos(asympdigits));
  d = stoi(glength(gammaV));
  td2 = gpow(t, gdivsg(-2, d), prec);
  res = EvaluateContFrac(fcf, ncf, td2, prec);
  res = gmul(gmul(gmul(res, asympconstf), gexp(gdiv(gneg(d), td2), prec)), gpow(t, expdifff, prec));
  setrealprecision(gtos(termdigits), &prec);
  res = gerepilecopy(ltop, res);
  return 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)                             \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
G(GEN t, GEN ss, GEN der, long prec)
{
  GEN nn = gen_0;
  GEN p1 = gen_0, p2 = gen_0;	  /* vec */
  GEN p3 = gen_0;
  if (!der) der = gen_0;
  p1 = cgetg(3, t_VEC);
  gel(p1, 1) = gcopy(ss);
  gel(p1, 2) = gcopy(der);
  if (!gequal(lastGs, p1) || (typ(Ginfterms) != t_VEC))
  {
    initGinf(ss, der, prec);
    p2 = cgetg(3, t_VEC);
    gel(p2, 1) = gcopy(ss);
    gel(p2, 2) = gcopy(der);
    lastGs = gclone(p2);
  }
  if (gcmp(t, GCaseBound) < 0) p3 = G0(t, ss, der, prec);
  else
  {
    nn = gmings(gaddsg(1, gfloor(gdiv(gabs(t, prec), termstep))), glength(Ginfterms));
    p3 = Ginf(t, ss, der, gel(Ginfterms, gtos(nn)), prec);
  }
  return p3;
}

GEN
LogInt(GEN i, GEN j, GEN logt, GEN der, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0, p2 = gen_0, p3 = gen_0;
  long sder = gtos(der);
  if (gcmp(gabs(i, prec), gpow(stoi(10), gsubsg(2, termdigits), prec)) < 0)
    p1 = gen_0;
  else
  {
    p2 = gsubgs(j, 1);
    {
      long k;
      p3 = gen_0;
      for (k = 0; gcmpsg(k, p2) <= 0; ++k)
        p3 = gadd(p3, gdiv(gmul(gmul(binomial(gsubsg(k, j), sder), mpfact(sder)), gpowgs(gmul(gneg(i), logt), k)), mpfact(k)));
    }
    p1 = gdiv(p3, gpow(i, gadd(j, der), prec));
  }
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

GEN
MakeLogSum(GEN ss, GEN der, long prec)
{
  GEN V = gen_0, logsum = gen_0, p1 = gen_0;
  long l2;
  GEN lt = pol_x(fetch_user_var("lt"));
  GEN p3 = gen_0;	  /* vec */
  long j;
  if (gequalsg(glength(LogSum), phiVnn))
  {
    p1 = gaddgs(gfloor(gdivgs(phiVnn, 10)), 1);
    for (j = 1; gcmpsg(j, p1) <= 0; ++j) RecursephiV();
  }
  phiV = gclone(phiV); phiVnn = gclone(phiVnn); phiVser = gclone(phiVser);
  l2 = glength(LogSum) + 1;
  {
    long nn;
    GEN p5 = gen_0, p6 = gen_0;	  /* vec */
    for (nn = l2; gcmpsg(nn, phiVnn) <= 0; ++nn)
    {
      /* generate logsums */
      V = gcopy(gel(phiV, nn));
      {
        long j, p7;
        GEN p8 = gen_0;
        p5 = cgetg(gtos(numpoles)+1, t_VEC);
        for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
        {
          p7 = gtos(gel(PoleOrders, j));
          {
	    long k;
            p8 = gen_0;
            for (k = 1; k <= p7; ++k)
              p8 = gadd(p8, gmul(gcoeff(V, j, k), LogInt(gadd(gaddgs(gel(poles, j), 2*(nn - 1)), ss), stoi(k), lt, der, prec)));
          }
          gel(p5, j) = p8;
        }
      }
      logsum = gtrans(p5);
      p6 = cgetg(2, t_VEC); gel(p6, 1) = gcopy(logsum);
      LogSum = concat(LogSum, p6);
    }
  }
  LogSum = gclone(LogSum);
  p3 = cgetg(3, t_VEC);
  gel(p3, 1) = gcopy(ss);
  gel(p3, 2) = gcopy(der);
  lastLSs = gclone(p3);
  return lastLSs;
}

GEN
G0(GEN t, GEN ss, GEN der, long prec)
{
  GEN t2 = gen_0, LT = gen_0, TPower = gen_0, res = gen_0, term = gen_0, gmser = gen_0, gmcf = gen_0;
  GEN p1 = gen_0, p2 = gen_0;	  /* vec */
  GEN lt = pol_x(fetch_user_var("lt")), S = pol_x(fetch_user_var("S"));
  long nn;
  setrealprecision(gtos(taylordigits), &prec);
  ss = precision0(ss, gtos(taylordigits));
  p1 = cgetg(3, t_VEC);
  gel(p1, 1) = gcopy(ss);
  gel(p1, 2) = gcopy(der);
  if (!gequal(p1, lastLSs)) LogSum = cgetg(1, t_VEC);
  t = precision0(t, gtos(taylordigits));
  t2 = gsqr(t);
  /* time */
  LT = glog(t, prec);
  {
    long j;
    p2 = cgetg(gtos(numpoles)+1, t_VEC);
    for (j = 1; gcmpsg(j, numpoles) <= 0; ++j)
      gel(p2, j) = gpow(t, gel(poles, j), prec);
  }
  /* = money */
  TPower = p2;
  res = gen_0;
  nn = 0;
  term = gen_1;
  {
    do
    {
      ++nn;
      if (nn > glength(LogSum)) MakeLogSum(ss, der, prec);
      term = gmul(TPower, gsubst(gel(LogSum, nn), gvar(lt), LT));
      res = gadd(res, term);
      TPower = gmul(TPower, t2);
    } while(!((nn > 3) && (gcmp(gabs(term, prec), gpow(stoi(10), gneg(gaddgs(termdigits, 1)), prec)) < 0)));
  }
  gmser = gdiv(fullgammaseries(ss, der, prec), gpow(t, gadd(S, ss), prec));
  gmcf = gmul(polcoeff0(gmser, gtos(der), gvar(S)), mpfact(gtos(der)));
  res = gsub(gmcf, res);
  setrealprecision(gtos(termdigits), &prec);
  return res;
}

GEN
Ginf(GEN t, GEN ss, GEN der, GEN ncf, long prec)
{
  pari_sp ltop = avma;
  GEN res = gen_0, d = gen_0, tt = gen_0;
  if (!ncf)
    ncf = gen_m1;
  /* at infinity and associated continued fraction */
  setrealprecision(gtos(asympdigits), &prec);
  ss = precision0(ss, gtos(asympdigits));
  t = precision0(t, gtos(asympdigits));
  if (gequalm1(ncf)) ncf = gcopy(gncf);
  d = stoi(glength(gammaV));
  tt = gpow(t, gdivsg(-2, d), prec);
  res = EvaluateContFrac(gcf_T, ncf, tt, prec);
  res = gmul(gmul(gmul(gmul(asympconstg, gexp(gdiv(gneg(d), tt), prec)), gpow(t, expdiffg, prec)), gpow(tt, der, prec)), res);
  setrealprecision(gtos(termdigits), &prec);
  res = gerepilecopy(ltop, res);
  return res;
}

long
initGinf(GEN ss, GEN der, long prec)
{
  GEN gvec = gen_0, gncf = gen_0, terms = gen_0, t1 = gen_0, s = pol_x(fetch_user_var("s")), p1 = gen_0;
  GEN p2 = gen_0;	  /* vec */
  GEN p3 = gen_0;
  GEN p4 = gen_0;	  /* vec */
  long l5, l6;
  long d, k;
  setrealprecision(gtos(asympdigits), &prec);
  ss = precision0(ss, gtos(asympdigits));
  d = glength(gammaV);
  gvec = gcopy(Gvec);
  for (k = 1; gcmpsg(k, der) <= 0; ++k)
  {
    gvec = deriv(gvec, gvar(s));
    gvec = concat(vecright(gvec, stoi(glength(gvec) - 1)), gen_1);
  }
  p1 = gaddgs(ncoeff, 1);
  {
    long k;
    p2 = cgetg(gtos(p1)+1, t_VEC);
    for (k = 1; gcmpsg(k, p1) <= 0; ++k)
      gel(p2, k) = gsubst(gel(gvec, d + k - 1), gvar(s), ss);
  }
  gcf_T = gclone(SeriesToContFrac(p2, prec));
  gncf = gclone(stoi(glength(gcf_T)/2));
  p3 = gaddgs(ground(gdiv(lastt, termstep)), 1);
  p4 = cgetg(gtos(p3)+1, t_VEC);
  for (k = 1; gcmpsg(k, p3) <= 0; ++k) gel(p4, k) = gen_m1;
  Ginfterms = p4;
  terms = gncf;
  GCaseBound = gen_0;
  l5 = glength(Ginfterms);
  {
    GEN p7 = gen_0;
    for (k = 1; k <= l5; ++k)
    {
      t1 = gmulsg(k - 1, termstep);
      {
        pari_sp btop = avma, st_lim = stack_lim(btop, 1);
        while (((k > 1) && (gcmpgs(terms, 1) > 0)) && (gcmp(gabs(gsub(Ginf(t1, ss, der, gsubgs(terms, 1), prec), Ginf(t1, ss, der, terms, prec)), prec), gpow(stoi(10), gsubgs(gneg(termdigits), 2), prec)) < 0))
        {
          terms = gsubgs(terms, 1);
          if (low_stack(st_lim, stack_lim(btop, 1)))
            terms = gerepilecopy(btop, terms);
        }
      }
      long j;
      p7 = gen_0;
      for (j = 1; gcmpsg(j, terms) <= 0; ++j) p7 = gadd(p7, gel(gcf_T, 2*j));
      if (gcmp(p7, ncoeff) < 0) gel(Ginfterms, k) = gcopy(terms);
      else GCaseBound = gclone(gmulsg(k, termstep));
    }
  }
  Ginfterms = gclone(Ginfterms);
  l6 = setrealprecision(gtos(termdigits), &prec);
  return l6;
}

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

GEN
ltheta(GEN t, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0, p2;
  GEN p3 = gen_0;
  long k;
  t = precision0(gdiv(t, vA), gtos(taylordigits));
  p1 = gmings(gfloor(gaddgs(gdiv(lastt, t), 1)), glength(cfvec));
  p2 = gen_0;
  for (k = 1; gcmpsg(k, p1) <= 0; ++k)
  {
    if (!gequal0(gel(cfvec, k))) p3 = phi(gmulsg(k, t), prec);
    else p3 = gen_0;
    p2 =gadd(p2, gmul(gel(cfvec, k), p3));
  }
  p2 = gerepilecopy(ltop, p2);
  return p2;
}

GEN
ldualtheta(GEN t, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0, p2;
  GEN p3 = gen_0;
  long k;
  t = precision0(gdiv(t, vA), gtos(taylordigits));
  p1 = gmings(gfloor(gaddgs(gdiv(lastt, t), 1)), glength(cfvec));
  p2 = gen_0;
  for (k = 1; gcmpsg(k, p1) <= 0; ++k)
  {
    if (!gequal0(gel(cfdualvec, k))) p3 = phi(gmulsg(k, t), prec);
    else p3 = gen_0;
    p2 = gadd(p2, gmul(gel(cfdualvec, k), p3));
  }
  p2 = gerepilecopy(ltop, p2);
  return p2;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* 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             \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
checkfeq(GEN t, long prec)
{
  GEN lpx = gen_0, lpv = gen_0, lpm = gen_0, res = gen_0, tq = gen_0, ss = gen_0, automatic = pol_x(fetch_user_var("automatic"));
  GEN p1 = gen_0, p2 = gen_0, p3 = gen_0;	  /* vec */
  long l4, nlp;
  GEN p5 = gen_0;
  if (!t) t = gdivgs(stoi(6), 5);
  /* .. and check the functional equation */
  if (gequal(Lresidues, automatic) && !gequal(Lpoles, cgetg(1, t_VEC)))
  {
    nlp = glength(Lpoles);
    {
      long k;
      p1 = cgetg(nlp+1, t_VEC);
      for (k = 1; k <= nlp; ++k)
        gel(p1, k) = gadd(strtor("1.15", prec), gdivgs(stoi(k - 1), 10));
    }
    lpx = p1;
    {
      long k;
      p2 = cgetg(nlp+1, t_VEC);
      for (k = 1; k <= nlp; ++k)
      {
        tq = gcopy(gel(lpx, k));
        gel(p2, k) = gsub(gmul(gmul(sgn, gpow(tq, weight, prec)), ldualtheta(tq, prec)), ltheta(ginv(tq), prec));
      }
    }
    lpv = p2;
    {
      long k, j;
      p3 = cgetg(nlp+1, t_MAT);
      for (j = 1; j <= nlp; ++j)
      {
        gel(p3, j) = cgetg(nlp+1, t_COL);
        for (k = 1; k <= nlp; ++k)
        {
          tq = gcopy(gel(lpx, k));
          ss = gcopy(gel(Lpoles, j));
          gcoeff(p3, k, j) = gsub(gpow(tq, ss, prec), gmul(sgn, gpow(tq, gsub(weight, ss), prec)));
        }
      }
    }
    lpm = p3;
    Lresidues = gclone(gtrans(gauss(lpm, gtrans(lpv))));
  }
  l4 = glength(Lpoles);
  {
    long k;
    p5 = gen_0;
    for (k = 1; k <= l4; ++k)
    {
      p5 = gadd(p5, gmul(gel(Lresidues, k), gsub(gpow(t, gneg(gel(Lpoles, k)), prec), gmul(sgn, gpow(t, gadd(gneg(weight), gel(Lpoles, k)), prec)))));
    }
  }
  /*for (k=1,nlp,print("Residue at ",Lpoles[k]," = ",Lresidues[k])); */
  res = gadd(gsub(ltheta(t, prec), gmul(gmul(sgn, gpow(t, gneg(weight), prec)), ldualtheta(ginv(t), prec))), p5);
  setrealprecision(gtos(answerdigits), &prec);
  return 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                             \\ */
/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */

GEN
L(GEN ss, GEN cutoff, GEN der, long prec)
{
  GEN p1 = gen_0;
  if (!cutoff) cutoff = gen_1;
  if (!der) der = gen_0;
  p1 = gmul(polcoeff0(Lseries(ss, cutoff, der, prec), gtos(der), -1), mpfact(gtos(der)));
  return p1;
}

/* Lseries(s,,der) = L(s)+L'(s)S+L''(s)*S^2/2!+... ,first der terms */
GEN
Lseries(GEN ss, GEN cutoff, GEN der, long prec)
{
  GEN FGSeries = gen_0, res = gen_0, LSSeries = gen_0, PrecisionLoss = gen_0, S = pol_x(fetch_user_var("S"));
  long l1;
  GEN p2 = gen_0;	  /* vec */
  long l3;
  GEN p4 = gen_0;	  /* vec */
  GEN p5 = gen_0;
  if (!cutoff) cutoff = gen_1;
  if (!der) der = gen_0;
  setrealprecision(gtos(lfundigits), &prec);
  FGSeries = gmul(fullgammaseries(ss, der, prec), gpow(vA, gadd(ss, S), prec));
  l1 = glength(Lpoles);
  {
    long k;
    p2 = cgetg(l1+1, t_VEC);
    for (k = 1; k <= l1; ++k)
      gel(p2, k) = gsub(gel(Lpoles, k), ss);
  }
  l3 = glength(Lpoles);
  {
    long k;
    p4 = cgetg(l3+1, t_VEC);
    for (k = 1; k <= l3; ++k)
      gel(p4, k) = gsub(gsub(weight, gel(Lpoles, k)), ss);
  }
  if (glength(Lpoles) && ((gcmp(vecmin(gabs(p2, prec)), gpow(stoi(10), gneg(answerdigits), prec)) < 0) || (gcmp(vecmin(gabs(p4, prec)), gpow(stoi(10), gneg(answerdigits), prec)) < 0)))
    pari_err(user, "L*(s) has a pole at s=%Ps", 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 = gsub(gceil(gdiv(gneg(glog(vecmax(gabs(gtovec(FGSeries), prec)), prec)), glog(stoi(10), prec))), gsub(lfundigits, answerdigits));
  if (gcmpgs(PrecisionLoss, 1) > 0)
    pari_printf("Warning: Loss of %Ps digits due to cancellation\n", PrecisionLoss);
  {
    GEN k = gen_0;
    p5 = gen_0;
    for (k = gen_0; gcmp(k, der) <= 0; k = gaddgs(k, 1))
    {
      p5 = gadd(p5, gdiv(gmul(Lstar(ss, cutoff, k, prec), gpow(S, k, prec)), mpfact(gtos(k))));
    }
  }
  LSSeries = gadd(p5, ggrando(S, gtos(gaddgs(der, 1))));
  res = gdiv(LSSeries, FGSeries);
  setrealprecision(gtos(answerdigits), &prec);
  return res;
}

GEN
Lstar(GEN ss, GEN cutoff, GEN der, long prec)
{
  GEN res = gen_0, ncf1 = gen_0, ncf2 = gen_0;
  long l1;
  GEN p2 = gen_0;
  long l3;
  GEN p4 = gen_0, p5 = gen_0, p6 = gen_0;
  if (!cutoff) cutoff = gen_1;
  if (!der) der = gen_0;
  if (!gequal0(der) && !gequalgs(cutoff, 1))
    pari_err(user, "L(s,cutoff,k>0) is only implemented for cutoff=1");
  ss = precision0(ss, gtos(taylordigits));
  cutoff = precision0(cutoff, gtos(taylordigits));
  ncf1 = gmings(ground(gmul(gmul(lastt, vA), cutoff)), glength(cfvec));
  ncf2 = gmings(ground(gdiv(gmul(lastt, vA), cutoff)), glength(cfvec));
  setrealprecision(gtos(termdigits), &prec);
  l1 = glength(Lpoles);
  {
    pari_sp btop = avma, st_lim = stack_lim(btop, 1);
    long k;
    p2 = gen_0;
    for (k = 1; k <= l1; ++k)
    {
      p2 = gadd(p2, gmul(gdiv(gmul(gmul(gpow(gen_m1, der, prec), mpfact(gtos(der))), gel(Lresidues, k)), gpow(gsub(ss, gel(Lpoles, k)), gaddgs(der, 1), prec)), gpow(cutoff, gneg(gel(Lpoles, k)), prec)));
      if (low_stack(st_lim, stack_lim(btop, 1)))
        p2 = gerepilecopy(btop, p2);
    }
  }
  l3 = glength(Lpoles);
  {
    long k;
    p4 = gen_0;
    for (k = 1; k <= l3; ++k)
    {
      p4 = gadd(p4, gmul(gdiv(gmul(gel(Lresidues, k), mpfact(gtos(der))), gpow(gsub(gsub(weight, gel(Lpoles, k)), ss), gaddgs(der, 1), prec)), gpow(cutoff, gadd(gneg(weight), gel(Lpoles, k)), prec)));
    }
  }
  {
    GEN k = gen_0, p7 = gen_0;
    p5 = gen_0;
    for (k = gen_1; gcmp(k, ncf1) <= 0; k = gaddgs(k, 1))
    {
      if (!gequal0(gel(cfdualvec, gtos(k))))
        p7 = gmul(gmul(gel(cfdualvec, gtos(k)), gpow(gen_m1, der, prec)), G(gdiv(gdiv(k, vA), cutoff), gsub(weight, ss), der, prec));
      else p7 = gen_0;
      p5 = gadd(p5, gdiv(p7, gpow(cutoff, weight, prec)));
    }
  }
  {
    GEN p8 = gen_0;
    long k;
    p6 = gen_0;
    for (k = 1; gcmpsg(k, ncf2) <= 0; ++k)
    {
      if (!gequal0(gel(cfvec, k)))
        p8 = gmul(gel(cfvec, k), G(gdiv(gmulsg(k, cutoff), vA), ss, der, prec));
      else p8 = gen_0;
      p6 = gadd(p6, p8);
    }
  }
  res = gmul(gadd(gadd(gsub(gneg(p2), gmul(sgn, p4)), gmul(sgn, p5)), p6), gpow(cutoff, ss, prec));
  return res;
}

