/********************************************************************/
/**                                                                **/
/**                   TRANSCENDENTAL FONCTIONS                     **/
/**                          (part 3)                              **/
/**                                                                **/
/********************************************************************/
/* $Id: trans3.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
#include "pari.h"

/***********************************************************************/
/**                                                                   **/
/**                       K BESSEL FUNCTION                           **/
/**                                                                   **/
/***********************************************************************/

GEN
kbessel0(GEN nu, GEN gx, long flag, long prec)
{
  switch(flag)
  {
    case 0: return kbessel(nu,gx,prec);
    case 1: return kbessel2(nu,gx,prec);
    default: err(flagerr,"besselk");
  }
  return NULL; /* not reached */
}

GEN
kbessel(GEN nu, GEN gx, long prec)
{
  GEN x,y,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w;
  long l,lbin,av,av1,k,k2,l1,n2,n;

  if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
  if (typ(gx)!=t_REAL)
    { l=prec; k=1; }
  else
    { l=lg(gx); k=0; x=gx; }
  y=cgetr(l); l1=l+1;
  av=avma; if (k) gaffect(gx,x=cgetr(l));
  u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1);
  e=cgetr(l1); f=cgetr(l1);
  nu2=gmulgs(gsqr(nu),-4);
  n = (long) (bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gnorm(nu)))) / 2;
  n2=(n<<1); pitemp=mppi(l1);
  /* this 10 should really be a 5, but then kbessel(15.99) enters oo loop */
  lbin = 10 - bit_accuracy(l); av1=avma;
  if (gcmpgs(x,n)<0)
  {
    zf=gsqrt(gdivgs(pitemp,n2),prec);
    zz=cgetr(l1); gaffect(ginv(stoi(n2<<2)), zz);
    s=gun; t=gzero;
    for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
    {
      if (k2 & ~(MAXHALFULONG>>1))
        p1 = gadd(mulss(k2,k2),nu2);
      else
        p1 = gaddsg(k2*k2,nu2);
      ak=gdivgs(gmul(p1,zz),-k);
      s=gadd(gun,gmul(ak,s));
      t=gaddsg(k2,gmul(ak,t));
    }
    gmulz(s,zf,u); t=gmul2n(t,-1);
    gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1;
    affsr(n2,q=cgetr(l1));
    r=gmul2n(x,1); av1=avma;
    for(;;)
    {
      p1=divsr(5,q); if (expo(p1) >= -1) p1=ghalf;
      p2=subsr(1,divrr(r,q));
      if (gcmp(p1,p2)>0) p1=p2;
      gnegz(p1,c); gaffsg(1,d); affrr(u,e); affrr(v,f);
      for (k=1; ; k++)
      {
	w=gadd(gmul(gsubsg(k,ghalf),u), gmul(subrs(q,k),v));
	w=gadd(w, gmul(nu,gsub(u,gmul2n(v,1))));
	gdivgsz(gmul(q,v),k,u);
	gdivgsz(w,k,v);
	gmulz(d,c,d);
	gaddz(e,gmul(d,u),e); p1=gmul(d,v);
	gaddz(f,p1,f);
	if (gexpo(p1)-gexpo(f) <= lbin) break;
	avma=av1;
      }
      p1=u; u=e; e=p1;
      p1=v; v=f; f=p1;
      gmulz(q,gaddsg(1,c),q);
      if (expo(subrr(q,r)) <= lbin) break;
    }
    gmulz(u,gpui(gdivgs(x,n),nu,prec),y);
  }
  else
  {
    p2=gmul2n(x,1);
    zf=gsqrt(gdiv(pitemp,p2),prec);
    zz=ginv(gmul2n(p2,2)); s=gun;
    for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
    {
      if (k2 & ~(MAXHALFULONG>>1))
        p1=gadd(mulss(k2,k2),nu2);
      else
        p1=gaddsg(k2*k2,nu2);
      ak=gdivgs(gmul(p1,zz),k);
      s=gsub(gun,gmul(ak,s));
    }
    gmulz(s,zf,y);
  }
  gdivz(y,mpexp(x),y); avma=av; return y;
}

/***********************************************************************/
/*                                                                    **/
/**                    FONCTION U(a,b,z) GENERALE                     **/
/**                         ET CAS PARTICULIERS                       **/
/**                                                                   **/
/***********************************************************************/

/* Assume gx > 0 and a,b real */
GEN
hyperu(GEN a, GEN b, GEN gx, long prec)
{
  GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn;
  long l,lbin,av,av1,av2,k,l1,n,ex;

  if (typ(gx)!=t_REAL)
    { l=prec; k=1; }
  else
    { l=lg(gx); k=0; x=gx; }
  ex = (iscomplex(a) || iscomplex(b));
  if (ex) { y=cgetg(3,t_COMPLEX); y[1]=lgetr(l); y[2]=lgetr(l); }
  else y=cgetr(l);
  l1=l+1; av=avma; if (k) gaffect(gx,x=cgetr(l));
  a1=gaddsg(1,gsub(a,b));
  if (ex)
  {
    u=cgetg(3,t_COMPLEX); u[1]=lgetr(l1); u[2]=lgetr(l1);
    v=cgetg(3,t_COMPLEX); v[1]=lgetr(l1); v[2]=lgetr(l1);
    c=cgetg(3,t_COMPLEX); c[1]=lgetr(l1); c[2]=lgetr(l1);
    d=cgetg(3,t_COMPLEX); d[1]=lgetr(l1); d[2]=lgetr(l1);
    e=cgetg(3,t_COMPLEX); e[1]=lgetr(l1); e[2]=lgetr(l1);
    f=cgetg(3,t_COMPLEX); f[1]=lgetr(l1); f[2]=lgetr(l1);
  }
  else
  {
    u=cgetr(l1); v=cgetr(l1); c=cgetr(l1);
    d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
  }
  n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
  lbin = 10-bit_accuracy(l); av1=avma;
  if (gcmpgs(x,n)<0)
  {
    gn=stoi(n); zf=gpui(gn,gneg_i(a),l1);
    zz=gdivsg(-1,gn); s=gun; t=gzero;
    for (k=n-1; k>=0; k--)
    {
      p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
      s=gaddsg(1,gmul(p1,s));
      t=gadd(gaddsg(k,a),gmul(p1,t));
    }
    gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1;
    q=cgetr(l1); affsr(n,q); r=x; av1=avma;
    do
    {
      p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf;
      p2=subsr(1,divrr(r,q)); if (gcmp(p1,p2)>0) p1=p2;
      gnegz(p1,c); avma=av1;
      k=0; gaffsg(1,d);
      gaffect(u,e); gaffect(v,f);
      p3=gsub(q,b); av2=avma;
      for(;;)
      {
	w=gadd(gmul(gaddsg(k,a),u),gmul(gaddsg(-k,p3),v));
	k++; gdivgsz(gmul(q,v),k,u);
	gdivgsz(w,k,v);
	gmulz(d,c,d);
	gaddz(e,gmul(d,u),e); p1=gmul(d,v);
	gaddz(f,p1,f);
	if (gexpo(p1)-gexpo(f) <= lbin) break;
	avma=av2;
      }
      p1=u; u=e; e=p1;
      p1=v; v=f; f=p1;
      gmulz(q,gadd(gun,c),q);
      p1=subrr(q,r); ex=expo(p1); avma=av1;
    }
    while (ex>lbin);
  }
  else
  {
    zf=gpui(x,gneg_i(a),l1);
    zz=gdivsg(-1,x); s=gun;
    for (k=n-1; k>=0; k--)
    {
      p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
      s=gadd(gun,gmul(p1,s));
    }
    u = gmul(s,zf);
  }
  gaffect(u,y); avma=av; return y;
}

GEN
kbessel2(GEN nu, GEN x, long prec)
{
  long av = avma,tetpil;
  GEN p1,p2,x2,a,pitemp;

  if (typ(x)==t_REAL) prec = lg(x);
  x2=gshift(x,1); pitemp=mppi(prec);
  if (gcmp0(gimag(nu))) a=cgetr(prec);
  else
  { a=cgetg(3,t_COMPLEX); a[1]=lgetr(prec); a[2]=lgetr(prec); }
  gaddz(gun,gshift(nu,1),a);
  p1=hyperu(gshift(a,-1),a,x2,prec);
  p1=gmul(gmul(p1,gpui(x2,nu,prec)),mpsqrt(pitemp));
  p2=gexp(x,prec); tetpil=avma;
  return gerepile(av,tetpil,gdiv(p1,p2));
}

GEN
incgam(GEN a, GEN x, long prec)
{
  GEN p1,z = cgetr(prec);
  long av = avma;

  if (typ(x)!=t_REAL) { gaffect(x,z); x=z; }
  if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
    p1 = incgam2(a,x,prec);
  else
    p1 = gsub(ggamma(a,prec), incgam3(a,x,prec));
  affrr(p1,z); avma = av; return z;
}

/* = incgam2(0, x, prec). typ(x) = t_REAL. Optimized for eint1 */
static GEN
incgam2_0(GEN x)
{
  long l,n,i;
  GEN p1;
  double m,mx;

  l = lg(x); mx = rtodbl(x);
  m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
  p1 = divsr(-n, addsr(n<<1,x));
  for (i=n-1; i >= 1; i--)
    p1 = divsr(-i, addrr(addsr(i<<1,x), mulsr(i,p1)));
  return mulrr(divrr(mpexp(negr(x)), x), addrr(realun(l),p1));
}

GEN
incgam1(GEN a, GEN x, long prec)
{
  GEN p2,p3,y, z = cgetr(prec);
  long av=avma, av1, l,n,i;
  double m,mx;

  if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
  l=lg(x); mx=rtodbl(x);
  m=(long) bit_accuracy(l)*LOG2; n=(long)(m/(log(m)-(1+log(mx))));
  p2 = cgetr(l); affrr(addir(gun,gsub(x,a)), p2);
  p3 = subrs(p2, n+1); av1 = avma;
  for (i=n; i>=1; i--)
  {
    affrr(addrr(subrs(p2,i), divrr(mulsr(i,x),p3)), p3);
    avma = av1;
  }
  y = mulrr(mpexp(negr(x)),gpui(x,a,prec));
  affrr(divrr(y,p3), z);
  avma = av; return z;
}

GEN
incgam2(GEN a, GEN x, long prec)
{
  GEN b,p1,p2,p3,y, z = cgetr(prec);
  long av = avma,av1,l,n,i;
  double m,mx;

  if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
  if (gcmp0(a)) { affrr(incgam2_0(x), z); avma = av; return z; }
  l=lg(x); mx=rtodbl(x);
  m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
  i = typ(a);
  if (i == t_REAL) b = addsr(-1,a);
  else
  { /* keep b integral : final powering more efficient */
    gaffect(a,p1=cgetr(prec));
    b = (i == t_INT)? addsi(-1,a): addsr(-1,p1);
    a = p1;
  }
  p2 = cgetr(l); gaffect(subrr(x,a),p2);
  p3 = divrr(addsr(-n,a), addrs(p2,n<<1));
  av1 = avma;
  for (i=n-1; i>=1; i--)
  {
    affrr(divrr(addsr(-i,a), addrr(addrs(p2,i<<1),mulsr(i,p3))), p3);
    avma = av1;
  }
  y = gmul(mpexp(negr(x)), gpui(x,b,prec));
  affrr(mulrr(y, addsr(1,p3)), z);
  avma = av; return z;
}

GEN
incgam3(GEN a, GEN x, long prec)
{
  GEN b,p1,p2,p3,y, z = cgetr(prec);
  long av = avma, av1,l,n,i;

  if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
  l=lg(x); n = -bit_accuracy(l)-1;
  p3 = realun(l);
  p2 = rcopy(p3);
  i = typ(a);
  if (i == t_REAL) b = a;
  else
  {
    gaffect(a,p1=cgetr(prec));
    b = (i == t_INT)? a: p1;
    a = p1;
  }
  av1 = avma;
  for (i=1; expo(p2) >= n; i++)
  {
    affrr(divrr(mulrr(x,p2), addsr(i,a)), p2);
    affrr(addrr(p2,p3), p3); avma = av1;
  }
  y = gdiv(mulrr(mpexp(negr(x)), gpui(x,b,prec)), a);
  affrr(mulrr(y,p3), z);
  avma = av; return z;
}

/* Assumes that g=gamma(a,prec). Unchecked */
GEN
incgam4(GEN a, GEN x, GEN g, long prec)
{
  GEN p1, z = cgetr(prec);
  long av = avma;

  if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
  if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
    p1 = incgam2(a,x,prec);
  else
    p1 = gsub(g, incgam3(a,x,prec));
  affrr(p1, z);
  avma = av; return z;
}

GEN
incgam0(GEN a, GEN x, GEN z,long prec)
{
  return z? incgam4(a,x,z,prec): incgam(a,x,prec);
}

GEN
eint1(GEN x, long prec)
{
  long av = avma,tetpil,l,n,i;
  GEN p1,p2,p3,p4,run,y;

  if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;}
  if (signe(x) >= 0)
  {
    if (expo(x) >= 4)
      return gerepileupto(av, incgam2_0(x));

    l = lg(x); consteuler(l);
    n = -bit_accuracy(l)-1;
    
    run = realun(l);
    p4 = p3 = p2 = p1 = run;
    for (i=2; expo(p2)>=n; i++)
    {
      p1 = addrr(p1,divrs(run,i)); /* p1 = sum_{i=1} 1/i */
      p4 = divrs(mulrr(x,p4),i);   /* p4 = sum_{i=1} x^(i-1)/i */
      p2 = mulrr(p4,p1);
      p3 = addrr(p2,p3);
    }
    p3 = mulrr(x,mulrr(mpexp(negr(x)),p3));
    p1 = addrr(mplog(x),geuler);
    return gerepileupto(av, subrr(p3,p1));
  }
  else
  { /* written by Manfred Radimersky */
    l  = lg(x);
    n  = bit_accuracy(l);
    /* IS: line split to avoid a Workshop cc-5.0 bug (Sun BugID #4254959) */
    n  = 3 * n / 4;
    y  = negr(x);
    if(gcmpgs(y, n) < 0) {
      p1 = p2 = p3 = y;
      p4 = gzero;
      i  = 2;
      while(gcmp(p3, p4)) {
        p4 = p3;
        p1 = gmul(p1, gdivgs(y, i));
        p2 = gdivgs(p1, i);
        p3 = gadd(p3, p2);
        i++;
      }
      consteuler(l);
      p1 = gadd(mplog(y), geuler);
      y  = gadd(p3, p1);
    } else {
      p1 = gdivsg(1, y);
      p2 = realun(l);
      p3 = p2;
      p4 = gzero;
      i  = 1;
      while(gcmp(p3, p4)) {
        p4 = p3;
        p2 = gmulgs(p2, i);
        p2 = gmul(p2, p1);
        p3 = gadd(p3, p2);
        i++;
      }
      p1 = gdiv(mpexp(y), y);
      y  = gmul(p3, p1);
    }
    tetpil = avma;
    y  = gerepile(av, tetpil, negr(y));
  }
  return y;
}

GEN
veceint1(GEN C, GEN nmax, long prec)
{
  long av,av1,k,n,nstop,i,cd,nmin,G,a;
  GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit;

  if (!nmax) return eint1(C,prec);

  if (signe(nmax)<=0) return cgetg(1,t_VEC);
  if (DEBUGLEVEL>1) fprintferr("Entering veceint1:\n");
  if (typ(C) != t_REAL || lg(C) > prec)
    { y = cgetr(prec); gaffect(C,y); C = y; }
  if (signe(C) <= 0) err(talker,"negative or zero constant in veceint1");

  G = -bit_accuracy(prec);
  n=itos(nmax); y=cgetg(n+1,t_VEC);
  for(i=1; i<=n; i++) y[i]=lgetr(prec);
  av=avma;

  nstop = itos(gceil(divsr(4,C)));
  if (nstop<1) nstop=1;
  if (nstop>n) nstop=n;
  nmin=n-10; if (nmin<nstop) nmin=nstop;
  if(DEBUGLEVEL>1) fprintferr("nstop = %ld\n",nstop);

  e1=mpexp(negr(mulsr(n,C)));
  e2=mpexp(mulsr(10,C));
  unr = realun(prec);
  zeror=realzero(prec); deninit=negr(unr);
  f=cgetg(3,t_COL); M2=cgetg(3,t_VEC); av1=avma;

  F0=(GEN)y[n];
  affrr(eint1(mulsr(n,C),prec), F0);
  do
  {
    if (DEBUGLEVEL>1) fprintferr("%ld ",n);
    minvn=divrs(unr,-n); mcn=mulrr(C,minvn);
    M2[1] = (long)zeror; M2[2] = lsubrr(minvn,C);
    f[1]=(long)zeror; f[2]=ldivrs(e1,-n);
    affrr(mulrr(e1,e2), e1);
    vdiff=cgetg(2,t_VEC); vdiff[1]=f[2];
    for (cd=a=1,n--; n>=nmin; n--,a++)
    {
      F = F0;
      ap = stoi(a); den = deninit;
      for (k=1;;)
      {
	if (k>cd)
        {
          cd++; p1 = (GEN)f[2];
          f[2] = lmul(M2,f);
          f[1] = (long)p1;
          M2[1] = laddrr((GEN)M2[1],mcn);
          M2[2] = laddrr((GEN)M2[2],minvn);
          vdiff = concatsp(vdiff,(GEN)f[2]);
        }
	p1 = mulrr(mulri(den,ap),(GEN)vdiff[k]);
        if (expo(p1) < G) { affrr(F,(GEN)y[n]); break; }
	F = addrr(F,p1); ap = mulis(ap,a);
	k++; den = divrs(den,-k);
      }
    }
    avma=av1; n++; F0=(GEN)y[n]; nmin -= 10;
    if (nmin < nstop) nmin=nstop;
  }
  while(n>nstop);
  for(i=1; i<=nstop; i++)
    affrr(eint1(mulsr(i,C),prec), (GEN)y[i]);
  if (DEBUGLEVEL>1) fprintferr("\n");
  avma=av; return y;
}

GEN
gerfc(GEN x, long prec)
{
  long av=avma;
  GEN p1,p2;

  if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; }
  av = avma; p1 = incgam(ghalf,gsqr(x),prec);
  p2 = mpsqrt(mppi(lg(x)));
  p1 = divrr(p1,p2);
  if (signe(x) < 0) p1 = subsr(2,p1);
  return gerepileupto(av,p1);
}

/***********************************************************************/
/**								      **/
/**		        FONCTION ZETA DE RIEMANN                      **/
/**								      **/
/***********************************************************************/

static GEN
czeta(GEN s, long prec)
{
  long av,n,p,n1,l,flag1,flag2,flag3,i,i2;
  double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf;
  GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp;

  l=precision(s);
  if (typ(s)==t_COMPLEX)
  {
    if (!l) l=prec;
    res=cgetg(3,t_COMPLEX); res[1]=lgetr(l); res[2]=lgetr(l);
    av=avma;
    p1=cgetg(3,t_COMPLEX); p1[1]=lgetr(l+1); p1[2]=lgetr(l+1);
    gaffect(s,p1); s=p1; sig=(GEN)s[1];
  }
  else
  {
    res = cgetr(l); av=avma;
    p1=cgetr(l+1); affrr(s,p1); sig=s=p1;
  }

  if (signe(sig)>0 && expo(sig)>-2) flag1 = 0;
  else
  {
    if (gcmp0(gimag(s)) && gcmp0(gfrac(gmul2n(sig,-1))))
    {
      if (gcmp0(sig)) gaffect(gneg_i(ghalf),res); else gaffsg(0,res);
      avma=av; return res;
    }
    flag1=1; s=gsub(gun,s); sig=greal(s);
  }
  ssig=rtodbl(sig); st=fabs(rtodbl(gimag(s))); maxbeta = pow(3.0,-2.5);
  if (st)
  {
    ns = ssig*ssig + st*st;
    alpha=pariC2*(prec-2)-0.39-0.5*(ssig-1.0)*log(ns)-log(ssig)+ssig*2*pariC1;
    beta=(alpha+ssig)/st-atan(ssig/st);
    if (beta<=0)
    {
      if (ssig>=1.0)
      {
	p=0; sn=sqrt(ns);
	n=(long)(ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig)));
      }
      else
      {
	p=1; sn=ssig+1; n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
      }
    }
    else
    {
      if (beta<maxbeta) xinf=beta+pow(3*beta,1.0/3.0);
      else
      {
	double eps=0.0087, x00 = beta+PI/2.0, y00,x11;
        for(;;)
	{
	  y00=x00*x00; x11=(beta+atan(x00))*(1+y00)/y00-1/x00;
	  if (x00-x11 < eps) break;
	  x00 = x11;
	}
	xinf=x11;
      }
      sp=1.0-ssig+st*xinf;
      if (sp>0)
      {
	p=(long)ceil(sp/2.0); sn=ssig+2*p-1;
	n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
      }
      else
      {
	p=0; sn=sqrt(ns);
	n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
      }
    }
  }
  else
  {
    beta=pariC2*(prec-2)+0.61+ssig*2*pariC1-ssig*log(ssig);
    if (beta>0)
    {
      p=(long)ceil(beta/2.0); sn=ssig+2*p-1;
      n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
    }
    else
    {
      p=0; sn=sqrt(ssig*ssig+st*st);
      n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
    }
  }
  if (n < 46340) { flag2=1; n1=n*n; } else flag2=0;
  y=gun; ms=gneg_i(s); p1=cgetr(prec+1);
  for (i=2; i<=n; i++)
  {
    affsr(i,p1); p2 = gexp(gmul(ms,mplog(p1)), prec+1);
    y = gadd(y,p2);
  }
  flag3 = (2*p < 46340);
  mpbern(p,prec+1); p31=cgetr(prec+1); z=gzero;
  for (i=p; i>=1; i--)
  {
    i2=i<<1;
    p1=gmul(gaddsg(i2-1,s),gaddsg(i2,s));
    p1=flag3? gdivgs(p1,i2*(i2+1)): gdivgs(gdivgs(p1,i2),i2+1);
    p1=flag2? gdivgs(p1,n1): gdivgs(gdivgs(p1,n),n);
    p3 = bern(i);
    if (bernzone[2]>prec+1) { affrr(p3,p31); p3=p31; }
    z=gadd(divrs(p3,i),gmul(p1,z));
  }
  p1=gsub(gdivsg(n,gsubgs(s,1)),ghalf);
  z=gmul(gadd(p1,gmul(s,gdivgs(z,n<<1))),p2);
  y = gadd(y,z);
  if (flag1)
  {
    pitemp=mppi(prec+1); setexpo(pitemp,2);
    y=gmul(gmul(y,ggamma(s,prec+1)),gpui(pitemp,ms,prec+1));
    setexpo(pitemp,0);
    y=gmul2n(gmul(y,gcos(gmul(pitemp,s),prec+1)),1);
  }
  gaffect(y,res); avma=av; return res;
}

/* y = binomial(n,k-2). Return binomial(n,k) */
static GEN
next_bin(GEN y, long n, long k)
{
  y = divrs(mulrs(y, n-k+2), k-1);
  return divrs(mulrs(y, n-k+1), k);
}

static GEN
izeta(long k, long prec)
{
  long av=avma,av2,tetpil,kk,n,li, limit;
  GEN y,p1,pitemp,qn,z,q,binom;

  if (!k) return gneg(ghalf);
  if (k<0)
  {
    if ((k&1) == 0) return gzero;
    y=bernreal(1-k,prec); tetpil=avma;
    return gerepile(av,tetpil,divrs(y,k-1));
  }
  if (k > bit_accuracy(prec)+1) return realun(prec);
  pitemp=mppi(prec); setexpo(pitemp,2);
  if ((k&1) == 0)
  {
    p1 = mulrr(gpuigs(pitemp,k),absr(bernreal(k,prec)));
    y = mpfactr(k,prec); tetpil=avma;
    y=divrr(p1,y); setexpo(y,expo(y)-1);
    return gerepile(av,tetpil,y);
  }
  binom = realun(prec+1);
  q = mpexp(pitemp); kk = k+1;
  li = -(1+bit_accuracy(prec));
  if ((k&3)==3)
  {
    for (n=0; n <= kk>>1; n+=2)
    {
      p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
      if (n) { binom = next_bin(binom,kk,n); setlg(binom,prec+1); }
      p1 = mulrr(binom,p1);
      if (n == kk>>1) setexpo(p1, expo(p1)-1);
      if ((n>>1)&1) setsigne(p1,-signe(p1));
      y = n? addrr(y,p1): p1;
    }
    y=mulrr(divrr(gpuigs(pitemp,k),mpfactr(kk,prec)),y);

    av2 = avma; limit = stack_lim(av2,1);
    qn=gsqr(q); z=divsr(1,addrs(q,-1));
    for (n=2; ; n++)
    {
      p1=divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1)));

      z=addrr(z,p1); if (expo(p1)< li) break;
      qn=mulrr(qn,q);
      if (low_stack(limit,stack_lim(av2,1)))
      {
        GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
        if (DEBUGMEM>1) err(warnmem,"izeta");
        gerepilemany(av2,gptr,2);
      }
    }
    setexpo(z,expo(z)+1); tetpil = avma;
    y = addrr(y,z); setsigne(y,-signe(y));
  }
  else
  {
    GEN p2 = divrs(pitemp, k-1);
    for (n=0; n <= k>>1; n+=2)
    {
      p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
      if (n) binom = next_bin(binom,kk,n);
      p1 = mulrr(binom,p1);
      p1 = mulsr(kk-(n<<1),p1);
      if ((n>>1)&1) setsigne(p1,-signe(p1));
      y = n? addrr(y,p1): p1;
    }
    y = mulrr(divrr(gpuigs(pitemp,k),mpfactr(kk,prec)),y);
    y = divrs(y,k-1);
    av2 = avma; limit = stack_lim(av2,1);
    qn = q; z=gzero;
    for (n=1; ; n++)
    {
      p1=mulir(gpuigs(stoi(n),k),gsqr(addrs(qn,-1)));
      p1=divrr(addrs(mulrr(qn,addsr(1,mulsr(n<<1,p2))),-1),p1);

      z=addrr(z,p1); if (expo(p1) < li) break;
      qn=mulrr(qn,q);
      if (low_stack(limit,stack_lim(av2,1)))
      {
        GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
        if (DEBUGMEM>1) err(warnmem,"izeta");
        gerepilemany(av2,gptr,2);
      }
    }
    setexpo(z,expo(z)+1); tetpil = avma;
    y = subrr(y,z);
  }
  return gerepile(av,tetpil,y);
}

GEN
gzeta(GEN x, long prec)
{
  if (gcmp1(x)) err(talker,"argument equal to one in zeta");
  switch(typ(x))
  {
    case t_INT:
      return izeta(itos(x),prec);

    case t_REAL: case t_COMPLEX:
      return czeta(x,prec);

    case t_INTMOD: case t_PADIC:
      err(typeer,"gzeta");
    case t_SER:
      err(impl,"zeta of power series");
  }
  return transc(gzeta,x,prec);
}

void
gzetaz(GEN x, GEN y)
{
  long av=avma, prec = precision(y);

  if (!prec) err(infprecer,"gzetaz");
  gaffect(gzeta(x,prec),y); avma=av;
}

/***********************************************************************/
/**                                                                   **/
/**                    FONCTIONS POLYLOGARITHME                       **/
/**                                                                   **/
/***********************************************************************/

/* validity domain contains .005 < |x| < 230 */
static GEN
cxpolylog(long m, GEN x, long prec)
{
  long av=avma,li,i,n,bern_upto;
  GEN p1,z,h,q,s;

  if (gcmp1(x)) return izeta(m,prec);
  z=glog(x,prec); h=gneg_i(glog(gneg_i(z),prec));
  for (i=1; i<m; i++) h = gadd(h,ginv(stoi(i)));

  bern_upto=m+50; mpbern(bern_upto,prec);
  q=gun; s=izeta(m,prec);
  for (n=1; n<=m+1; n++)
  {
    q=gdivgs(gmul(q,z),n);
    s=gadd(s,gmul((n==(m-1))? h: izeta(m-n,prec),q));
  }

  n = m+3; z=gsqr(z); li = -(bit_accuracy(prec)+1);

  for(;;)
  {
    q = gdivgs(gmul(q,z),(n-1)*n);
    p1 = gmul(izeta(m-n,prec),q);
    s = gadd(s,p1);
    if (gexpo(p1) < li) break;
    n+=2;
    if (n>=bern_upto) { bern_upto += 50; mpbern(bern_upto,prec); }
  }
  return gerepileupto(av,s);
}

GEN
polylog(long m, GEN x, long prec)
{
  long av,av1,limpile,tetpil,l,e,sx,i;
  GEN p1,p2,n,y,logx,logx2;
  GEN *gptr[4];

  if (m<0) err(talker,"negative index in polylog");
  if (!m) return gneg(ghalf);
  if (gcmp0(x)) return gcopy(x);
  av=avma;
  if (m==1)
  {
    p1=glog(gsub(gun,x),prec); tetpil=avma;
    return gerepile(av,tetpil,gneg(p1));
  }
  l=precision(x);
  if (!l) { l=prec; x=gmul(x, realun(l)); }
  e=gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
  if (e>0)
  {
    logx=glog(x,l); sx=gsigne(gimag(x));
    if (!sx)
    {
      if (m&1)
        sx = gsigne(greal(gsub(gun,x)));
      else
        sx = -gsigne(greal(x));
    }
    x=ginv(x);
  }
  av1=avma; limpile=stack_lim(av1,1);
  y=x; n=gun; p1=x;
  do
  {
    n=addsi(1,n); p1=gmul(x,p1); p2=gdiv(p1,gpuigs(n,m));
    y=gadd(y,p2);
    if (low_stack(limpile, stack_lim(av1,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"polylog");
      gptr[0]=&y; gptr[1]=&n; gptr[2]=&p1; gptr[3]=&p2;
      gerepilemany(av1,gptr,4);
    }
  }
  while (gexpo(p2) > -bit_accuracy(l));
  tetpil=avma;
  if (e<=0) return gerepile(av,tetpil,gcopy(y));
  logx2=gsqr(logx); p1=gneg_i(ghalf);
  for (i=m-2; i>=0; i-=2)
  {
    p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
  }
  if (m&1) p1 = gmul(logx,p1);
  p2=cgetg(3,t_COMPLEX); p2[1]=zero;
  p2[2]=ldiv(gmulsg(sx,mppi(l)),mpfact(m-1));
  p1=gadd(gmul2n(p1,1), gmul(p2,gpuigs(logx,m-1)));
  if ((m&1) == 0) y = gneg_i(y);
  tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));
}

GEN
dilog(GEN x, long prec)
{
  GEN p1,p2,p3,y;
  long av=avma,tetpil;

  if (gcmpgs(gnorm(x),1)<1)
  {
    tetpil=avma; return gerepile(av,tetpil,polylog(2,x,prec));
  }
  y=gneg_i(polylog(2,ginv(x),prec)); p3=mppi(prec); p2=gdivgs(gsqr(p3),6);
  p1=glog(gneg_i(x),prec); p1=gadd(p2,gmul2n(gsqr(p1),-1));
  p1 = gneg_i(p1); tetpil=avma;
  return gerepile(av,tetpil,gadd(y,p1));
}

GEN
polylogd0(long m, GEN x, long flag, long prec)
{
  long k,l,av,tetpil,fl,m2;
  GEN p1,p2,p3,y;

  m2=m&1; av=avma;
  if (gcmp0(x)) return gcopy(x);
  if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
  l=precision(x);
  if (!l) { l=prec; x=gmul(x,realun(l)); }
  p1=gabs(x,prec); fl=0;
  if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }

  p1=gneg_i(glog(p1,prec)); p2=gun;
  y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
  for (k=1; k<m; k++)
  {
    p2=gdivgs(gmul(p2,p1),k);
    p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
    tetpil=avma; y=gadd(y,gmul(p2,p3));
  }
  if (m2)
  {
    if (flag)
      p2 = gdivgs(gmul(p2,p1),-2*m);
    else
      p2 = gdivgs(gmul(glog(gnorm(gsub(gun,x)),prec),p2),2*m);
    tetpil=avma; y=gadd(y,p2);
  }
  if (fl) { tetpil=avma; return gerepile(av,tetpil,gneg(y)); }
  return gerepile(av,tetpil,y);
}

GEN
polylogd(long m, GEN x, long prec)
{
  return polylogd0(m,x,0,prec);
}

GEN
polylogdold(long m, GEN x, long prec)
{
  return polylogd0(m,x,1,prec);
}

GEN
polylogp(long m, GEN x, long prec)
{
  long k,l,av,tetpil,fl,m2;
  GEN p1,y;

  m2=m&1; av=avma;
  if (gcmp0(x)) return gcopy(x);
  if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
  l=precision(x);
  if (!l) { l=prec; x=gmul(x,realun(l)); }
  p1=gabs(x,prec); fl=0;
  if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }

  p1=gmul2n(glog(p1,prec),1); mpbern(m>>1,prec);
  y=polylog(m,x,prec); y=m2?greal(y):gimag(y);

  if (m==1)
  {
    p1=gmul2n(p1,-2); tetpil=avma; y=gadd(y,p1);
  }
  else
  {
    GEN p2=gun, p3, p4, p5, p51=cgetr(prec);

    for (k=1; k<m; k++)
    {
      p2=gdivgs(gmul(p2,p1),k);
      if (!(k&1) || k==1)
      {
	if (k!=1)
	{
	  p5=(GEN)bern(k>>1);
	  if (bernzone[2]>prec) { affrr(p5,p51); p5=p51; }
	  p4=gmul(p2,p5);
	}
	else p4=gneg_i(gmul2n(p2,-1));
	p3=polylog(m-k,x,prec); p3=m2?greal(p3):gimag(p3);
	tetpil=avma; y=gadd(y,gmul(p4,p3));
      }
    }
  }
  if (fl) { tetpil=avma; return gerepile(av,tetpil,gneg(y)); }
  return gerepile(av,tetpil,y);
}

GEN
gpolylog(long m, GEN x, long prec)
{
  long i,lx,av=avma,tetpil,v,n;
  GEN y,p1,p2;

  if (m<=0)
  {
    p1=polx[0]; p2=gsub(gun,p1);
    for (i=1; i<=(-m); i++)
      p1=gmul(polx[0],gadd(gmul(p2,derivpol(p1)),gmulsg(i,p1)));
    p1=gdiv(p1,gpuigs(p2,1-m)); tetpil=avma;
    return gerepile(av,tetpil,poleval(p1,x));
  }

  switch(typ(x))
  {
    case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
    case t_COMPLEX: case t_QUAD:
      return polylog(m,x,prec);

    case t_INTMOD: case t_PADIC: case t_POLMOD:
      p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
      for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
      tetpil=avma; y=cgetg(lx,t_COL);
      for (i=1; i<lx; i++) y[i]=(long)polylog(m,(GEN)p2[i],prec);
      return gerepile(av,tetpil,y);

    case t_POL: case t_RFRAC: case t_RFRACN:
      p1=tayl(x,gvar(x),precdl); tetpil=avma;
      return gerepile(av,tetpil,gpolylog(m,p1,prec));

    case t_SER:
      if (!m) return gneg(ghalf);
      if (m==1)
      {
	p1=glog(gsub(gun,x),prec); tetpil=avma;
	return gerepile(av,tetpil,gneg(p1));
      }
      if (valp(x)<=0) err(impl,"polylog around a!=0");
      v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2);
      for (i=n; i>=1; i--)
      {
	p1=gadd(gpuigs(stoi(i),-m),y); tetpil=avma;
	y=gmul(x,p1);
      }
      return gerepile(av,tetpil,y);

    case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,typ(x));
      for (i=1; i<lx; i++)
	y[i]=(long)gpolylog(m,(GEN)x[i],prec);
      return y;
  }
  err(typeer,"gpolylog");
  return NULL; /* not reached */
}

void
gpolylogz(long m, GEN x, GEN y)
{
  long av=avma, prec = precision(y);

  if (!prec) err(infprecer,"gpolylogz");
  gaffect(gpolylog(m,x,prec),y); avma=av;
}

GEN
polylog0(long m, GEN x, long flag, long prec)
{
  switch(flag)
  {
    case 0: return gpolylog(m,x,prec);
    case 1: return polylogd(m,x,prec);
    case 2: return polylogdold(m,x,prec);
    case 3: return polylogp(m,x,prec);
    default: err(flagerr,"polylog");
  }
  return NULL; /* not reached */
}

static GEN
qq(GEN x, long prec)
{
  long tx=typ(x);

  if (tx==t_PADIC) return x;
  if (is_scalar_t(tx))
  {
    long l=precision(x);
    GEN p1,p2,q;

    if (tx != t_COMPLEX || gcmp((GEN)x[2],gzero)<=0)
      err(talker,"argument must belong to upper half-plane");

    if (!l) l=prec; p1=mppi(l); setexpo(p1,2);
    p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
    q=gmul(x,p2); return gexp(q,prec);
  }
  if (tx != t_POL && tx != t_SER)
    err(talker,"bad argument for modular function");

  return tayl(x,gvar(x),precdl);
}

static GEN
inteta(GEN q)
{
  long tx=typ(q);
  GEN p1,ps,qn,y;

  y=gun; qn=gun; ps=gun;
  if (tx==t_PADIC)
  {
    for(;;)
    {
      p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
      y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
      p1=y; y=gadd(y,ps); if (gegal(p1,y)) return y;
    }
  }
  else
  {
    long l,v, av = avma, lim = stack_lim(av,3);

    if (is_scalar_t(tx)) l = -bit_accuracy(precision(q));
    else
      { v=gvar(q); tx = 0; }
    for(;;)
    {
      p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
      y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
      y=gadd(y,ps);
      if (tx)
        { if (gexpo(ps)-gexpo(y) < l) return y; }
      else
        { if (gval(ps,v)>=precdl) return y; }
      if (low_stack(lim, stack_lim(av,3)))
      {
        GEN *gptr[3];
        if(DEBUGMEM>1) err(warnmem,"inteta");
        gptr[0]=&y; gptr[1]=&qn; gptr[2]=&ps;
        gerepilemany(av,gptr,3);
      }
    }
  }
}

GEN
eta(GEN x, long prec)
{
  long av = avma;
  GEN q = qq(x,prec);
  return gerepileupto(av,inteta(q));
}

/* returns the true value of eta(x) for Im(x) > 0, using reduction */
GEN
trueeta(GEN x, long prec)
{
  long tx=typ(x), av=avma, tetpil,l;
  GEN p1,p2,q,q24,n,z,m,unapprox;

  if (!is_scalar_t(tx)) err(typeer,"trueeta");
  if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0)
    err(talker,"argument must belong to upper half-plane");
  l=precision(x); if (l) prec=l;
  p1=mppi(prec); setexpo(p1,2);
  p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
  z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */
  unapprox=gsub(gun,gpuigs(stoi(10),-8));
  m=gun;
  for(;;)
  {
    n=ground(greal(x));
    if (signe(n)) {x=gsub(x,n); m=gmul(m,powgi(z,n));}
    if (gcmp(gnorm(x), unapprox)>=0) break;
    m=gmul(m,gsqrt(gdiv(gi,x),prec)); x=gdivsg(-1,x);
  }
  q=gmul(p2,x);
  q24=gexp(gdivgs(q,24),prec); q=gpuigs(q24,24);
  p1=gmul(m,q24); q = inteta(q); tetpil=avma;
  return gerepile(av,tetpil,gmul(p1,q));
}

GEN
eta0(GEN x, long flag,long prec)
{
  return flag? trueeta(x,prec): eta(x,prec);
}

GEN
jell(GEN x, long prec)
{
  long av=avma,tetpil,tx = typ(x);
  GEN p1;

  if (!is_scalar_t(tx))
  {
    GEN p2,q = qq(x,prec);
    p1=gdiv(inteta(gsqr(q)), inteta(q));
    p1=gmul2n(gsqr(p1),1);
    p1=gmul(q,gpuigs(p1,12));
    p2=gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
    p1=gmulsg(48,p1); tetpil=avma;
    return gerepile(av,tetpil,gadd(p2,p1));
  }
  p1=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
  p1=gsqr(gsqr(gsqr(p1)));
  p1=gadd(gmul2n(gsqr(p1),8), ginv(p1)); tetpil=avma;
  return gerepile(av,tetpil,gpuigs(p1,3));
}

GEN
wf2(GEN x, long prec)
{
  long av=avma,tetpil;
  GEN p1,p2;

  p1=gsqrt(gdeux,prec);
  p2=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
  tetpil=avma;
  return gerepile(av,tetpil,gmul(p1,p2));
}

GEN
wf1(GEN x, long prec)
{
  long av=avma,tetpil;
  GEN p1,p2;

  p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec);
  tetpil=avma;
  return gerepile(av,tetpil,gdiv(p1,p2));
}

GEN
wf(GEN x, long prec)
{
  long av=avma,tetpil;
  GEN p1,p2;

  p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec));
  p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=ldivrs(mppi(prec),-24);
  p2=gexp(p2,prec); tetpil=avma;
  return gerepile(av,tetpil,gmul(p1,p2));
}

GEN
weber0(GEN x, long flag,long prec)
{
  switch(flag)
  {
    case 0: return wf(x,prec);
    case 1: return wf1(x,prec);
    case 2: return wf2(x,prec);
    default: err(flagerr,"weber");
  }
  return NULL; /* not reached */
}

static GEN
sagm(GEN x, long prec)
{
  GEN p1,a,b,a1,b1,y;
  long av,tetpil,l,ep;

  if (gcmp0(x)) return gcopy(x);
  switch(typ(x))
  {
    case t_REAL:
      l = precision(x); y = cgetr(l); av=avma;
      a1 = x; b1 = realun(l);
      l = 5-bit_accuracy(prec);
      do
      {
	a=a1; b=b1; a1 = addrr(a,b);
        setexpo(a1,expo(a1)-1);
	b1=mpsqrt(mulrr(a,b));
      }
      while (expo(subrr(b1,a1))-expo(b1) >= l);
      affrr(a1,y); avma=av; return y;

    case t_COMPLEX:
      if (gcmp0((GEN)x[2]))
        return transc(sagm,(GEN)x[1],prec);
      av=avma; l=precision(x); if (l) prec=l;
      a1=x; b1=gun; l = 5-bit_accuracy(prec);
      do
      {
	a=a1; b=b1;
	a1=gmul2n(gadd(a,b),-1);
	b1=gsqrt(gmul(a,b),prec);
      }
      while (gexpo(gsub(b1,a1))-gexpo(b1) >= l);
      tetpil=avma; return gerepile(av,tetpil,gcopy(a1));

    case t_PADIC:
      av=avma; a1=x; b1=gun; l=precp(x);
      do
      {
	a=a1; b=b1;
	a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
	p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
	if (ep<=0) { b1=gneg_i(b1); p1=gsub(b1,a1); ep=valp(p1)-valp(b1); }
      }
      while (ep<l && !gcmp0(p1));
      tetpil=avma; return gerepile(av,tetpil,gcopy(a1));

    case t_SER:
      av=avma; a1=x; b1=gun; l=lg(x)-2;
      do
      {
	a=a1; b=b1;
	a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
	p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
      }
      while (ep<l && !gcmp0(p1));
      tetpil=avma; return gerepile(av,tetpil,gcopy(a1));

    case t_INTMOD:
      err(impl,"agm of mod");
  }
  return transc(sagm,x,prec);
}

GEN
agm(GEN x, GEN y, long prec)
{
  long av,tetpil, ty=typ(y);
  GEN z;

  if (is_matvec_t(ty))
  {
    ty=typ(x);
    if (is_matvec_t(ty)) err(talker,"agm of two vector/matrices");
    z=x; x=y; y=z;
  }
  if (gcmp0(y)) return gcopy(y);
  av=avma; z=sagm(gdiv(x,y),prec); tetpil=avma;
  return gerepile(av,tetpil,gmul(y,z));
}

GEN
logagm(GEN q)
{
  long av=avma,prec=lg(q),tetpil,s,n,lim;
  GEN y,q4,q1;

  if (typ(q)!=t_REAL) err(typeer,"logagm");
  if (signe(q)<=0) err(talker,"non positive argument in logagm");
  if (expo(q)<0) s=1; else { q=ginv(q); s=0; }
  lim = - (bit_accuracy(prec) >> 1);
  for (n=0; expo(q)>=lim ; n++) { q1=q; q=gsqr(q); }
  q4=gmul2n(q,2);
  if (!n) q1=gsqrt(q,prec);
  y=divrr(mppi(prec), agm(addsr(1,q4),gmul2n(q1,2),prec));
  tetpil=avma; y=gmul2n(y,-n); if (s) setsigne(y,-1);
  return gerepile(av,tetpil,y);
}

GEN
glogagm(GEN x, long prec)
{
  long av,tetpil;
  GEN y,p1,p2;

  switch(typ(x))
  {
    case t_REAL:
      if (signe(x)>=0) return logagm(x);

      y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
      setsigne(x,1); y[1]=(long)logagm(x);
      setsigne(x,-1); return y;

    case t_COMPLEX:
      y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
      av=avma; p1=glogagm(gnorm(x),prec); tetpil=avma;
      y[1]=lpile(av,tetpil,gmul2n(p1,-1));
      return y;

    case t_PADIC:
      return palog(x);

    case t_SER:
      av=avma; if (valp(x)) err(negexper,"glogagm");
      p1=gdiv(derivser(x),x);
      p1=integ(p1,varn(x));
      if (gcmp1((GEN)x[2])) return gerepileupto(av,p1);
      p2=glogagm((GEN)x[2],prec); tetpil=avma;
      return gerepile(av,tetpil,gadd(p1,p2));

    case t_INTMOD:
      err(typeer,"glogagm");
  }
  return transc(glogagm,x,prec);
}

GEN
theta(GEN q, GEN z, long prec)
{
  long av=avma,tetpil,l,n;
  GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold;

  if (gexpo(q)>=0) err(thetaer1);
  l=precision(q); if (l) prec=l;
  p1=realun(prec); z=gmul(p1,z);
  if (!l) q=gmul(p1,q);
  zy = gimag(z);
  if (gcmp0(zy)) k=gzero;
  else
  {
    lq=glog(q,prec); k=ground(gdiv(zy,greal(lq)));
    if (!gcmp0(k)) { zold=z; z=gadd(z,gdiv(gmul(lq,k),gi)); }
  }
  y=gsin(z,prec); n=0; qn=gun;
  ps2=gsqr(q); ps=gneg_i(ps2);
  do
  {
    n++; p1=gsin(gmulsg(2*n+1,z),prec);
    qnold=qn; qn=gmul(qn,ps);
    ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
  }
  while (gexpo(qnold) >= -bit_accuracy(prec));
  if (signe(k))
  {
    y=gmul(y,gmul(gpui(q,gsqr(k),prec),
                  gexp(gmul2n(gmul(gmul(gi,zold),k),1),prec)));
    if (mod2(k)) y=gneg_i(y);
  }
  p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
  return gerepile(av,tetpil,gmul(p1,y));
}

GEN
thetanullk(GEN q, long k, long prec)
{
  long av=avma,tetpil,l,n;
  GEN p1,ps,qn,y,ps2;

  if (gexpo(q)>=0) err(thetaer1);
  if (!(k&1)) return gzero;
  n=0; qn=gun; ps2=gsqr(q); ps=gneg_i(ps2);
  y=gun; l=precision(q);
  if (!l) { l=prec; q=gmul(q,realun(l)); }

  do
  {
    n++; p1=gpuigs(stoi(n+n+1),k); qn=gmul(qn,ps);
    ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
  }
  while (gexpo(p1) >= -bit_accuracy(l));
  p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
  if (k&2) { p1=gneg_i(p1); tetpil=avma; }
  return gerepile(av,tetpil,gmul(p1,y));
}

GEN
jbesselh(GEN n, GEN z, long prec)
{
  long av,tetpil,k,l,i,lz;
  GEN y,p1,p2,zinv,p0,s,c;

  if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh");
  k=itos(n); if (k<0) err(impl,"ybesselh");

  switch(typ(z))
  {
    case t_REAL: case t_COMPLEX:
      if (gcmp0(z)) return gzero;
      av=avma; zinv=ginv(z);
      l=precision(z); if (l>prec) prec=l;
      gsincos(z,&s,&c,prec);
      p1=gmul(zinv,s);
      if (k)
      {
        p0=p1; p1=gmul(zinv,gsub(p0,c));
	for (i=2; i<=k; i++)
	{
	  p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0);
          p0=p1; p1=p2;
	}
      }
      p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec);
      tetpil=avma; return gerepile(av,tetpil,gmul(p2,p1));

    case t_VEC: case t_COL: case t_MAT:
      lz=lg(z); y=cgetg(lz,typ(z));
      for (i=1; i<lz; i++)
	y[i]=(long)jbesselh(n,(GEN)z[i],prec);
      return y;

    case t_INT: case t_FRAC: case t_FRACN:
      av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
      return gerepile(av,tetpil,jbesselh(n,p1,prec));

    case t_QUAD:
      av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
      return gerepile(av,tetpil,jbesselh(n,p1,prec));

    case t_POL: case t_RFRAC: case t_RFRACN:
      av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
      return gerepile(av,tetpil,jbesselh(n,p1,prec));

    case t_POLMOD:
      av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
      for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
      tetpil=avma; y=cgetg(lz,t_COL);
      for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec);
      return gerepile(av,tetpil,y);

    case t_PADIC:
      err(impl,"p-adic jbessel function");
    case t_SER:
      err(impl,"jbessel of power series");
  }
  err(typeer,"jbesselh");
  return NULL; /* not reached */
}