/********************************************************************/
/**                                                                **/
/**                       ELLIPTIC CURVES                          **/
/**                                                                **/
/********************************************************************/
/* $Id: elliptic.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
#include "pari.h"

void
checkpt(GEN z)
{
  if (typ(z)!=t_VEC) err(elliper1);
}

long
checkell(GEN e)
{
  long lx=lg(e);
  if (typ(e)!=t_VEC || lx<14) err(elliper1);
  return lx;
}

void
checkbell(GEN e)
{
  if (typ(e)!=t_VEC || lg(e)<20) err(elliper1);
}

void
checksell(GEN e)
{
  if (typ(e)!=t_VEC || lg(e)<6) err(elliper1);
}

static void
checkch(GEN z)
{
  if (typ(z)!=t_VEC || lg(z)!=5) err(elliper1);
}

/* 4 X^3 + b2 X^2 + 2b4 X + b6 */
static GEN
RHSpol(GEN e)
{
  GEN z = cgetg(6, t_POL); z[1] = evalsigne(1)|evallgef(6);
  z[2] = e[8];
  z[3] = lmul2n((GEN)e[7],1);
  z[4] = e[6];
  z[5] = lstoi(4); return z;
}

/* x^3 + a2 x^2 + a4 x + a6 */
static GEN
ellRHS(GEN e, GEN x)
{
  GEN p1;
  p1 = gadd((GEN)e[2],x);
  p1 = gadd((GEN)e[4], gmul(x,p1));
  p1 = gadd((GEN)e[5], gmul(x,p1));
  return p1;
}

/* a1 x + a3 */
static GEN
ellLHS0(GEN e, GEN x)
{
  return gcmp0((GEN)e[1])? (GEN)e[3]: gadd((GEN)e[3], gmul(x,(GEN)e[1]));
}

static GEN
ellLHS0_i(GEN e, GEN x)
{
  return signe(e[1])? addii((GEN)e[3], mulii(x, (GEN)e[1])): (GEN)e[3];
}

/* y^2 + a1 xy + a3 y */
static GEN
ellLHS(GEN e, GEN z)
{
  GEN y = (GEN)z[2];
  return gmul(y, gadd(y, ellLHS0(e,(GEN)z[1])));
}

/* 2y + a1 x + a3 */
static GEN
d_ellLHS(GEN e, GEN z)
{
  return gadd(ellLHS0(e, (GEN)z[1]), gmul2n((GEN)z[2],1));
}


static void
smallinitell0(GEN x, GEN y)
{
  GEN b2,b4,b6,b8,d,j,a11,a13,a33,a64,b81,b22,c4,c6;
  long i;

  checksell(x); for (i=1; i<=5; i++) y[i]=x[i];

  b2=gadd(a11=gsqr((GEN)y[1]),gmul2n((GEN)y[2],2));
  y[6]=(long)b2;

  b4=gadd(a13=gmul((GEN)y[1],(GEN)y[3]),gmul2n((GEN)y[4],1));
  y[7]=(long)b4;

  b6=gadd(a33=gsqr((GEN)y[3]),a64=gmul2n((GEN)y[5],2));
  y[8]=(long)b6;

  b81=gadd(gadd(gmul(a11,(GEN)y[5]),gmul(a64,(GEN)y[2])),gmul((GEN)y[2],a33));
  b8=gsub(b81,gmul((GEN)y[4],gadd((GEN)y[4],a13)));
  y[9]=(long)b8;

  c4=gadd(b22=gsqr(b2),gmulsg(-24,b4));
  y[10]=(long)c4;

  c6=gadd(gmul(b2,gsub(gmulsg(36,b4),b22)),gmulsg(-216,b6));
  y[11]=(long)c6;

  b81=gadd(gmul(b22,b8),gmulsg(27,gsqr(b6)));
  d=gsub(gmul(b4,gadd(gmulsg(9,gmul(b2,b6)),gmulsg(-8,gsqr(b4)))),b81);
  y[12]=(long)d;

  if (gcmp0(d)) err(talker,"singular curve in ellinit");

  j = gdiv(gmul(gsqr(c4),c4),d);
  y[13]=(long)j;
}

GEN
smallinitell(GEN x)
{
  GEN y;
  long av,tetpil;

  av=avma; y=cgetg(14,t_VEC);
  smallinitell0(x,y); tetpil=avma;
  return gerepile(av,tetpil,gcopy(y));
}

GEN
ellinit0(GEN x, long flag,long prec)
{
  switch(flag)
  {
    case 0: return initell(x,prec);
    case 1: return smallinitell(x);
    default: err(flagerr,"ellinit");
  }
  return NULL; /* not reached */
}

void
ellprint(GEN e)
{
  long av = avma;
  long vx = fetch_var();
  long vy = fetch_var();
  GEN z = cgetg(3,t_VEC);
  if (typ(e) != t_VEC || lg(e) < 6)
    err(talker, "not an elliptic curve in ellprint");
  z[1] = lpolx[vx]; name_var(vx, "X");
  z[2] = lpolx[vy]; name_var(vy, "Y");
  fprintferr("%Z = %Z\n", ellLHS(e, z), ellRHS(e, polx[vx]));
  (void)delete_var();
  (void)delete_var(); avma = av;
}

static GEN
do_agm(GEN *ptx1, GEN a1, GEN b1, long prec, long sw)
{
  GEN p1,r1,a,b,x,x1;
  long G;

  x1 = gmul2n(gsub(a1,b1),-2);
  if (gcmp0(x1))
    err(talker,"precision too low in initell");
  G = 6 - bit_accuracy(prec);
  for(;;)
  {
    a=a1; b=b1; x=x1;
    b1=gsqrt(gmul(a,b),prec); setsigne(b1,sw);
    a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2);
    r1=gsub(a1,b1);
    p1=gsqrt(gdiv(gadd(x,r1),x),prec);
    x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
    if (gexpo(r1) <= G + gexpo(b1)) break;
  }
  if (gprecision(x1)*2 <= (prec+2))
    err(talker,"precision too low in initell");
  *ptx1 = x1; return ginv(gmul2n(a1,2));
}

static GEN
do_padic_agm(GEN *ptx1, GEN a1, GEN b1, GEN p)
{
  GEN p1,r1,a,b,x,bmod1, bmod = modii((GEN)b1[4],p), x1 = *ptx1;

  if (!x1) x1 = gmul2n(gsub(a1,b1),-2);
  for(;;)
  {
    a=a1; b=b1; x=x1;
    b1=gsqrt(gmul(a,b),0); bmod1=modii((GEN)b1[4],p);
    if (!egalii(bmod1,bmod)) b1 = gneg_i(b1);
    a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2);
    r1=gsub(a1,b1);
    p1=gsqrt(gdiv(gadd(x,r1),x),0);
    if (! gcmp1(modii((GEN)p1[4],p))) p1 = gneg_i(p1);
    x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
    if (gcmp0(r1)) break;
  }
  *ptx1 = x1; return ginv(gmul2n(a1,2));
}

static GEN
padic_initell(GEN y, GEN p, long prec)
{
  GEN b2,b4,c4,c6,p1,p2,w,pv,a1,b1,x1,u2,q,e0,e1;
  long i,alpha;

  if (valp(y[13]) >= 0) /* p | j */
    err(talker,"valuation of j must be negative in p-adic ellinit");
  if (egalii(p,gdeux))
    err(impl,"initell for 2-adic numbers"); /* pv=stoi(4); */

  pv=p; q=ggrandocp(p,prec);
  for (i=1; i<=5; i++) y[i]=ladd(q,(GEN)y[i]);
  b2= (GEN)y[6];
  b4= (GEN)y[7];
  c4= (GEN)y[10];
  c6= (GEN)y[11];
  alpha=valp(c4)>>1;
  setvalp(c4,0);
  setvalp(c6,0); e1=gdivgs(gdiv(c6,c4),6);
  c4=gdivgs(c4,48); c6=gdivgs(c6,864);
  do
  {
    e0=e1; p2=gsqr(e0);
    e1=gdiv(gadd(gmul2n(gmul(e0,p2),1),c6), gsub(gmulsg(3,p2),c4));
  }
  while (!gegal(e0,e1));
  setvalp(e1,valp(e1)+alpha);

  e1=gsub(e1,gdivgs(b2,12));
  w=gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),0);

  p1=gaddgs(gdiv(gmulsg(3,e0),w),1);
  if (valp(p1)<=0) w=gneg_i(w);
  y[18]=(long)w;

  a1=gmul2n(gsub(w,gadd(gmulsg(3,e1),gmul2n(b2,-2))),-2);
  b1=gmul2n(w,-1); x1=NULL;
  u2 = do_padic_agm(&x1,a1,b1,pv);

  w=gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
  q=ginv(gadd(w,gsqrt(gaddgs(gsqr(w),-1),0)));
  if (valp(q)<0) q=ginv(q);

  p1=cgetg(2,t_VEC); p1[1]=(long)e1;
  y[14]=(long)p1;
  y[15]=(long)u2;
  y[16] = (kronecker((GEN)u2[4],p) <= 0 || (valp(u2)&1))? zero: lsqrt(u2,0);
  y[17]=(long)q;
  y[19]=zero; return y;
}

static int
invcmp(GEN x, GEN y) { return -gcmp(x,y); }

static GEN
initell0(GEN x, long prec)
{
  GEN y,b2,b4,D,p1,p2,p,u,w,a1,b1,x1,u2,q,e1,pi,pi2,tau,w2;
  long ty,i,e,sw;

  y=cgetg(20,t_VEC);
  smallinitell0(x,y);

  e = BIGINT; p = NULL;
  for (i=1; i<=5; i++)
  {
    q = (GEN)y[i];
    if (typ(q)==t_PADIC)
    {
      long e2 = signe(q[4])? precp(q)+valp(q): valp(q);
      if (e2 < e) e = e2;
      if (!p) p = (GEN)q[2];
      else if (!egalii(p,(GEN)q[2]))
        err(talker,"incompatible p-adic numbers in initell");
    }
  }
  if (e<BIGINT) return padic_initell(y,p,e);

  b2= (GEN)y[6];
  b4= (GEN)y[7];
  D = (GEN)y[12]; ty = typ(D);
  if (!prec || !is_const_t(ty) || ty==t_INTMOD)
    { y[14]=y[15]=y[16]=y[17]=y[18]=y[19]=zero; return y; }

  pi=mppi(prec); pi2=gmul2n(pi,1);
  p1 = roots(RHSpol(y),prec);
  if (gsigne(D) < 0) p1[1] = lreal((GEN)p1[1]);
  else /* sort roots in decreasing order */
    p1 = gen_sort(greal(p1), 0, invcmp);
  y[14]=(long)p1;

  e1=(GEN)p1[1];
  w = gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
  p2=gadd(gmulsg(3,e1), gmul2n(b2,-2));
  if (gsigne(p2)>0) w = gneg_i(w);
  a1=gmul2n(gsub(w,p2),-2);
  b1=gmul2n(w,-1); sw=signe(w);
  u2 = do_agm(&x1,a1,b1,prec,sw);

  w = gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
  q = gsqrt(gaddgs(gsqr(w),-1),prec);
  if (gsigne(greal(w))>0)
    q = ginv(gadd(w,q));
  else
    q = gsub(w,q);
  if (gexpo(q)>=0) q=ginv(q);
  tau=gmul(gdiv(glog(q,prec),pi2), gneg_i(gi));

  y[19]=lmul(gmul(gsqr(pi2),gabs(u2,prec)),gimag(tau));
  u=gmul(pi2,gsqrt(gneg_i(u2),prec)); w2=gmul(tau,u);
  if (sw<0)
  {
    y[15]=(long)u;
    q=gsqrt(q,prec);
  }
  else
  {
    y[15]=lmul2n(gabs((GEN)w2[1],prec),1);
    q=gexp(gmul2n(gmul(gmul(pi2,gi),gdiv(w2,(GEN)y[15])), -1), prec);
  }
  y[16]=(long)w2;
  p1 = gdiv(gsqr(pi),gmulsg(6,(GEN)y[15]));
  p2 = thetanullk(q,1,prec);
  if (gcmp0(p2)) err(talker,"precision too low in initell");
  y[17]=lmul(p1,gdiv(thetanullk(q,3,prec),p2));
  y[18]=ldiv(gsub(gmul((GEN)y[17],(GEN)y[16]),gmul(gi,pi)),(GEN)y[15]);
  return y;
}

GEN
initell(GEN x, long prec)
{
  long av=avma;
  return gerepileupto(av, gcopy(initell0(x,prec)));
}

GEN
coordch(GEN e, GEN ch)
{
  GEN y,p1,p2,v,v2,v3,v4,v6,r,s,t,u;
  long av,tetpil,i,lx = checkell(e);

  checkch(ch);
  u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];
  av=avma; y=cgetg(lx,t_VEC);
  v=ginv(u); v2=gsqr(v); v3=gmul(v,v2);v4=gsqr(v2); v6=gsqr(v3);
  y[1] = lmul(v,gadd((GEN)e[1],gmul2n(s,1)));
  y[2] = lmul(v2,gsub(gadd((GEN)e[2],gmulsg(3,r)),gmul(s,gadd((GEN)e[1],s))));
  p2 = ellLHS0(e,r);
  p1 = gadd(gmul2n(t,1), p2);
  y[3] = lmul(v3,p1);
  p1 = gsub((GEN)e[4],gadd(gmul(t,(GEN)e[1]),gmul(s,p1)));
  y[4] = lmul(v4,gadd(p1,gmul(r,gadd(gmul2n((GEN)e[2],1),gmulsg(3,r)))));
  p2 = gmul(t,gadd(t, p2));
  y[5] = lmul(v6,gsub(ellRHS(e,r),p2));
  y[6] = lmul(v2,gadd((GEN)e[6],gmulsg(12,r)));
  y[7] = lmul(v4,gadd((GEN)e[7],gmul(r,gadd((GEN)e[6],gmulsg(6,r)))));
  y[8] = lmul(v6,gadd((GEN)e[8],gmul(r,gadd(gmul2n((GEN)e[7],1),gmul(r,gadd((GEN)e[6],gmul2n(r,2)))))));
  p1 = gadd(gmulsg(3,(GEN)e[7]),gmul(r,gadd((GEN)e[6],gmulsg(3,r))));
  y[9] = lmul(gsqr(v4),gadd((GEN)e[9],gmul(r,gadd(gmulsg(3,(GEN)e[8]),gmul(r,p1)))));
  y[10] = lmul(v4,(GEN)e[10]);
  y[11] = lmul(v6,(GEN)e[11]);
  y[12] = lmul(gsqr(v6),(GEN)e[12]);
  y[13] = e[13];
  if (lx>14)
  {
    p1=(GEN)e[14];
    if (gcmp0(p1))
    {
      y[14] = y[15] = y[16] = y[17] = y[18] = y[19] = zero;
    }
    else
    {
      if (typ(e[1])==t_PADIC)
      {
        p2=cgetg(2,t_VEC); p2[1]=lmul(v2,gsub((GEN)p1[1],r));
        y[14]=(long)p2;
        y[15]=lmul(gsqr(u),(GEN)e[15]);
        y[16]=lmul(u,(GEN)e[16]);
/* FIXME: how do q and w change ??? */
        y[17]=e[17];
        y[18]=e[18];
        y[19]=zero;
      }
      else
      {
        p2=cgetg(4,t_COL);
        for (i=1; i<=3; i++) p2[i]=lmul(v2,gsub((GEN)p1[i],r));
        y[14]=(long)p2;
        y[15]=lmul(u,(GEN)e[15]);
        y[16]=lmul(u,(GEN)e[16]);
        y[17]=ldiv((GEN)e[17],u);
        y[18]=ldiv((GEN)e[18],u);
        y[19]=lmul(gsqr(u),(GEN)e[19]);
      }
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

static GEN
pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t)
{
  GEN p1,z;

  if (lg(x) < 3) return x;

  z = cgetg(3,t_VEC); p1=gadd((GEN)x[1],mor);
  z[1] = lmul(v2,p1);
  z[2] = lmul(v3,gsub((GEN)x[2],gadd(gmul(s,p1),t)));
  return z;
}

GEN
pointch(GEN x, GEN ch)
{
  GEN y,v,v2,v3,mor,r,s,t,u;
  long av,tetpil,tx,lx=lg(x),i;

  checkpt(x); checkch(ch);
  if (lx < 2) return gcopy(x);
  av=avma; u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];
  tx=typ(x[1]); v=ginv(u); v2=gsqr(v); v3=gmul(v,v2); mor=gneg_i(r);
  if (is_matvec_t(tx))
  {
    y=cgetg(lx,tx);
    for (i=1; i<lx; i++)
      y[i]=(long) pointch0((GEN)x[i],v2,v3,mor,s,t);
  }
  else
    y = pointch0(x,v2,v3,mor,s,t);
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

static long
ellexpo(GEN e)
{
  long i,k2, k = gexpo((GEN)e[1]);
  for (i=2; i<6; i++)
  {
    k2 = gexpo((GEN)e[i]);
    if (k<k2) k = k2;
  }
  return k;
}

/* Exactness of lhs and rhs in the following depends in non-obvious ways
   on the coeffs of the curve as well as on the components of the point z.
   Thus if e is exact, with a1==0, and z has exact y coordinate only, the
   lhs will be exact but the rhs won't. */
int
oncurve(GEN e, GEN z)
{
  GEN p1,p2,x;
  long av=avma,p,q;

  checksell(e); checkpt(z); if (lg(z)<3) return 1; /* oo */
  p1 = ellLHS(e,z);
  p2 = ellRHS(e,(GEN)z[1]); x = gsub(p1,p2);
  if (gcmp0(x)) { avma=av; return 1; }
  p = precision(p1);
  q = precision(p2);
  if (!p && !q) { avma=av; return 0; } /* both of p1, p2 are exact */
  if (!q || (p && p < q)) q = p; /* min among nonzero elts of {p,q} */
  /* the constant 0.93 is arbitrary */
  q = (gexpo(x)-ellexpo(e) < - bit_accuracy(q) * 0.93);
  avma = av; return q;
}

GEN
addell(GEN e, GEN z1, GEN z2)
{
  GEN p1,p2,x,y,x1,x2,y1,y2;
  long av=avma,tetpil;

  checksell(e); checkpt(z1); checkpt(z2);
  if (lg(z1)<3) return gcopy(z2);
  if (lg(z2)<3) return gcopy(z1);

  x1=(GEN)z1[1]; y1=(GEN)z1[2];
  x2=(GEN)z2[1]; y2=(GEN)z2[2];
  if (x1 == x2 || gegal(x1,x2))
  { /* y1 = y2 or -LHS0-y2 */
    if (y1 != y2)
    {
      int eq;
      if (precision(y1) || precision(y2))
        eq = (gexpo(gadd(ellLHS0(e,x1),gadd(y1,y2))) >= gexpo(y1));
      else
        eq = gegal(y1,y2);
      if (!eq) { avma=av; y=cgetg(2,t_VEC); y[1]=zero; return y; }
    }
    p2 = d_ellLHS(e,z1);
    if (gcmp0(p2)) { avma=av; y=cgetg(2,t_VEC); y[1]=zero; return y; }
    p1 = gadd(gsub((GEN)e[4],gmul((GEN)e[1],y1)),
              gmul(x1,gadd(gmul2n((GEN)e[2],1),gmulsg(3,x1))));
  }
  else { p1=gsub(y2,y1); p2=gsub(x2,x1); }
  p1 = gdiv(p1,p2);
  x = gsub(gmul(p1,gadd(p1,(GEN)e[1])), gadd(gadd(x1,x2),(GEN)e[2]));
  y = gadd(gadd(y1, ellLHS0(e,x)), gmul(p1,gsub(x,x1)));
  tetpil=avma; p1=cgetg(3,t_VEC); p1[1]=lcopy(x); p1[2]=lneg(y);
  return gerepile(av,tetpil,p1);
}

static GEN
invell(GEN e, GEN z)
{
  GEN p1;

  if (lg(z)<3) return z;
  p1=cgetg(3,t_VEC); p1[1]=z[1];
  p1[2]=(long)gneg_i(gadd((GEN)z[2], ellLHS0(e,(GEN)z[1])));
  return p1;
}

GEN
subell(GEN e, GEN z1, GEN z2)
{
  long av=avma,tetpil;

  checksell(e); checkpt(z2);
  z2=invell(e,z2); tetpil=avma;
  return gerepile(av,tetpil,addell(e,z1,z2));
}

GEN
ordell(GEN e, GEN x, long prec)
{
  long av=avma,td,i,lx,tx=typ(x);
  GEN D,a,b,d,pd,p1,y;

  checksell(e);
  if (is_matvec_t(tx))
  {
    lx=lg(x); y=cgetg(lx,tx);
    for (i=1; i<lx; i++) y[i]=(long)ordell(e,(GEN)x[i],prec);
    return y;
  }

  a=ellRHS(e,x);
  b=ellLHS0(e,x); /* y*(y+b) = a */
  D=gadd(gsqr(b),gmul2n(a,2)); td=typ(D);
  if (gcmp0(D))
  {
    b = gneg_i(b);
    y = cgetg(2,t_VEC);
    if (td == t_INTMOD && egalii((GEN)D[1], gdeux))
      y[1] = (long)gmodulss(gcmp0(a)?0:1, 2);
    else
      y[1] = lmul2n(b,-1);
    return gerepileupto(av,y);
  }

  if (td==t_INT || is_frac_t(td))
  {
    pd = (td==t_INT)? NULL: (GEN)D[2];
    if (pd) D = mulii((GEN)D[1],pd);
    if (!carrecomplet(D,&d)) { avma=av; return cgetg(1,t_VEC); }
    if (pd) d = gdiv(d,pd);
  }
  else
  {
    if (td==t_INTMOD)
    {
      if (egalii((GEN)D[1],gdeux))
      {
        avma=av;
        if (!gcmp0(a)) return cgetg(1,t_VEC);
        y = cgetg(3,t_VEC);
        y[1] = (long)gmodulss(0,2);
        y[2] = (long)gmodulss(1,2); return y;
      }
      if (kronecker((GEN)D[2],(GEN)D[1]) == -1)
        { avma=av; return cgetg(1,t_VEC); }
    }
    d = gsqrt(D,prec);
  }
  p1=gsub(d,b); y = cgetg(3,t_VEC);
  y[1] = lmul2n(p1,-1);
  y[2] = lsub((GEN)y[1],d);
  return gerepileupto(av,y);
}

static GEN
CM_powell(GEN e, GEN z, GEN n)
{
  GEN x,y,p0,p1,q0,q1,p2,q2,z1,z2,pol,grdx;
  long av=avma,tetpil,ln,ep,vn;

  if (lg(z)<3) return gcopy(z);
  pol=(GEN)n[1];
  if (signe(discsr(pol))>=0)
    err(talker,"not a negative quadratic discriminant in CM");
  if (!gcmp1(denom((GEN)n[2])) || !gcmp1(denom((GEN)n[3])))
    err(impl,"powell for nonintegral CM exponent");

  p1=gaddgs(gmul2n(gnorm(n),2),4);
  if (gcmpgs(p1,(((ulong)MAXULONG)>>1)) > 0)
    err(talker,"norm too large in CM");
  ln=itos(p1); vn=(ln-4)>>2;
  z1 = weipell(e,ln);
  z2 = gsubst(z1,0,gmul(n,polx[0]));
  grdx=gadd((GEN)z[1],gdivgs((GEN)e[6],12));
  p0=gzero; p1=gun;
  q0=gun; q1=gzero;
  do
  {
    GEN ss=gzero;
    do
    {
      ep=(-valp(z2))>>1; ss=gadd(ss,gmul((GEN)z2[2],gpuigs(polx[0],ep)));
      z2=gsub(z2,gmul((GEN)z2[2],gpuigs(z1,ep)));
    }
    while (valp(z2)<=0);
    p2=gadd(p0,gmul(ss,p1)); p0=p1; p1=p2;
    q2=gadd(q0,gmul(ss,q1)); q0=q1; q1=q2;
    if (!signe(z2)) break;
    z2=ginv(z2);
  }
  while (lgef(p1)-3 < vn);
  if (lgef(p1)-3 > vn || signe(z2))
    err(talker,"not a complex multiplication in powell");
  x=gdiv(p1,q1); y=gdiv(deriv(x,0),n);
  x=gsub(poleval(x,grdx), gdivgs((GEN)e[6],12));
  y=gsub(gmul(d_ellLHS(e,z),poleval(y,grdx)), ellLHS0(e,x));
  tetpil=avma; z=cgetg(3,t_VEC); z[1]=lcopy(x); z[2]=lmul2n(y,-1);
  return gerepile(av,tetpil,z);
}

GEN
powell(GEN e, GEN z, GEN n)
{
  GEN y;
  long av=avma,i,j,tetpil,s;
  ulong m;

  checksell(e); checkpt(z);
  if (typ(n)==t_QUAD) return CM_powell(e,z,n);
  if (typ(n)!=t_INT)
    err(impl,"powell for nonintegral or non CM exponents");
  if (lg(z)<3) return gcopy(z);
  s=signe(n);
  if (!s) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
  if (s < 0) { n=negi(n); z = invell(e,z); }
  if (is_pm1(n)) { tetpil=avma; return gerepile(av,tetpil,gcopy(z)); }

  y=cgetg(2,t_VEC); y[1]=zero;
  for (i=lgefint(n)-1; i>2; i--)
    for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
    {
      if (m&1) y = addell(e,y,z);
      z = addell(e,z,z);
    }
  for (m=n[2]; m>1; m>>=1)
  {
    if (m&1) y = addell(e,y,z);
    z = addell(e,z,z);
  }
  tetpil=avma; return gerepile(av,tetpil,addell(e,y,z));
}

GEN
mathell(GEN e, GEN x, long prec)
{
  GEN y,p1,p2, *pdiag;
  long av=avma,tetpil,lx=lg(x),i,j,tx=typ(x);

  if (!is_vec_t(tx)) err(elliper1);
  lx=lg(x); y=cgetg(lx,t_MAT); pdiag=(GEN*) new_chunk(lx);
  for (i=1; i<lx; i++)
  {
    pdiag[i]=ghell(e,(GEN)x[i],prec);
    y[i]=lgetg(lx,t_COL);
  }
  for (i=1; i<lx; i++)
  {
    p1=(GEN)y[i]; p1[i]=lmul2n(pdiag[i],1);
    for (j=i+1; j<lx; j++)
    {
      p2=ghell(e,addell(e,(GEN)x[i],(GEN)x[j]),prec);
      p2=gsub(p2, gadd(pdiag[i],pdiag[j]));
      p1[j]=(long)p2; coeff(y,i,j)=(long)p2;
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

static GEN
bilhells(GEN e, GEN z1, GEN z2, GEN h2, long prec)
{
  long lz1=lg(z1),tx,av=avma,tetpil,i;
  GEN y,p1,p2;

  if (lz1==1) return cgetg(1,typ(z1));

  tx=typ(z1[1]);
  if (!is_matvec_t(tx))
  {
    p1 = ghell(e,addell(e,z1,z2),prec);
    p2 = gadd(ghell(e,z1,prec),h2);
    tetpil=avma; return gerepile(av,tetpil,gsub(p1,p2));
  }
  y=cgetg(lz1,typ(z1));
  for (i=1; i<lz1; i++)
    y[i]=(long)bilhells(e,(GEN)z1[i],z2,h2,prec);
  return y;
}

GEN
bilhell(GEN e, GEN z1, GEN z2, long prec)
{
  GEN p1,h2;
  long av=avma,tetpil,tz1=typ(z1),tz2=typ(z2);

  if (!is_matvec_t(tz1) || !is_matvec_t(tz2)) err(elliper1);
  if (lg(z1)==1) return cgetg(1,tz1);
  if (lg(z2)==1) return cgetg(1,tz2);

  tz1=typ(z1[1]); tz2=typ(z2[1]);
  if (is_matvec_t(tz2))
  {
    if (is_matvec_t(tz1))
      err(talker,"two vector/matrix types in bilhell");
    p1=z1; z1=z2; z2=p1;
  }
  h2=ghell(e,z2,prec); tetpil=avma;
  return gerepile(av,tetpil,bilhells(e,z1,z2,h2,prec));
}

static GEN
new_coords(GEN e, GEN x, GEN *pta, GEN *ptb, long prec)
{
  GEN a,b,r0,r1,p1,p2,w, e1 = gmael(e,14,1), b2 = (GEN)e[6];
  long ty = typ(e[12]);

  r0 = gmul2n(b2,-2);
  p2 = gadd(gmulsg(3,e1),r0);
  if (ty == t_PADIC)
    w = (GEN)e[18];
  else
  {
    GEN b4 = (GEN)e[7];
    if (!is_const_t(ty)) err(typeer,"zell");

    w = gsqrt(gmul2n(gadd(b4, gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
    if (gsigne(greal(p2)) > 0) w = gneg_i(w);
  }
  a = gmul2n(gsub(w,p2),-2);
  b = gmul2n(w,-1);
  r1 = gsub(a,b);
  p1 = gadd(x, gmul2n(gadd(e1,r0),-1));
  if (gcmp0(p1)) p1=gsqrt(gneg_i(gmul(a,r1)),prec);
  else
  {
    GEN delta=gdiv(gmul(a,r1),gsqr(p1));
    p1=gmul2n(gmul(p1,gaddsg(1,gsqrt(gsubsg(1,gmul2n(delta,2)),prec))),-1);
  }
  *pta = a; *ptb = b;
  return gmul(p1,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(p1,r1),p1),prec)),-1)));
}

GEN
zell(GEN e, GEN z, long prec)
{
  long av=avma,ty,sw,fl;
  GEN t,u,p1,p2,r1,a,b,x1,u2,D = (GEN)e[12];

  checkbell(e);
  if (!oncurve(e,z)) err(heller1);
  ty=typ(D);
  if (ty==t_INTMOD) err(typeer,"zell");
  if (lg(z)<3) return (ty==t_PADIC)? gun: gzero;

  x1 = new_coords(e,(GEN)z[1],&a,&b,prec);
  if (ty==t_PADIC)
  {
    u2 = do_padic_agm(&x1,a,b,(GEN)D[2]);
    if (!gcmp0((GEN)e[16]))
    {
      t=gsqrt(gaddsg(1,gdiv(x1,a)),prec);
      t=gdiv(gaddsg(-1,t),gaddsg(1,t));
    }
    else t=gaddsg(2,ginv(gmul(u2,x1)));
    return gerepileupto(av,t);
  }

  sw=gsigne(greal(b)); fl=0;
  for(;;) /* agm */
  {
    GEN a0=a, b0=b, x0=x1;

    b=gsqrt(gmul(a0,b0),prec);
    if (gsigne(greal(b)) != sw) b = gneg_i(b);
    a=gmul2n(gadd(gadd(a0,b0),gmul2n(b,1)),-2);
    r1=gsub(a,b);
    p1=gsqrt(gdiv(gadd(x0,r1),x0),prec);
    x1=gmul(x0,gsqr(gmul2n(gaddsg(1,p1),-1)));
    if (gexpo(gsub(x1,x0)) < gexpo(x1) - bit_accuracy(prec) + 5)
    {
      if (fl) break;
      fl = 1;
    }
    else fl = 0;
  }
  u=gdiv(x1,a); t=gaddsg(1,u);
  if (gexpo(t) <  5 - bit_accuracy(prec))
    t = negi(gun);
  else
    t = gdiv(u,gsqr(gaddsg(1,gsqrt(t,prec))));
  u=gsqrt(ginv(gmul2n(a,2)),prec);
  t=glog(t,prec); t=gmul(t,u);

  /* which square root? test the reciprocal function (pointell) */
  if (!gcmp0(t))
  {
    GEN x1;
    long bad;

    u = pointell(e,t,3); /* we don't need much precision */
    /* Either z = u (ok: keep t), or z = invell(e,u) (bad: t <-- -t) */
    x1 = gsub(z,u); bad = (gexpo((GEN)x1[1]) >= gexpo((GEN)u[1])
                        || gexpo((GEN)x1[2]) >= gexpo((GEN)u[2]));
    if (bad) t = gneg(t);
    if (DEBUGLEVEL)
    {
      if (DEBUGLEVEL>4)
      {
        fprintferr("  z  = %Z\n",z);
        fprintferr("  u  = %Z\n",u);
        fprintferr("  x1 = %Z\n",x1);
      }
      fprintferr("ellpointtoz: %s square root\n", bad? "bad": "good");
      flusherr();
    }
  }
  /* send t to the fundamental domain if necessary */
  p2 = gdiv(gimag(t),gmael(e,16,2));
  p1 = gsub(p2, gmul2n(gun,-2));
  if (gcmp(gabs(p1,prec),ghalf) >= 0)
    t = gsub(t, gmul((GEN)e[16],gfloor(gadd(p2,dbltor(0.1)))));
  if (gsigne(greal(t)) < 0) t = gadd(t,(GEN)e[15]);
  return gerepileupto(av,t);
}

/* compute gamma in SL_2(Z) and t'=gamma(t) so that t' is in the usual
   fundamental domain. Internal function no check, no garbage. */
static GEN
getgamma(GEN t)
{
  GEN a,b,c,d,n,m,y,p1,unapprox;

  unapprox=gsub(gun,gpuigs(stoi(10),-8));
  a=d=gun;b=c=gzero;
  for(;;)
  {
    n=ground(greal(t));
    if (signe(n))
    {
      n=negi(n); t=gadd(t,n);
      a=addii(a,mulii(n,c));
      b=addii(b,mulii(n,d));
    }
    m=gnorm(t); if (gcmp(m,unapprox)>=0) break;
    t=gneg_i(gdiv(gconj(t),m));
    p1=negi(c); c=a; a=p1;
    p1=negi(d); d=b; b=p1;
  }
  y=cgetg(3,t_VEC);
  m=cgetg(3,t_MAT); y[1]=(long)m;
  p1=cgetg(3,t_COL); m[1]=(long)p1;
  p1[1]=(long)a; p1[2]=(long)c;
  p1=cgetg(3,t_COL); m[2]=(long)p1;
  p1[1]=(long)b; p1[2]=(long)d;
  y[2]=(long)t;
  return y;
}

static int
get_periods(GEN e, GEN *om1, GEN *om2)
{
  long tx = typ(e);
  if (is_vec_t(tx))
    switch(lg(e))
    {
      case  3: *om1=(GEN)e[1];  *om2=(GEN)e[2]; return 1;
      case 20: *om1=(GEN)e[16]; *om2=(GEN)e[15]; return 1;
    }
  return 0;
}

/* computes the numerical values of eisenstein series. k is equal to a positive
   even integer. If k=4 or 6, computes g2 or g3. If k=2, or k>6 even,
   compute (2iPi/om2)^k*(1+2/zeta(1-k)*sum(n>=1,n^(k-1)q^n/(1-q^n)) with no
   constant factor. */
GEN
elleisnum(GEN om, long k, long flag, long prec)
{
  long av=avma,lim,av1,fl,si;
  GEN om1,om2,p1,pii2,tau,q,y,qn,v,ga,court,asub;

  if (k%2 || k<=0) err(talker,"k not a positive even integer in elleisnum");
  if (!get_periods(om, &om1, &om2)) err(typeer,"elleisnum");
  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  tau=gdiv(om1,om2); si=gsigne(gimag(tau));
  if (si==0)
    err(talker,"omega1 and omega2 are R-linearly dependent in elleisnum");
  if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
  v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
  if (k==2) asub=gdiv(gmul(pii2,gmulsg(12,gcoeff(ga,2,1))),om2);
  om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
  if (k==2) asub=gdiv(asub,om2);
  q=gexp(gmul(pii2,tau),prec);
  y=gzero; court=stoi(3);
  av1=avma; lim=stack_lim(av1,1); qn=gun; court[2]=0;
  do
  {
    court[2]++; qn=gmul(q,qn);
    p1=gdiv(gmul(gpuigs(court,k-1),qn),gsub(gun,qn));
    y=gadd(y,p1);
    fl=(gexpo(p1) > - bit_accuracy(prec) - 5);
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
      if(DEBUGMEM>1) err(warnmem,"elleisnum");
      gerepilemany(av1,gptr,2);
    }
  }
  while (fl);

  y=gadd(gun,gmul(gdiv(gdeux,gzeta(stoi(1-k),prec)),y));
  p1=gpuigs(gdiv(pii2,om2),k);
  y = gmul(p1,y);
  if (k==2) y=gsub(y,asub);
  else if (k==4 && flag) y=gdivgs(y,12);
  else if (k==6 && flag) y=gdivgs(y,216);
  return gerepileupto(av,y);
}

/* compute eta1, eta2 */

GEN
elleta(GEN om, long prec)
{
  long av=avma,tetpil;
  GEN e2,p1,pii2,y1,y2,y;

  e2=gdivgs(elleisnum(om,2,0,prec),12);
  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  y2=gmul((GEN)om[2],e2);
  y1=gadd(gdiv(pii2,(GEN)om[2]),gmul((GEN)om[1],e2));
  tetpil=avma;
  y=cgetg(3,t_VEC); y[1]=lneg(y1); y[2]=lneg(y2);
  return gerepile(av,tetpil,y);
}

/* computes the numerical value of wp(z | om1 Z + om2 Z),
   If flall=1, compute also wp'. Reduce to the fundamental domain first. */
static GEN
weipellnumall(GEN om1, GEN om2, GEN z, long flall, long prec)
{
  long av=avma,tetpil,lim,av1,si,toadd;
  GEN p1,pii2,pii4,a,tau,q,u,y,yp,u1,u2,qn,v,ga;

  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  tau=gdiv(om1,om2); si=gsigne(gimag(tau));
  if (si==0)
    err(talker,"omega1 and omega2 are R-linearly dependent in ellwpnum");
  if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
  v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
  om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
  z=gdiv(z,om2);
  a=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(a,tau));
  a=ground(greal(z)); z=gsub(z,a);
  if (gexpo(z) < 5 - bit_accuracy(prec))
  {
    avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v;
  }

  q=gexp(gmul(pii2,tau),prec);
  u=gexp(gmul(pii2,z),prec);
  u1=gsub(gun,u); u2=gsqr(u1);
  y=gadd(gdivgs(gun,12),gdiv(u,u2));
  if (flall) yp=gdiv(gadd(gun,u),gmul(u1,u2));
  toadd=(long)ceil(9.065*gtodouble(gimag(z)));

  av1=avma; lim=stack_lim(av1,1); qn=q;
  do
  {
    GEN p2,qnu,qnu1,qnu2,qnu3,qnu4;

    qnu=gmul(qn,u); qnu1=gsub(gun,qnu); qnu2=gsqr(qnu1);
    qnu3=gsub(qn,u); qnu4=gsqr(qnu3);
    p1=gsub(gmul(u,gadd(ginv(qnu2),ginv(qnu4))),
	    gmul2n(ginv(gsqr(gsub(gun,qn))),1));
    p1=gmul(qn,p1);
    y=gadd(y,p1);
    if (flall)
    {
    p2=gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)),
	    gdiv(gadd(qn,u),gmul(qnu3,qnu4)));
    p2=gmul(qn,p2);
      yp=gadd(yp,p2);
    }
    qn=gmul(q,qn);
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[3]; gptr[0]=&y; gptr[1]=&qn; gptr[2]=&yp;
      if(DEBUGMEM>1) err(warnmem,"weipellnum");
      gerepilemany(av1,gptr,flall?3:2);
    }
  }
  while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);

  pii2=gdiv(pii2,om2);
  pii4=gsqr(pii2);
  y = gmul(pii4,y);
  if (flall) yp=gmul(u,gmul(gmul(pii4,pii2),yp));
  tetpil=avma;
  if (flall) { v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lmul2n(yp,-1); }
  else v=gcopy(y);
  return gerepile(av,tetpil,v);
}

GEN
ellzeta(GEN om, GEN z, long prec)
{
  long av=avma,tetpil,lim,av1,si,toadd;
  GEN zinit,om1,om2,p1,pii2,tau,q,u,y,u1,qn,v,ga,x1,x2,et;

  if (!get_periods(om, &om1, &om2)) err(typeer,"ellzeta");
  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  tau=gdiv(om1,om2); si=gsigne(gimag(tau));
  if (si==0)
    err(talker,"omega1 and omega2 are R-linearly dependent in ellzeta");
  if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
  v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
  om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
  om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
  z=gdiv(z,om2);

  x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
  x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
  et=elleta(om,prec);
  et=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
  if (gexpo(z) < 5 - bit_accuracy(prec))
  {
    p1=ginv(zinit); tetpil=avma; return gerepile(av,tetpil,gadd(p1,et));
  }
  q=gexp(gmul(pii2,tau),prec);
  u=gexp(gmul(pii2,z),prec);
  u1=gsub(u,gun);
  y=gdiv(gmul(gsqr(om2),elleisnum(om,2,0,prec)),pii2);
  y=gadd(ghalf,gdivgs(gmul(z,y),-12));
  y=gadd(y,ginv(u1));
  toadd=(long)ceil(9.065*gtodouble(gimag(z)));
  av1=avma; lim=stack_lim(av1,1); qn=q;
  do
  {
    p1=gadd(gdiv(u,gsub(gmul(qn,u),gun)),ginv(gsub(u,qn)));
    p1=gmul(qn,p1);
    y=gadd(y,p1);
    qn=gmul(q,qn);
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
      if(DEBUGMEM>1) err(warnmem,"ellzeta");
      gerepilemany(av1,gptr,2);
    }
  }
  while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);

  y=gmul(gdiv(pii2,om2),y);
  tetpil=avma;
  return gerepile(av,tetpil,gadd(y,et));
}

/* if flag=0, return ellsigma, otherwise return log(ellsigma) */

static GEN
ellsigmasum(GEN om, GEN z, long flag, long prec)
{
  long av=avma,tetpil,lim,av1,si,toadd,n;
  GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,qn,qn2,urn,urninv,v,ga;
  GEN uinv,x1,x2,et,etnew,uhalf,q8;

  if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmasum");
  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  tau=gdiv(om1,om2); si=gsigne(gimag(tau));
  if (si==0)
    err(talker,"omega1 and omega2 are R-linearly dependent in ellsigma");
  if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
  v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
  om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
  om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
  z=gdiv(z,om2);

  x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
  x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
  et=elleta(om,prec);
  etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
  etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit));
  if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1));
  if (gexpo(z) < 5 - bit_accuracy(prec))
  {
    if (flag)
    {
      y=glog(zinit,prec);
      tetpil=avma;
      return gerepile(av,tetpil,gadd(etnew,y));
    }
    else
    {
      et=gexp(et,prec);
      tetpil=avma;
      return gerepile(av,tetpil,gmul(etnew,zinit));
    }
  }

  y1=gadd(etnew,gmul2n(gmul(gmul(z,zinit),(GEN)et[2]),-1));

  q8=gexp(gmul2n(gmul(pii2,tau),-3),prec);
  q=gpuigs(q8,8);
  uhalf=gexp(gmul(gmul2n(pii2,-1),z),prec);
  u=gneg_i(gsqr(uhalf)); uinv=ginv(u);
  y=gzero;
  toadd=(long)ceil(9.065*gtodouble(gabs(gimag(z),prec)));
/* 9.065 = 2*Pi/log(2) */
  av1=avma; lim=stack_lim(av1,1); qn=q; qn2=gun;
  urn=uhalf; urninv=ginv(uhalf); n=0;
  do
  {
    y=gadd(y,gmul(qn2,gsub(urn,urninv)));
    qn2=gmul(qn,qn2);
    qn=gmul(q,qn);
    urn=gmul(urn,u); urninv=gmul(urninv,uinv);
    n++;
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[5]; gptr[0]=&y; gptr[1]=&qn; gptr[2]=&qn2; gptr[3]=&urn;
      gptr[4]=&urninv;
      if(DEBUGMEM>1) err(warnmem,"ellsigma");
      gerepilemany(av1,gptr,5);
    }
  }
  while (gexpo(qn2) + (n-1)*toadd > - bit_accuracy(prec) - 5);

  p1=gmul(q8,gmul(gdiv(gdiv((GEN)om[2],pii2),gpuigs(trueeta(tau,prec),3)),y));
  if (flag)
  {
    p1=glog(p1,prec); tetpil=avma;
    return gerepile(av,tetpil,gadd(y1,p1));
  }
  else
  {
    y=gexp(y1,prec); tetpil=avma;
    return gerepile(av,tetpil,gmul(p1,y));
  }
}

/* if flag=0, return ellsigma, otherwise return log(ellsigma) */

static GEN
ellsigmaprod(GEN om, GEN z, long flag, long prec)
{
  long av=avma,tetpil,lim,av1,si,toadd;
  GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,u1,qn,v,ga,negu,uinv,x1,x2,et,etnew,uhalf;

  if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmaprod");
  p1=mppi(prec); setexpo(p1,2);
  pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
  tau=gdiv(om1,om2); si=gsigne(gimag(tau));
  if (si==0)
    err(talker,"omega1 and omega2 are R-linearly dependent in ellsigma");
  if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
  v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
  om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
  om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
  z=gdiv(z,om2);

  x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
  x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
  et=elleta(om,prec);
  etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
  etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit));
  if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1));
  if (gexpo(z) < 5 - bit_accuracy(prec))
  {
    if (flag)
    {
      y=glog(zinit,prec);
      tetpil=avma;
      return gerepile(av,tetpil,gadd(etnew,y));
    }
    else
    {
      et=gexp(et,prec);
      tetpil=avma;
      return gerepile(av,tetpil,gmul(etnew,zinit));
    }
  }

  y1=gadd(etnew,gmul2n(gmul(gmul(z,zinit),(GEN)et[2]),-1));

  q=gexp(gmul(pii2,tau),prec);
  uhalf=gexp(gmul(gmul2n(pii2,-1),z),prec); u=gsqr(uhalf);
  uinv=ginv(u);
  u1=gsub(uhalf,ginv(uhalf));
  y=gdiv(gmul(om2,u1),pii2);
  toadd=(long)ceil(9.065*gtodouble(gabs(gimag(z),prec)));
/* 9.065 = 2*Pi/log(2) */
  av1=avma; lim=stack_lim(av1,1); qn=q;
  negu=stoi(-1);
  do
  {
    p1=gmul(gadd(gmul(qn,u),negu),gadd(gmul(qn,uinv),negu));
    p1=gdiv(p1,gsqr(gadd(qn,negu)));
    y=gmul(y,p1);
    qn=gmul(q,qn);
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
      if(DEBUGMEM>1) err(warnmem,"ellsigma");
      gerepilemany(av1,gptr,2);
    }
  }
  while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);

  if (flag)
  {
    p1=glog(y,prec);
    tetpil=avma;
    return gerepile(av,tetpil,gadd(y1,p1));
  }
  else
  {
    p1=gexp(y1,prec);
    tetpil=avma;
    return gerepile(av,tetpil,gmul(p1,y));
  }
}

GEN
ellsigma(GEN om, GEN z, long flag, long prec)
{
  if (flag>=2) return ellsigmaprod(om,z,flag&1,prec);
  else return ellsigmasum(om,z,flag,prec);
}

GEN
pointell(GEN e, GEN z, long prec)
{
  long av=avma,tetpil;
  GEN y,yp,v,p1;

  checkbell(e);
  p1=weipellnumall((GEN)e[16],(GEN)e[15],z,1,prec);
  if (lg(p1)==2) { avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v; }
  y = gsub((GEN)p1[1], gdivgs((GEN)e[6],12));
  yp = gsub((GEN)p1[2], gmul2n(ellLHS0(e,y),-1));
  tetpil=avma; v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lcopy(yp);
  return gerepile(av,tetpil,v);
}

GEN
weipell(GEN e, long prec)
{
  long av1,tetpil,precres,i,k,l;
  GEN res,p1,s,t;

  checkell(e); precres = 2*prec+2;
  res=cgetg(precres,t_SER);
  res[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0);
  if (!prec) { setsigne(res,0); return res; }
  for (i=3; i<precres; i+=2) res[i]=zero;
  switch(prec)
  {
    default: res[8]=ldivgs((GEN)e[11],6048);
    case 3: res[6]=ldivgs((GEN)e[10],240);
    case 2: res[4]=zero;
    case 1: res[2]=un;
    case 0: break;
  }
  for (k=4; k<prec; k++)
  {
    av1 = avma;
    s = k&1? gzero: gsqr((GEN)res[k+2]);
    t = gzero;
    for (l=2; l+l<k; l++)
      t = gadd(t, gmul((GEN)res[(l+1)<<1],(GEN)res[(k-l+1)<<1]));
    p1=gmulsg(3,gadd(s,gmul2n(t,1)));
    tetpil=avma;
    p1=gdivgs(p1,(k-3)*(2*k+1));
    res[(k+1)<<1] = lpile(av1,tetpil,p1);
  }
  return res;
}

GEN
ellwp0(GEN om, GEN z, long flag, long prec, long PREC)
{
  GEN v,om1,om2;
  long av = avma;

  if (z==NULL) return weipell(om,PREC);
  if (typ(z)==t_POL)
  {
    if (lgef(z) != 4 || !gcmp0((GEN)z[2]) || !gcmp1((GEN)z[3]))
      err(talker,"expecting a simple variable in ellwp");
    v = weipell(om,PREC); setvarn(v, varn(z));
    return v;
  }
  if (!get_periods(om, &om1, &om2)) err(typeer,"ellwp");
  switch(flag)
  {
    case 0: v=weipellnumall(om1,om2,z,0,prec);
      if (typ(v)==t_VEC && lg(v)==2) { avma=av; v=gpuigs(z,-2); }
      return v;
    case 1: v=weipellnumall(om1,om2,z,1,prec);
      if (typ(v)==t_VEC && lg(v)==2)
      {
        GEN p1 = gmul2n(gpuigs(z,3),1);
        long tetpil=avma;
        v=cgetg(3,t_VEC);
	v[1]=lpuigs(z,-2);
	v[2]=lneg(p1); return gerepile(av,tetpil,v);
      }
      return v;
    case 2: return pointell(om,z,prec);
    default: err(flagerr,"ellwp"); return NULL;
  }
}

/* compute a_2 using Jacobi sum */
static GEN
_a_2(GEN e)
{
  long av = avma;
  GEN unmodp = gmodulss(1,8);
  ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]);
  ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]);
  ulong e72= itos((GEN)gmul(unmodp,gmul2n((GEN)e[7],1))[2]);
  long s = kross(e8, 2) + kross(e8 + e72 + e6 + 4, 2);
  avma = av; return stoi(-s);
}

/* a_p using Jacobi sums */
static GEN
apell2_intern(GEN e, ulong p)
{
  if (p == 2) return _a_2(e);
  else
  {
    long av=avma,i;
    GEN unmodp = gmodulss(1,p);
    ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]);
    ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]);
    ulong e72= itos((GEN)gmul(unmodp,(GEN)e[7])[2]) << 1;
    long s = kross(e8, p);

    if (p < 757)
      for (i=1; i<p; i++)
        s += kross(e8 + i*(e72 + i*(e6 + (i<<2))), p);
    else
      for (i=1; i<p; i++)
        s += kross(e8 + mulssmod(i, e72 + mulssmod(i, e6 + (i<<2), p), p), p);
    avma=av; return stoi(-s);
  }
}

GEN
apell2(GEN e, GEN pp)
{
  checkell(e); if (typ(pp)!=t_INT) err(elliper1);
  if (expi(pp) > 29)
    err(talker,"prime too large in jacobi apell2, use apell instead");

  return apell2_intern(e, (ulong)pp[2]);
}

GEN ellap0(GEN e, GEN p, long flag)
{
  return flag? apell2(e,p): apell(e,p);
}

/* invert all elements of x mod p using Montgomery's trick */
GEN
multi_invmod(GEN x, GEN p)
{
  long i, lx = lg(x);
  GEN u,y = cgetg(lx, t_VEC);

  y[1] = x[1];
  for (i=2; i<lx; i++)
    y[i] = lresii(mulii((GEN)y[i-1], (GEN)x[i]), p);

  u = mpinvmod((GEN)y[--i], p);
  for ( ; i > 1; i--)
  {
    y[i] = lresii(mulii(u, (GEN)y[i-1]), p);
    u = resii(mulii(u, (GEN)x[i]), p); /* u = 1 / (x[1] ... x[i-1]) */
  }
  y[1] = (long)u; return y;
}

static GEN
addsell(GEN e, GEN z1, GEN z2, GEN p)
{
  GEN p1,p2,x,x1,x2,y,y1,y2;
  long av = avma;

  if (!z1) return z2;
  if (!z2) return z1;

  x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
  x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
  p2 = subii(x2, x1);
  if (p2 == gzero)
  {
    if (!signe(y1) || !egalii(y1,y2)) return NULL;
    p2 = shifti(y1,1);
    p1 = addii(e, mulii(x1,mulsi(3,x1)));
    p1 = resii(p1, p);
  }
  else p1 = subii(y2,y1);
  p1 = mulii(p1, mpinvmod(p2, p));
  p1 = resii(p1, p);
  x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p);
  y = negi(addii(y1, mulii(p1,subii(x,x1))));
  avma = av; p1 = cgetg(3,t_VEC);
  p1[1] = licopy(x);
  p1[2] = lmodii(y, p); return p1;
}

/* z1 <-- z1 + z2 */
static void
addsell_part2(GEN e, GEN z1, GEN z2, GEN p, GEN p2inv)
{
  GEN p1,x,x1,x2,y,y1,y2;

  x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
  x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
  if (x1 == x2)
  {
    p1 = addii(e, mulii(x1,mulsi(3,x1)));
    p1 = resii(p1, p);
  }
  else p1 = subii(y2,y1);

  p1 = mulii(p1, p2inv);
  p1 = resii(p1, p);
  x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p);
  y = negi(addii(y1, mulii(p1,subii(x,x1)))); y = modii(y,p);
  affii(x, x1);
  affii(y, y1);
}

static GEN
powsell(GEN e, GEN z, GEN n, GEN p)
{
  GEN y;
  long s=signe(n),i,j;
  ulong m;

  if (!s || !z) return NULL;
  if (s < 0)
  {
    n = negi(n); y = cgetg(3,t_VEC);
    y[2] = lnegi((GEN)z[2]);
    y[1] = z[1]; z = y;
  }
  if (is_pm1(n)) return z;

  y = NULL;
  for (i=lgefint(n)-1; i>2; i--)
    for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
    {
      if (m&1) y = addsell(e,y,z,p);
      z = addsell(e,z,z,p);
    }
  for (m=n[2]; m>1; m>>=1)
  {
    if (m&1) y = addsell(e,y,z,p);
    z = addsell(e,z,z,p);
  }
  return addsell(e,y,z,p);
}

/* make sure *x has lgefint >= k */
static void
_fix(GEN x, long k)
{
  GEN y = (GEN)*x;
  if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
}

/* low word of integer x */
#define _low(x) ((x)[lgefint(x)-1])

/* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 20, say */
GEN
apell1(GEN e, GEN p)
{
  long *tx, *ty, *ti, av = avma, av2,pfinal,i,j,j2,s,flc,flcc,x,nb;
  GEN p1,p2,p3,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,acon,bcon,c4,c6,cp4,pts;

  if (DEBUGLEVEL) timer2();
  p1 = gmodulsg(1,p);
  c4 = gdivgs(gmul(p1,(GEN)e[10]), -48);
  c6 = gdivgs(gmul(p1,(GEN)e[11]), -864);
  pordmin = gceil(gmul2n(gsqrt(p,DEFAULTPREC),2));
  p1p = addsi(1,p); p2p = shifti(p1p,1);
  x=0; flcc=0; flc = kronecker((GEN)c6[2],p);
  u=c6; acon=gzero; bcon=gun; h=p1p;
  for(;;)
  {
    while (flc==flcc || !flc)
    {
      x++;
      u = gadd(c6, gmulsg(x, gaddgs(c4,x*x)));
      flc = kronecker((GEN)u[2],p);
    }
    flcc = flc;
    f = cgetg(3,t_VEC);
    f[1] = (long)lift_intern(gmulsg(x,u));
    f[2] = (long)lift_intern(gsqr(u));
    cp4 = lift_intern(gmul(c4, (GEN)f[2]));
    fh = powsell(cp4,f,h,p);
    s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1;
    nb = min(128, s >> 1);
    if (bcon == gun)
    { /* first time: initialize */
      tx = newbloc(s+1); *tx = evaltyp(t_VECSMALL) | evallg(s+1);
      ty = newbloc(s+1);
      ti = newbloc(s+1);
    }
    else f = powsell(cp4,f,bcon,p);
    if (!fh) goto FOUND;

    p1 = gcopy(fh);
    pts = new_chunk(nb+1);
    j = lgefint(p);
    for (i=1; i<=nb; i++)
    { /* baby steps */
      pts[i] = (long)p1;
      _fix(p1+1, j); tx[i] = _low((GEN)p1[1]);
      _fix(p1+2, j); ty[i] = _low((GEN)p1[2]);
      p1 = addsell(cp4,p1,f,p); /* f^h * F^nb */
      if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; }
    }
    mfh = dummycopy(fh);
    mfh[2] = lnegi((GEN)mfh[2]);
    fg = addsell(cp4,p1,mfh,p); /* F^nb */
    if (!fg) { h = addii(h, mulsi(nb,bcon)); goto FOUND; }
    u = cgetg(nb+1, t_VEC);
    av2 = avma; /* more baby steps, nb points at a time */
    while (i <= s)
    {
      long maxj;
      for (j=1; j<=nb; j++)
      {
        p1 = (GEN)pts[j];
        u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
        if (u[j] == zero)
        {
          if (!signe(p1[2]) || !egalii((GEN)p1[2],(GEN)fg[2]))
            { h = addii(h, mulsi(i+j-1,bcon)); goto FOUND; }
          /* doubling never occurs */
          err(bugparier, "apell1: doubling?");
        }
      }
      v = multi_invmod(u, p);
      maxj = (i-1 + nb <= s)? nb: s % nb;
      for (j=1; j<=maxj; j++,i++)
      {
        p1 = (GEN)pts[j];
        addsell_part2(cp4,p1,fg,p, (GEN)v[j]);
        tx[i] = _low((GEN)p1[1]);
        ty[i] = _low((GEN)p1[2]);
      }
      avma = av2;
    }
    p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = f^(s-1) */
    if (DEBUGLEVEL) msgtimer("[apell1] baby steps, s = %ld",s);

    /* giant steps: fg = f^s */
    fg = addsell(cp4,p1,f,p);
    if (!fg) { h = addii(h, mulsi(s,bcon)); goto FOUND; }
    pfinal = _low(p); av2 = avma;

    p1 = gen_sort(tx, cmp_IND | cmp_C, NULL);
    for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
    for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
    for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
    if (DEBUGLEVEL) msgtimer("[apell1] sorting");
    avma = av2;

    gaffect(fg, (GEN)pts[1]);
    for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */
      gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]);
    /* replace fg by fg^nb since we do nb points at a time */
    avma = av2;
    fg = gcopy((GEN)pts[nb]);
    av2 = avma;

    for (i=1,j=1; ; i++)
    {
      GEN ftest = (GEN)pts[j];
      ulong m, l = 1, r = s+1;
      long k, k2;

      avma = av2;
      k = _low((GEN)ftest[1]);
      while (l<r)
      {
        m = (l+r) >> 1;
        if (tx[m] < k) l = m+1; else r = m;
      }
      if (r <= s && tx[r] == k)
      {
        while (tx[r] == k && r) r--;
        k2 = _low((GEN)ftest[2]);
        for (r++; tx[r] == k && r <= s; r++)
          if (ty[r] == k2 || ty[r] == pfinal - k2)
          { /* [h+j2] f == ± ftest (= [i.s] f)? */
            if (DEBUGLEVEL) msgtimer("[apell1] giant steps, i = %ld",i);
            j2 = ti[r] - 1;
            p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p);
            if (egalii((GEN)p1[1], (GEN)ftest[1]))
            {
              h  = addii(h, mulsi(j2,bcon));
              if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i;
              h = addii(h, mulsi(s, mulsi(i, bcon)));
              goto FOUND;
            }
          }
      }
      if (++j > nb)
      { /* compute next nb points */
        long save;
        for (j=1; j<=nb; j++)
        {
          p1 = (GEN)pts[j];
          u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
          if (u[j] == zero) /* occurs once: i = j = nb, p1 == fg */
          {
            u[j] = lshifti((GEN)p1[2],1);
            save = fg[1]; fg[1] = p1[1];
          }
        }
        v = multi_invmod(u, p);
        for (j=1; j<=nb; j++)
          addsell_part2(cp4, (GEN)pts[j],fg,p, (GEN)v[j]);
        if (i == nb) { fg[1] = save; }
        j = 1;
      }
    }

FOUND: /* success, found a point of exponent h */
    p2 = decomp(h); p1=(GEN)p2[1]; p2=(GEN)p2[2];
    for (i=1; i<lg(p1); i++)
      for (j=itos((GEN)p2[i]); j; j--)
      {
        p3 = divii(h,(GEN)p1[i]);
        if (powsell(cp4,f,p3,p)) break;
	h = p3;
      }
    /* now h is the exact order */
    if (bcon == gun) bcon = h;
    else
    {
      p1 = chinois(gmodulcp(acon,bcon), gmodulsg(0,h));
      acon = (GEN)p1[2];
      bcon = (GEN)p1[1];
    }

    i = (cmpii(bcon, pordmin) < 0);
    if (i) acon = centermod(subii(p2p,acon), bcon);
    p1 = ground(gdiv(gsub(p1p,acon),bcon));
    h = addii(acon, mulii(p1,bcon));
    if (!i) break;
  }
  gunclone(tx);
  gunclone(ty);
  gunclone(ti);
  p1 = (flc==1)? subii(p1p,h): subii(h,p1p);
  return gerepileupto(av,p1);
}

typedef struct
{
  int isnull;
  long x,y;
} sellpt;

/* set p1 <-- p1 + p2, safe with p1 = p2 */
static void
addsell1(long e, long p, sellpt *p1, sellpt *p2)
{
  long num, den, lambda;

  if (p1->isnull) { *p1 = *p2; return; }
  if (p2->isnull) return;
  if (p1->x == p2->x)
  {
    if (! p1->y || p1->y != p2->y) { p1->isnull = 1; return; }
    num = addssmod(e, mulssmod(3, mulssmod(p1->x, p1->x, p), p), p);
    den = addssmod(p1->y, p1->y, p);
  }
  else
  {
    num = subssmod(p1->y, p2->y, p);
    den = subssmod(p1->x, p2->x, p);
  }
  lambda = divssmod(num, den, p);
  num = subssmod(mulssmod(lambda, lambda, p), addssmod(p1->x, p2->x, p), p);
  p1->y = subssmod(mulssmod(lambda, subssmod(p2->x, num, p), p), p2->y, p);
  p1->x = num; /* necessary in case p1 = p2: we need p2->x above */
}

static void
powssell1(long e, long p, long n, sellpt *p1, sellpt *p2)
{
  sellpt p3 = *p1;

  if (n < 0) { n = -n; if (p3.y) p3.y = p - p3.y; }
  p2->isnull = 1;
  for(;;)
  {
    if (n&1) addsell1(e, p, p2, &p3);
    n>>=1; if (!n) return;
    addsell1(e, p, &p3, &p3);
  }
}

typedef struct
{
  long x,y,i;
} multiple;

static int
compare_multiples(multiple *a, multiple *b)
{
  return a->x - b->x;
}

/* assume e has good reduction at p. Should use Montgomery. */
static GEN
apell0(GEN e, long p)
{
  GEN p1,p2;
  sellpt f,fh,fg,ftest,f2;
  long pordmin,u,p1p,p2p,acon,bcon,c4,c6,cp4;
  long av,i,j,s,flc,flcc,x,q,h,p3,l,r,m;
  multiple *table;

  if (p < 99) return apell2_intern(e,(ulong)p);

  av = avma; p1 = gmodulss(1,p);
  c4 = itos((GEN)gdivgs(gmul(p1,(GEN)e[10]), -48)[2]);
  c6 = itos((GEN)gdivgs(gmul(p1,(GEN)e[11]), -864)[2]);
  pordmin = (long)(1 + 4*sqrt((float)p));
  p1p = p+1; p2p = p1p << 1;
  x=0; flcc=0; flc = kross(c6, p);
  u=c6; acon=0; bcon=1; h=p1p;
  for(;;)
  {
    while (flc==flcc || !flc)
    {
      x++;
      u = addssmod(c6, mulssmod(x, c4+mulssmod(x,x,p), p), p);
      flc = kross(u,p);
    }
    flcc = flc;
    f.isnull = 0;
    f.x = mulssmod(x, u, p);
    f.y = mulssmod(u, u, p);
    cp4 = mulssmod(c4, f.y, p);
    powssell1(cp4, p, h, &f, &fh);
    s = (long) (sqrt(((float)pordmin)/bcon) / 2);
    if (!s) s=1;
    if (bcon==1)
    {
      table = (multiple *) gpmalloc((s+1)*sizeof(multiple));
      f2 = f;
    }
    else powssell1(cp4, p, bcon, &f, &f2);
    for (i=0; i < s; i++)
    {
      if (fh.isnull) { h += bcon*i; goto FOUND; }
      table[i].x = fh.x;
      table[i].y = fh.y;
      table[i].i = i;
      addsell1(cp4, p, &fh, &f2);
    }
    qsort(table,s,sizeof(multiple),(QSCOMP)compare_multiples);
    powssell1(cp4, p, s, &f2, &fg); ftest = fg;
    for (i=1; ; i++)
    {
      if (ftest.isnull) err(bugparier,"apell (f^(i*s) = 1)");
      l=0; r=s;
      while (l<r)
      {
	m = (l+r) >> 1;
	if (table[m].x < ftest.x) l=m+1; else r=m;
      }
      if (r < s && table[r].x == ftest.x) break;
      addsell1(cp4, p, &ftest, &fg);
    }
    h += table[r].i * bcon;
    if (table[r].y == ftest.y) i = -i;
    h += s * i * bcon;

FOUND:
    p2=decomp(stoi(h)); p1=(GEN)p2[1]; p2=(GEN)p2[2];
    for (i=1; i < lg(p1); i++)
      for (j = mael(p2,i,2); j; j--)
      {
	p3 = h / mael(p1,i,2);
	powssell1(cp4, p, p3, &f, &fh);
	if (!fh.isnull) break;
	h = p3;
      }
    if (bcon == 1) bcon = h;
    else
    {
      p1 = chinois(gmodulss(acon,bcon), gmodulss(0,h));
      acon = itos((GEN)p1[2]);
      if (is_bigint(p1[1])) { h = acon; break; }
      bcon = itos((GEN)p1[1]);
    }

    i = (bcon < pordmin);
    if (i)
    {
      acon = (p2p - acon) % bcon;
      if ((acon << 1) > bcon) acon -= bcon;
    }
    q = ((ulong)(p2p + bcon - (acon << 1))) / (bcon << 1);
    h = acon + q*bcon;
    avma = av; if (!i) break;
  }
  free(table); return stoi((flc==1)? p1p-h: h-p1p);
}

GEN
apell(GEN e, GEN p)
{
  checkell(e);
  if (typ(p)!=t_INT || signe(p)<0) err(talker,"not a prime in apell");
  if (gdivise((GEN)e[12],p)) /* e[12] may be an intmod */
  {
    long av = avma,s;
    GEN c6 = gmul((GEN)e[11],gmodulsg(1,p));
    s = kronecker((GEN)c6[2],p); avma=av;
    switch(mod4(p))
    {
      case 0:
      case 3: s = -s;
    }
    return stoi(s);
  }
  if (cmpis(p, 0x3fffffff) > 0) return apell1(e, p);
  return apell0(e, itos(p));
}

/* TEMPC is the largest prime whose square is less than HIGHBIT */
#ifndef LONG_IS_64BIT
#  define TEMPC 46337
#  define TEMPMAX 16777215UL
#else
#  define TEMPC 3037000493
#  define TEMPMAX 4294967295UL
#endif

GEN
anell(GEN e, long n)
{
  long tab[4]={0,1,1,-1}; /* p prime; (-1/p) = tab[p&3]. tab[0] is not used */
  long p, i, m, av, tetpil;
  GEN p1,p2,an;

  checkell(e);
  if (n <= 0) return cgetg(1,t_VEC);
  if (n>TEMPMAX) err(impl,"anell for n>=2^24 (or 2^32 for 64 bit machines)");
  an = cgetg(n+1,t_VEC); an[1] = un;
  for (i=2; i <= n; i++) an[i] = 0;
  for (p=2; p<=n; p++)
    if (!an[p])
    {
      if (!smodis((GEN)e[12],p)) /* mauvaise reduction, p | e[12] */
	switch (tab[p&3] * krogs((GEN)e[11],p)) /* renvoie (-c6/p) */
	{
	  case -1:  /* non deployee */
	    for (m=p; m<=n; m+=p)
	      if (an[m/p]) an[m]=lneg((GEN)an[m/p]);
	    continue;
	  case 0:   /* additive */
	    for (m=p; m<=n; m+=p) an[m]=zero;
	    continue;
	  case 1:   /* deployee */
	    for (m=p; m<=n; m+=p)
	      if (an[m/p]) an[m]=lcopy((GEN)an[m/p]);
	}
      else /* bonne reduction */
      {
        GEN ap = apell0(e,p);

	if (p < TEMPC)
	{
	  ulong pk, oldpk = 1;
	  for (pk=p; pk <= n; oldpk=pk, pk *= p)
	  {
	    if (pk == p) an[pk] = (long) ap;
	    else
	    {
	      av = avma;
	      p1 = mulii(ap, (GEN)an[oldpk]);
	      p2 = mulsi(p, (GEN)an[oldpk/p]);
	      tetpil = avma;
	      an[pk] = lpile(av,tetpil,subii(p1,p2));
	    }
	    for (m = n/pk; m > 1; m--)
	      if (an[m] && m%p) an[m*pk] = lmulii((GEN)an[m], (GEN)an[pk]);
	  }
	}
	else
	{
	  an[p] = (long) ap;
	  for (m = n/p; m > 1; m--)
	    if (an[m] && m%p) an[m*p] = lmulii((GEN)an[m], (GEN)an[p]);
	}
      }
    }
  return an;
}

GEN
akell(GEN e, GEN n)
{
  long i,j,ex,av=avma;
  GEN p1,p2,ap,u,v,w,y,pl;

  checkell(e);
  if (typ(n)!=t_INT) err(talker,"not an integer type in akell");
  if (signe(n)<= 0) return gzero;
  y=gun; if (gcmp1(n)) return y;
  p2=auxdecomp(n,1); p1=(GEN)p2[1]; p2=(GEN)p2[2];
  for (i=1; i<lg(p1); i++)
  {
    pl=(GEN)p1[i];
    if (divise((GEN)e[12], pl)) /* mauvaise reduction */
    {
      j = (((mod4(pl)+1)&2)-1)*kronecker((GEN)e[11],pl);
      if (j<0 && mpodd((GEN)p2[i])) y = negi(y);
      if (!j) { avma=av; return gzero; }
    }
    else /* bonne reduction */
    {
      ap=apell(e,pl); ex=itos((GEN)p2[i]);
      u=ap; v=gun;
      for (j=2; j<=ex; j++)
      {
	w = subii(mulii(ap,u), mulii(pl,v));
	v=u; u=w;
      }
      y = mulii(u,y);
    }
  }
  return gerepileupto(av,y);
}

GEN
hell(GEN e, GEN a, long prec)
{
  long av=avma,tetpil,n;
  GEN p1,p2,y,z,q,pi2surw,pi2isurw,qn,ps;

  checkbell(e);
  pi2surw=gdiv(gmul2n(mppi(prec),1),(GEN)e[15]);
  pi2isurw=cgetg(3,t_COMPLEX); pi2isurw[1]=zero; pi2isurw[2]=(long)pi2surw;
  z=gmul(greal(zell(e,a,prec)),pi2surw);
  q=greal(gexp(gmul((GEN)e[16],pi2isurw),prec));
  y=gsin(z,prec); n=0; qn=gun; ps=gneg_i(q);
  do
  {
    n++; p1=gsin(gmulsg(2*n+1,z),prec); qn=gmul(qn,ps);
    ps=gmul(ps,q); p1=gmul(p1,qn); y=gadd(y,p1);
  }
  while (gexpo(qn) >= - bit_accuracy(prec));
  p1=gmul(gsqr(gdiv(gmul2n(y,1), d_ellLHS(e,a))),pi2surw);
  p2=gsqr(gsqr(gdiv(p1,gsqr(gsqr(denom((GEN)a[1]))))));
  p1=gdiv(gmul(p2,q),(GEN)e[12]);
  p1=gmul2n(glog(gabs(p1,prec),prec),-5);
  tetpil=avma; return gerepile(av,tetpil,gneg(p1));
}

static GEN
hells(GEN e, GEN x, long prec)
{
  GEN w,z,t,mu,e72,e82;
  long n,lim;

  t = gdiv(realun(prec),(GEN)x[1]);
  mu = gmul2n(glog(numer((GEN)x[1]),prec),-1);
  e72 = gmul2n((GEN)e[7],1);
  e82 = gmul2n((GEN)e[8],1);
  lim = 6 + (bit_accuracy(prec) >> 1);
  for (n=0; n<lim; n++)
  {
    w = gmul(t,gaddsg(4,gmul(t,gadd((GEN)e[6],gmul(t,gadd(e72,gmul(t,(GEN)e[8])))))));
    z = gsub(gun,gmul(gsqr(t),gadd((GEN)e[7],gmul(t,gadd(e82,gmul(t,(GEN)e[9]))))));
    mu = gadd(mu,gmul2n(glog(z,prec), -((n<<1)+3)));
    t = gdiv(w,z);
  }
  return mu;
}

GEN
hell2(GEN e, GEN x, long prec)
{
  GEN ep,e3,ro,p1,p2,mu,d,xp;
  long av=avma,tetpil,lx,lc,i,j,tx;

  if (!oncurve(e,x)) err(heller1);
  d=(GEN)e[12]; ro=(GEN)e[14]; e3=(gsigne(d) < 0)?(GEN)ro[1]:(GEN)ro[3];
  p1=cgetg(5,t_VEC); p1[1]=un; p1[2]=laddgs(gfloor(e3),-1);
  p1[3]=p1[4]=zero; ep=coordch(e,p1); xp=pointch(x,p1);
  tx=typ(x[1]); lx=lg(x);
  if (!is_matvec_t(tx))
  {
    if (lx<3) { avma=av; return gzero; }
    tetpil=avma; return gerepile(av,tetpil,hells(ep,xp,prec));
  }
  tx=typ(x);
  tetpil=avma; mu=cgetg(lx,tx);
  if (tx != t_MAT)
    for (i=1; i<lx; i++) mu[i]=(long)hells(ep,(GEN)xp[i],prec);
  else
  {
    lc=lg(x[1]);
    for (i=1; i<lx; i++)
    {
      p1=cgetg(lc,t_COL); mu[i]=(long)p1; p2=(GEN)xp[i];
      for (j=1; j<lc; j++) p1[j]=(long)hells(ep,(GEN)p2[j],prec);
    }
  }
  return gerepile(av,tetpil,mu);
}

GEN
hell0(GEN e, GEN z, long prec)
{
  GEN a,b,s,x,u,v,u1,p1,p2;	
  long n,i, ex = 5-bit_accuracy(prec);

  /* cf. zell mais ne marche pas. Comment corriger? K.B. */
  x = new_coords(e,(GEN)z[1],&a,&b,prec);

  u = gmul2n(gadd(a,b), -1);
  v = gsqrt(gmul(a,b), prec); s = gun;
  for(n=0; ; n++)
  {
    p1 = gmul2n(gsub(x, gsqr(v)), -1);
    p2 = gsqr(u);
    x = gadd(p1, gsqrt(gadd(gsqr(p1), gmul(x, p2)), prec));
    p2 = gadd(x, p2);
    for (i=1; i<=n; i++) p2 = gsqr(p2);
    s = gmul(s, p2);
    u1 = gmul2n(gadd(u,v), -1);
    if (gexpo(gsub(u,u1)) < ex) break;

    v = gsqrt(gmul(u,v), prec);
    u = u1;
  }
  return gmul2n(glog(gdiv(gsqr(p2), s), prec) ,-1);
}

/* On suppose que `e' est a coeffs entiers donnee par un modele minimal */
static GEN
ghell0(GEN e, GEN a, long flag, long prec)
{
  long av=avma,lx,i,n,n2,grandn,tx=typ(a);
  GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;

  checkbell(e); if (!is_matvec_t(tx)) err(elliper1);
  lx = lg(a); if (lx==1) return cgetg(1,tx);
  tx=typ(a[1]);
  if (is_matvec_t(tx))
  {
    z=cgetg(lx,tx);
    for (i=1; i<lx; i++) z[i]=(long)ghell0(e,(GEN)a[i],flag,prec);
    return z;
  }
  if (lg(a)<3) return gzero;
  if (!oncurve(e,a)) err(heller1);

  psi2=numer(d_ellLHS(e,a));
  if (!signe(psi2)) { avma=av; return gzero; }

  x=(GEN)a[1]; y=(GEN)a[2];
  p2=gadd(gmulsg(3,(GEN)e[7]),gmul(x,gadd((GEN)e[6],gmulsg(3,x))));
  psi3=numer(gadd((GEN)e[9],gmul(x,gadd(gmulsg(3,(GEN)e[8]),gmul(x,p2)))));
  if (!signe(psi3)) { avma=av; return gzero; }

  p1 = gmul(x,gadd(shifti((GEN)e[2],1),gmulsg(3,x)));
  phi2=numer(gsub(gadd((GEN)e[4],p1), gmul((GEN)e[1],y)));
  p1=(GEN)factor(mppgcd(psi2,phi2))[1]; lx=lg(p1);
  switch(flag)
  {
    case 0: z = hell2(e,a,prec); break; /* Tate 4^n */
    case 1: z = hell(e,a,prec);  break; /* Silverman's trick */
    case 2: z = hell0(e,a,prec); break; /* Mestre's trick */

  }
  for (i=1; i<lx; i++)
  {
    p=(GEN)p1[i];
    if (signe(resii((GEN)e[10],p)))
    {
      grandn=ggval((GEN)e[12],p);
      if (grandn)
      {
        n2=ggval(psi2,p); n=n2<<1;
        logdep=gneg_i(glog(p,prec));
	if (n>grandn) n=grandn;
	p2=divrs(mulsr(n*(grandn+grandn-n),logdep),grandn<<3);
	z=gadd(z,p2);
      }
    }
    else
    {
      n2=ggval(psi2,p);
      logdep=gneg_i(glog(p,prec));
      n=ggval(psi3,p);
      if (n>=3*n2) p2=gdivgs(mulsr(n2,logdep),3);
      else p2=gmul2n(mulsr(n,logdep),-3);
      z=gadd(z,p2);
    }
  }
  return gerepileupto(av,z);
}

GEN
ellheight0(GEN e, GEN a, long flag, long prec)
{
  switch(flag)
  {
    case 0: return ghell(e,a,prec);
    case 1: return ghell2(e,a,prec);
    case 2: return ghell0(e,a,2,prec);
  }
  err(flagerr,"ellheight");
  return NULL; /* not reached */
}

GEN
ghell2(GEN e, GEN a, long prec)
{
  return ghell0(e,a,0,prec);
}

GEN
ghell(GEN e, GEN a, long prec)
{
  return ghell0(e,a,1,prec);
}

static long ellrootno_all(GEN e, GEN p, GEN* ptcond);

GEN
lseriesell(GEN e, GEN s, GEN A, long prec)
{
  long av=avma,av1,tetpil,lim,l,n,eps,flun;
  GEN z,p1,p2,cg,cg1,v,cga,cgb,s2,ns,gs,N;

  if (!A) A = gun;
  else
  {
    if (gsigne(A)<=0)
      err(talker,"cut-off point must be positive in lseriesell");
    if (gcmpgs(A,1) < 0) A = ginv(A);
  }
  flun = gcmp1(A) && gcmp1(s);
  eps = ellrootno_all(e,gun,&N);
  if (flun && eps<0) { z=cgetr(prec); affsr(0,z); return z; }
  cg1=mppi(prec); setexpo(cg1,2); cg=divrr(cg1,gsqrt(N,prec));
  cga=gmul(cg,A); cgb=gdiv(cg,A);
  l=(long)((pariC2*(prec-2) + fabs(gtodouble(s)-1.)*log(rtodbl(cga)))
            / rtodbl(cgb)+1);
  v = anell(e, min(l,TEMPMAX));
  if (!flun) { s2=gsubsg(2,s); ns=gpui(cg,gsubgs(gmul2n(s,1),2),prec); }
  z=gzero;
  if (typ(s)==t_INT)
  {
    if (signe(s)<=0) { avma=av; return gzero; }
    gs=mpfactr(itos(s)-1,prec);
  }
  else gs=ggamma(s,prec);
  av1=avma; lim=stack_lim(av1,1);
  for (n=1; n<=l; n++)
  {
    p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpui(stoi(n),s,prec));
    p2=flun? p1: gdiv(gmul(ns,incgam(s2,gmulsg(n,cgb),prec)),
                      gpui(stoi(n),s2,prec));
    if (eps<0) p2=gneg_i(p2);
    z = gadd(z,gmul(gadd(p1,p2),(n<=TEMPMAX)? (GEN)v[n]: akell(e,stoi(n))));
    if (low_stack(lim, stack_lim(av1,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"lseriesell");
      tetpil=avma; z=gerepile(av1,tetpil,gcopy(z));
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gdiv(z,gs));
}

/********************************************************************/
/**                                                                **/
/**                 Tate's algorithm e (cf Anvers IV)              **/
/**               Kodaira types, global minimal model              **/
/**                                                                **/
/********************************************************************/

/* Given an integral elliptic curve in ellinit form, and a prime p, returns the
  type of the fiber at p of the Neron model, as well as the change of variables
  in the form [f, kod, v, c].

  * The integer f is the conductor's exponent.

  * The integer kod is the Kodaira type using the following notation:
    II , III , IV  -->  2, 3, 4
    I0  -->  1
    Inu --> 4+nu for nu > 0
  A '*' negates the code (e.g I* --> -2)

  * v is a quadruple [u, r, s, t] yielding a minimal model

  * c is the Tamagawa number.

  Uses Tate's algorithm (Anvers IV). Given the remarks at the bottom of
  page 46, the "long" algorithm is used for p < 4 only. */
static void cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t);
static void cumule1(GEN *vtotal, GEN *e, GEN v2);

static GEN
localreduction_result(long av, long f, long kod, long c, GEN v)
{
  long tetpil = avma;
  GEN result = cgetg(5, t_VEC);
  result[1] = lstoi(f); result[2] = lstoi(kod);
  result[3] = lcopy(v); result[4] = lstoi(c);
  return gerepile(av,tetpil, result);
}

/* ici, p1 != 2 et p1 != 3 */
static GEN
localreduction_carac_not23(GEN e, GEN p1)
{
  long av = avma, k, f, kod, c, nuj, nudelta;
  GEN pk, p2k, a2prime, a3prime;
  GEN p2, r = gzero, s = gzero, t = gzero, v;
  GEN c4, c6, delta, unmodp, xun, tri, var, p4k, p6k;

  nudelta = ggval((GEN)e[12], p1);
  v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;
  nuj = gcmp0((GEN)e[13]) ? 0 : - ggval((GEN)e[13], p1);
  k = (nuj > 0 ? nudelta - nuj : nudelta) / 12;
  c4 = (GEN)e[10]; c6 = (GEN)e[11]; delta = (GEN)e[12];
  if (k > 0) /* modele non minimal */
  {
    pk = gpuigs(p1, k);
    if (mpodd((GEN)e[1]))
      s = shifti(subii(pk, (GEN)e[1]), -1);
    else
      s = negi(shifti((GEN)e[1], -1));
    p2k = sqri(pk);
    p4k = sqri(p2k);
    p6k = mulii(p4k, p2k);

    a2prime = subii((GEN)e[2], mulii(s, addii((GEN)e[1], s)));
    switch(smodis(a2prime, 3))
    {
      case 0: r = negi(divis(a2prime, 3)); break;
      case 1: r = divis(subii(p2k, a2prime), 3); break;
      case 2: r = negi(divis(addii(a2prime, p2k), 3)); break;
    }
    a3prime = ellLHS0_i(e,r);
    if (mpodd(a3prime))
      t = shifti(subii(mulii(pk, p2k), a3prime), -1);
    else
      t = negi(shifti(a3prime, -1));
    v[1] = (long)pk; v[2] = (long)r; v[3] = (long)s; v[4] = (long)t;
    nudelta -= 12 * k;
    c4 = divii(c4, p4k); c6 = divii(c6, p6k);
    delta = divii(delta, sqri(p6k));
  }
  if (nuj > 0) switch(nudelta - nuj)
  {
    case 0: f = 1; kod = 4+nuj; /* Inu */
      switch(kronecker(negi(c6),p1))
      {
	case  1: c = nudelta; break;
	case -1: c = 2 - (nudelta % 2); break;
	default: err(tater1);
      }
      break;
    case 6: f = 2; kod = -4-nuj; /* Inu* */
      if (nuj & 1)
	c = 3 + kronecker(divii(mulii(c6, delta),gpuigs(p1, 9+nuj)), p1);
      else
	c = 3 + kronecker(divii(delta, gpuigs(p1, 6+nuj)), p1);
      break;
    default: err(tater1);
  }
  else switch(nudelta)
  {
    case  0: f = 0; kod = 1; c = 1; break; /* I0, regulier */
    case  2: f = 2; kod = 2; c = 1; break; /* II   */
    case  3: f = 2; kod = 3; c = 2; break; /* III  */
    case  4: f = 2; kod = 4; /* IV   */
      c = 2 + kronecker(gdiv(mulis(c6, -6), sqri(p1)), p1);
      break;
    case  6: f = 2; kod = -1; /* I0*  */
      p2 = sqri(p1);
      unmodp = gmodulsg(1,p1);
      var = gmul(unmodp,polx[0]);
      tri = gsub(gsqr(var),gmul(divii(gmulsg(3, c4), p2),unmodp));
      tri = gsub(gmul(tri, var),
		 gmul(divii(gmul2n(c6,1), mulii(p2,p1)),unmodp));
      xun = gmodulcp(var,tri);
      c = lgef(ggcd((GEN)(gsub(gpui(xun,p1,0),xun))[2], tri)) - 2;
      break;
    case  8: f = 2; kod = -4; /* IV*  */
      c = 2 + kronecker(gdiv(mulis(c6,-6), gpuigs(p1,4)), p1);
      break;
    case  9: f = 2; kod = -3; c = 2; break; /* III* */
    case 10: f = 2; kod = -2; c = 1; break; /* II*  */
    default: err(tater1);
  }
  return localreduction_result(av,f,kod,c,v);
}

/* renvoie a_{ k,l } avec les notations de Tate */
static int
aux(GEN ak, int p, int l)
{
  long av = avma, pl = p, res;
  while (--l) pl *= p;
  res = smodis(divis(ak, pl), p);
  avma = av; return res;
}

static int
aux2(GEN ak, int p, GEN pl)
{
  long av = avma, res;
  res = smodis(divii(ak, pl), p);
  avma = av;
  return res;
}

/* renvoie le nombre de racines distinctes du polynome XXX + aXX + bX + c
 * modulo p s'il y a une racine multiple, elle est renvoyee dans *mult
 */
static int
numroots3(int a, int b, int c, int p, int *mult)
{
  if (p == 2)
  {
    if ((c + a * b) & 1) return 3;
    else { *mult = b; return (a + b) & 1 ? 2 : 1; }
  }
  else
  {
    if (a % 3) { *mult = a * b; return (a * b * (1 - b) + c) % 3 ? 3 : 2; }
    else { *mult = -c; return b % 3 ? 3 : 1; }
  }
}

/* idem pour aXX +bX + c */
static int
numroots2(int a, int b, int c, int p, int *mult)
{
  if (p == 2) { *mult = c; return b & 1 ? 2 : 1; }
  else { *mult = a * b; return (b * b - a * c) % 3 ? 2 : 1; }
}

/* ici, p1 = 2 ou p1 = 3 */
static GEN
localreduction_carac_23(GEN e, GEN p1)
{
  long av = avma, p, c, nu, nudelta;
  int a21, a42, a63, a32, a64, theroot, al, be, ga;
  GEN pk, p2k, pk1, p4, p6;
  GEN p2, p3, r = gzero, s = gzero, t = gzero, v;

  nudelta = ggval((GEN)e[12], p1);
  v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;

  for(;;)
  {
    if (!nudelta)
      return localreduction_result(av, 0, 1, 1, v);
	/* I0   */
    p = itos(p1);
    if (!divise((GEN)e[6], p1))
    {
      if (smodis(negi((GEN)e[11]), p == 2 ? 8 : 3) == 1)
	c = nudelta;
      else
	c = 2 - (nudelta & 1);
      return localreduction_result(av, 1, 4 + nudelta, c, v);
    }
	/* Inu  */
    if (p == 2)
    {
      r = modis((GEN)e[4], 2);
      s = modis(addii(r, (GEN)e[2]), 2);
      if (signe(r)) t = modis(addii(addii((GEN)e[4], (GEN)e[5]), s), 2);
      else t = modis((GEN)e[5], 2);
    }
    else /* p == 3 */
    {
      r = negi(modis((GEN)e[8], 3));
      s = modis((GEN)e[1], 3);
      t = modis(ellLHS0_i(e,r), 3);
    }
    cumule(&v, &e, gun, r, s, t); /* p | a1, a2, a3, a4 et a6 */
    p2 = stoi(p*p);
    if (!divise((GEN)e[5], p2))
      return localreduction_result(av, nudelta, 2, 1, v);
	/* II   */
    p3 = stoi(p*p*p);
    if (!divise((GEN)e[9], p3))
      return localreduction_result(av, nudelta - 1, 3, 2, v);
	/* III  */
    if (!divise((GEN)e[8], p3))
    {
      if (smodis((GEN)e[8], (p==2)? 32: 27) == p*p)
	c = 3;
      else
	c = 1;
      return localreduction_result(av, nudelta - 2, 4, c, v);
    }
	/* IV   */

	/* now for the last five cases... */

    if (!divise((GEN)e[5], p3))
      cumule(&v, &e, gun, gzero, gzero, p == 2? gdeux: modis((GEN)e[3], 9));
	/* p | a1, a2; p^2  | a3, a4; p^3 | a6 */
    a21 = aux((GEN)e[2], p, 1); a42 = aux((GEN)e[4], p, 2);
    a63 = aux((GEN)e[5], p, 3);
    switch (numroots3(a21, a42, a63, p, &theroot))
    {
      case 3:
	if (p == 2)
	  c = 1 + (a63 == 0) + ((a21 + a42 + a63) & 1);
	else
	  c = 1 + (a63 == 0) + (((1 + a21 + a42 + a63) % 3) == 0)
	      + (((1 - a21 + a42 - a63) % 3) == 0);
	return localreduction_result(av, nudelta - 4, -1, c, v);
	    /* I0*  */
      case 2: /* calcul de nu */
	if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero);
	    /* p | a1; p^2  | a2, a3; p^3 | a4; p^4 | a6 */
	nu = 1;
	pk = p2;
	p2k = stoi(p * p * p * p);
	for(;;)
	{
	  if (numroots2(al = 1, be = aux2((GEN)e[3], p, pk),
			ga = -aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
	    break;
	  if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk));
	  pk1 = pk; pk = mulsi(p, pk); p2k = mulsi(p, p2k);
	  nu++;
	  if (numroots2(al = a21, be = aux2((GEN)e[4], p, pk),
			ga = aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
	    break;
	  if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero);
	  p2k = mulsi(p, p2k);
	  nu++;
	}
	if (p == 2)
	  c = 4 - 2 * (ga & 1);
	else
	  c = 3 + kross(be * be - al * ga, 3);
	return localreduction_result(av, nudelta - 4 - nu, -4 - nu, c, v);
	    /* Inu* */
      case 1:
	if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero);
	    /* p | a1; p^2  | a2, a3; p^3 | a4; p^4 | a6 */
	a32 = aux((GEN)e[3], p, 2); a64 = aux((GEN)e[5], p, 4);
	if (numroots2(1, a32, -a64, p, &theroot) == 2)
	{
	  if (p == 2)
	    c = 3 - 2 * a64;
	  else
	    c = 2 + kross(a32 * a32 + a64, 3);
	  return localreduction_result(av, nudelta - 6, -4, c, v);
	}
	    /* IV*  */
	if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p));
	    /* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */
	p4 = sqri(p2);
	if (!divise((GEN)e[4], p4))
	  return localreduction_result(av, nudelta - 7, -3, 2, v);
	    /* III* */
	p6 = mulii(p4, p2);
	if (!divise((GEN)e[5], p6))
	  return localreduction_result(av, nudelta - 8, -2, 1, v);
	    /* II*  */
	cumule(&v, &e, p1, gzero, gzero, gzero); /* non minimal, on repart
						     pour un tour */
	nudelta -= 12;
    }
  }
  /* Not reached */
}

GEN
localreduction(GEN e, GEN p1)
{
  checkell(e);
  if (typ(e[12]) != t_INT)
    err(talker,"not an integral curve in localreduction");
  if (gcmpgs(p1, 3) > 0)       /* p different de 2 ou 3 */
    return localreduction_carac_not23(e,p1);
  else
    return localreduction_carac_23(e,p1);
}

#if 0
/*  Calcul de toutes les fibres non elliptiques d'une courbe sur Z.
 *  Etant donne une courbe elliptique sous forme longue e, dont les coefficients
 *  sont entiers, renvoie une matrice dont les lignes sont de la forme
 *  [p, fp, kodp, cp]. Il y a une ligne par diviseur premier du discriminant.
 */
GEN
globaltatealgo(GEN e)
{
  long k, l,av;
  GEN p1, p2, p3, p4, prims, result;

  checkell(e);
  prims = decomp((GEN)e[12]);
  l = lg(p1 = (GEN)prims[1]);
  p2 = (GEN)prims[2];
  if ((long)prims == avma) cgiv(prims);
  result = cgetg(5, t_MAT);
  result[1] = (long)p1;
  result[2] = (long)p2;
  result[3] = (long)(p3 = cgetg(l, t_COL));
  for (k = 1; k < l; k++) p3[k] = lgeti(3);
  result[4] = (long)(p4 = cgetg(l, t_COL));
  for (k = 1; k < l; k++) p4[k] = lgeti(3);
  av = avma;
  for (k = 1; k < l; k++)
  {
    GEN q = localreduction(e, (GEN)p1[k]);
    affii((GEN)q[1],(GEN)p2[k]);
    affii((GEN)q[2],(GEN)p3[k]);
    affii((GEN)q[4],(GEN)p4[k]);
    avma = av;
  }
  return result;
}
#endif

/* Algorithme de reduction d'une courbe sur Q a sa forme standard.  Etant
 * donne une courbe elliptique sous forme longue e, dont les coefficients
 * sont rationnels, renvoie son [N, [u, r, s, t], c], ou N est le conducteur
 * arithmetique de e, [u, r, s, t] est le changement de variables qui reduit
 * e a sa forme minimale globale dans laquelle a1 et a3 valent 0 ou 1, et a2
 * vaut -1, 0 ou 1 et tel que u est un rationnel positif. Enfin c est le
 * produit des nombres de Tamagawa locaux cp.
 */
GEN
globalreduction(GEN e1)
{
  long i, k, l, m, tetpil, av = avma;
  GEN p1, c = gun, prims, result, N = gun, u = gun, r, s, t;
  GEN v = cgetg(5, t_VEC), a = cgetg(7, t_VEC), e = cgetg(20, t_VEC);

  checkell(e1);
  for (i = 1; i < 5; i++) a[i] = e1[i]; a[5] = zero; a[6] = e1[5];
  prims = decomp(denom(a));
  p1 = (GEN)prims[1]; l = lg(p1);
  for (k = 1; k < l; k++)
  {
    int n = 0;
    for (i = 1; i < 7; i++)
      if (!gcmp0((GEN)a[i]))
      {
	m = i * n + ggval((GEN)a[i], (GEN)p1[k]);
	while (m < 0) { n++; m += i; }
      }
    u = gmul(u, gpuigs((GEN)p1[k], n));
  }
  v[1] = linv(u); v[2] = v[3] = v[4] = zero;
  for (i = 1; i < 14; i++) e[i] = e1[i];
  for (; i < 20; i++) e[i] = zero;
  if (!gcmp1(u)) e = coordch(e, v);
  prims = decomp((GEN)e[12]);
  l = lg(p1 = (GEN)prims[1]);
  for (k = (signe(e[12]) < 0) + 1; k < l; k++)
  {
    GEN q = localreduction(e, (GEN)p1[k]);
    GEN v1 = (GEN)q[3];
    N = mulii(N, gpui((GEN)p1[k],(GEN)q[1],0));
    c = mulii(c, (GEN)q[4]);
    if (!gcmp1((GEN)v1[1])) cumule1(&v, &e, v1);
  }
  s = gdiventgs((GEN)e[1], -2);
  r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3);
  t = gdiventgs(ellLHS0(e,r), -2);
  cumule(&v, &e, gun, r, s, t);
  tetpil = avma;
  result = cgetg(4, t_VEC); result[1] = lcopy(N); result[2] = lcopy(v);
  result[3] = lcopy(c);
  return gerepile(av, tetpil, result);
}

/* cumule les effets de plusieurs chgts de variable. On traite a part les cas
 * particuliers frequents, tels que soit u = 1, soit r' = s' = t' = 0
 */
static void
cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t)
{
  long av = avma, tetpil;
  GEN temp, v = *vtotal, v3 = cgetg(5, t_VEC);
  if (gcmp1((GEN)v[1]))
  {
    v3[1] = lcopy(u);
    v3[2] = ladd((GEN)v[2], r);
    v3[3] = ladd((GEN)v[3], s);
    av = avma;
    temp = gadd((GEN)v[4], gmul((GEN)v[3], r));
    tetpil = avma;
    v3[4] = lpile(av, tetpil, gadd(temp, t));
  }
  else if (gcmp0(r) && gcmp0(s) && gcmp0(t))
  {
    v3[1] = lmul((GEN)v[1], u);
    v3[2] = lcopy((GEN)v[2]);
    v3[3] = lcopy((GEN)v[3]);
    v3[4] = lcopy((GEN)v[4]);
  }
  else /* cas general */
  {
    v3[1] = lmul((GEN)v[1], u);
    temp = gsqr((GEN)v[1]);
    v3[2] = ladd(gmul(temp, r), (GEN)v[2]);
    v3[3] = ladd(gmul((GEN)v[1], s), (GEN)v[3]);
    v3[4] = ladd((GEN)v[4], gmul(temp, gadd(gmul((GEN)v[1], t), gmul((GEN)v[3], r))));

    tetpil = avma;
    v3 = gerepile(av, tetpil, gcopy(v3));
  }
  *vtotal = v3;
}

static void
cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t)
{
  long av = avma, tetpil;
  GEN v2 = cgetg(5, t_VEC);
  v2[1] = (long)u; v2[2] = (long)r; v2[3] = (long)s; v2[4] = (long)t;
  tetpil = avma;
  *e = gerepile(av, tetpil, coordch(*e, v2));
  cumulev(vtotal, u, r, s, t);
}

static void
cumule1(GEN *vtotal, GEN *e, GEN v2)
{
  *e = coordch(*e, v2);
  cumulev(vtotal, (GEN)v2[1], (GEN)v2[2], (GEN)v2[3], (GEN)v2[4]);
}

/********************************************************************/
/**                                                                **/
/**                   Parametrisation modulaire                    **/
/**                                                                **/
/********************************************************************/

GEN
taniyama(GEN e)
{
  GEN v,w,c,d,s1,s2,s3;
  long n,m,av=avma,tetpil;

  checkell(e); v = cgetg(precdl+3,t_SER);
  v[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0);
  v[2] = un;
  c=gtoser(anell(e,precdl+1),0); setvalp(c,1);
  d=ginv(c); c=gsqr(d);
  for (n=-3; n<=precdl-4; n++)
  {
    if (n!=2)
    {
      s3=n?gzero:(GEN)e[7];
      if (n>-3) s3=gadd(s3,gmul((GEN)e[6],(GEN)v[n+4]));
      s2=gzero;
      for (m=-2; m<=n+1; m++)
	s2 = gadd(s2,gmulsg(m*(n+m),gmul((GEN)v[m+4],(GEN)c[n-m+4])));
      s2=gmul2n(s2,-1);
      s1=gzero;
      for (m=-1; m+m<=n; m++)
      {
	if (m+m==n)
          s1=gadd(s1,gsqr((GEN)v[m+4]));
	else
          s1=gadd(s1,gmul2n(gmul((GEN)v[m+4],(GEN)v[n-m+4]),1));
      }
      v[n+6]=ldivgs(gsub(gadd(gmulsg(6,s1),s3),s2),(n+2)*(n+1)-12);
    }
    else
    {
      setlg(v,9); v[8]=(long)polx[MAXVARN];
      w=deriv(v,0); setvalp(w,-2);
      s1=gadd((GEN)e[8],gmul(v,gadd(gmul2n((GEN)e[7],1),gmul(v,gadd((GEN)e[6],gmul2n(v,2))))));
      setlg(v,precdl+3);
      s2=gsub(s1,gmul(c,gsqr(w)));
      s2=gsubst((GEN)s2[2],MAXVARN,polx[0]);
      v[n+6]=lneg(gdiv((GEN)s2[2],(GEN)s2[3]));
    }
  }
  w=gsub(gmul(polx[0],gmul(d,deriv(v,0))), ellLHS0(e,v));
  tetpil=avma; s1=cgetg(3,t_VEC); s1[1]=lcopy(v); s1[2]=lmul2n(w,-1);
  return gerepile(av,tetpil,s1);
}

/********************************************************************/
/**                                                                **/
/**                       TORSION POINTS (over Q)                  **/
/**                                                                **/
/********************************************************************/
/* assume e is defined over Q (use Mazur's theorem) */
GEN
orderell(GEN e, GEN p)
{
  GEN p1;
  long av=avma,k;

  checkell(e); checkpt(p);
  k=typ(e[13]);
  if (k!=t_INT && !is_frac_t(k))
    err(impl,"orderell for nonrational elliptic curves");
  p1=p; k=1;
  for (k=1; k<16; k++)
  {
    if (lg(p1)<3) { avma=av; return stoi(k); }
    p1 = addell(e,p1,p);
  }
  avma=av; return gzero;
}

/* one can do much better by factoring denom(D) (see ellglobalred) */
static GEN
ellintegralmodel(GEN e)
{
  GEN a = cgetg(6,t_VEC), v;
  long i;

  for (i=1; i<6; i++) a[i]=e[i];
  a = denom(a); if (gcmp1(a)) return NULL;
  v = cgetg(5,t_VEC);
  v[1]=linv(a); v[2]=v[3]=v[4]=zero; return v;
}

/* Using Lutz-Nagell */

/* p is a polynomial of degree exactly 3 with integral coefficients
 * and leading term 4. Outputs the vector of rational roots of p
 */
static GEN
ratroot(GEN p)
{
  GEN v,a,ld;
  long i,t;

  i=2; while (!signe(p[i])) i++;
  if (i==5)
    { v=cgetg(2,t_VEC); v[1]=zero; return v; }
  if (i==4)
    { v=cgetg(3,t_VEC); v[1]=zero; v[2]=ldivgs((GEN)p[4],-4); return v; }

  v=cgetg(4,t_VEC); t=1;
  if (i==3) v[t++]=zero;
  ld=divisors(gmul2n((GEN)p[i],2));
  for (i=1; i<lg(ld); i++)
  {
    a = gmul2n((GEN)ld[i],-2);
    if (!gsigne(poleval(p,a))) v[t++]=(long)a;
    a = gneg_i(a);
    if (!gsigne(poleval(p,a))) v[t++]=(long)a;
  }
  setlg(v,t); return v;
}

static int
is_new_torsion(GEN e, GEN v, GEN p, long t2) {
  GEN pk = p, pkprec = NULL;
  long k,l;

  for (k=2; k<=6; k++)
  {
    pk=addell(e,pk,p);
    if (lg(pk)==2) return 1;

    for (l=2; l<=t2; l++)
      if (gegal((GEN)pk[1],gmael(v,l,1))) return 1;

    if (pkprec && k<=5)
      if (gegal((GEN)pk[1],(GEN)pkprec[1])) return 1;
    pkprec=pk;
  }
  return 0;
}

GEN
torsellnagelllutz(GEN e)
{
  GEN d,ld,pol,p1,lr,r,v,w,w2,w3;
  long i,j,nlr,t,t2,k,k2,av=avma;

  checkell(e);
  v = ellintegralmodel(e);
  if (v) e = coordch(e,v);
  pol = RHSpol(e);
  lr=ratroot(pol); nlr=lg(lr)-1;
  r=cgetg(17,t_VEC); p1=cgetg(2,t_VEC); p1[1]=zero; r[1]=(long)p1;
  for (t=1,i=1; i<=nlr; i++)
  {
    p1=cgetg(3,t_VEC);
    p1[1] = lr[i];
    p1[2] = lmul2n(gneg(ellLHS0(e,(GEN)lr[i])), -1);
    r[++t]=(long)p1;
  }
  ld = factor(gmul2n(absi((GEN)e[12]), 4));
  p1 = (GEN)ld[2]; k = lg(p1);
  for (i=1; i<k; i++) p1[i] = lshifti((GEN)p1[i], -1);
  ld = divisors(ld);
  for (t2=t,j=1; j<lg(ld); j++)
  {
    d=(GEN)ld[j]; lr=ratroot(gsub(pol,gsqr(d)));
    for (i=1; i<lg(lr); i++)
    {
      p1 = cgetg(3,t_VEC);
      p1[1] = lr[i];
      p1[2] = lmul2n(gsub(d,ellLHS0(e,(GEN)lr[i])), -1);

      if (is_new_torsion(e,r,p1,t2))
      {
        GEN p2 = cgetg(3,t_VEC);
        p2[1] = p1[1];
        p2[2] = lsub((GEN)p1[2],d);
	r[++t]=(long)p1;
        r[++t]=(long)p2;
      }
    }
  }
  if (t==1)
  {
    avma=av; w=cgetg(4,t_VEC);
    w[1] = un;
    w[2] = lgetg(1,t_VEC);
    w[3] = lgetg(1,t_VEC);
    return w;
  }

  if (nlr<3)
  {
    w2=cgetg(2,t_VEC); w2[1]=lstoi(t);
    for (k=2; k<=t; k++)
      if (itos(orderell(e,(GEN)r[k])) == t) break;
    if (k>t) err(bugparier,"torsell (bug1)");

    w3=cgetg(2,t_VEC); w3[1]=r[k];
  }
  else
  {
    if (t&3) err(bugparier,"torsell (bug2)");
    t2 = t>>1;
    w2=cgetg(3,t_VEC); w2[1]=lstoi(t2); w2[2]=(long)gdeux;
    for (k=2; k<=t; k++)
      if (itos(orderell(e,(GEN)r[k])) == t2) break;
    if (k>t) err(bugparier,"torsell (bug3)");

    p1 = powell(e,(GEN)r[k],stoi(t>>2));
    k2 = (lg(p1)==3 && gegal((GEN)r[2],p1))? 3: 2;
    w3=cgetg(3,t_VEC); w3[1]=r[k]; w3[2]=r[k2];
  }
  if (v)
  {
    v[1] = linv((GEN)v[1]);
    w3 = pointch(w3,v);
  }
  w=cgetg(4,t_VEC);
  w[1] = lstoi(t);
  w[2] = (long)w2;
  w[3] = (long)w3;
  return gerepileupto(av, gcopy(w));
}

/* Using Doud's algorithm */

/* Input e and n, finds a bound for #Tor */
static long
torsbound(GEN e, long n)
{
  long av = avma, m, b, c, d, prime = 2;
  byteptr p = diffptr;
  GEN D = (GEN)e[12];

  b = c = m = 0; p++;
  while (m<n)
  {
    d = *p++; if (!d) err(primer1);
    prime += d;
    if (ggval(D,stoi(prime)) == 0)
    {
      b = cgcd(b, prime+1 - itos(apell0(e,prime)));
      if (b==c) m++; else {c = b; m = 0;}
      avma = av;
    }
  }
  return b;
}

static GEN
_round(GEN x, long *e)
{
  GEN y = grndtoi(x,e);
  if (*e > -5 && bit_accuracy(gprecision(x)) < gexpo(y) - 10)
    err(talker, "ellinit data not accurate enough. Increase precision");
  return y;
}

/* Input the curve, a point, and an integer n, returns a point of order n
   on the curve, or NULL if q is not rational. */
static GEN
torspnt(GEN E, GEN q, long n)
{
  GEN p = cgetg(3,t_VEC);
  long e;
  p[1] = lmul2n(_round(gmul2n((GEN)q[1],2), &e),-2);
  if (e > -5) return NULL;
  p[2] = lmul2n(_round(gmul2n((GEN)q[2],3), &e),-3);
  if (e > -5) return NULL;
  return (gcmp0(gimag(p)) && oncurve(E,p)
      && lg(powell(E,p,stoi(n))) == 2
      && itos(orderell(E,p)) == n)? greal(p): NULL;
}

static int
smaller_x(GEN p, GEN q)
{
  int s = absi_cmp(denom(p), denom(q));
  return (s<0 || (s==0 && absi_cmp(numer(p),numer(q)) < 0));
}

/* best generator in cycle of length k */
static GEN
best_in_cycle(GEN e, GEN p, long k)
{
  GEN p0 = p,q = p;
  long i;

  for (i=2; i+i<k; i++)
  {
    q = addell(e,q,p0);
    if (cgcd(i,k)==1 && smaller_x((GEN)q[1], (GEN)p[1])) p = q;
  }
  return (gsigne(d_ellLHS(e,p)) < 0)? invell(e,p): p;
}

static GEN
tors(GEN e, long k, GEN p, GEN q, GEN v)
{
  GEN p1,r;
  if (q)
  {
    long n = k>>1;
    GEN p1, best = q, np = powell(e,p,stoi(n));
    if (n % 2 && smaller_x((GEN)np[1], (GEN)best[1])) best = np;
    p1 = addell(e,q,np);
    if (smaller_x((GEN)p1[1], (GEN)best[1])) q = p1;
    else if (best == np) { p = addell(e,p,q); q = np; }
    p = best_in_cycle(e,p,k);
    if (v)
    {
      v[1] = linv((GEN)v[1]);
      p = pointch(p,v);
      q = pointch(q,v);
    }
    r = cgetg(4,t_VEC);
    r[1] = lstoi(2*k); p1 = cgetg(3,t_VEC); p1[1] = lstoi(k); p1[2] = deux;
    r[2] = (long)p1; p1 = cgetg(3,t_VEC); p1[1] = lcopy(p); p1[2] = lcopy(q);
    r[3] = (long)p1;
  }
  else
  {
    if (p)
    {
      p = best_in_cycle(e,p,k);
      if (v)
      {
        v[1] = linv((GEN)v[1]);
        p = pointch(p,v);
      }
      r = cgetg(4,t_VEC);
      r[1] = lstoi(k); p1 = cgetg(2,t_VEC); p1[1] = r[1];
      r[2] = (long)p1; p1 = cgetg(2,t_VEC); p1[1] = lcopy(p);
      r[3] = (long)p1;
    }
    else
    {
      r = cgetg(4,t_VEC);
      r[1] = un;
      r[2] = lgetg(1,t_VEC);
      r[3] = lgetg(1,t_VEC);
    }
  }
  return r;
}

GEN
torselldoud(GEN e)
{
  long b,i,ord,av=avma,prec, k = 1;
  GEN v,w1,w22,w1j,w12,p,tor1,tor2;

  checkbell(e);
  v = ellintegralmodel(e);
  if (v) e = coordch(e,v);

  prec=precision((GEN)e[15]);
  prec=max(prec,MEDDEFAULTPREC);
  b = torsbound(e,3);
  if (b==1) { avma=av; return tors(e,1,NULL,NULL, v); }
  w22 = gmul2n((GEN)e[16],-1);
  w1 = (GEN)e[15];
  if (b % 4)
  {
    for (i=10; i>1; i--)
    {
      if (b%i==0)
      {
	w1j = gdivgs(w1,i);
	p = torspnt(e,pointell(e,w1j,prec),i);
	if (i%2==0 && p==NULL)
	{
	  p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);
	  if (p==NULL)
	    p = torspnt(e,pointell(e,gadd(w22,gmul2n(w1j,1)),prec),i);
	}
      }
      else p = NULL;
      if (p) {k = i; break; }
    }
    return gerepileupto(av, tors(e,k,p,NULL, v));
  }

  ord = 0; tor2 = NULL;
  w12 = gmul2n((GEN)e[15],-1);
  if ((p = torspnt(e,pointell(e,w12,prec),2)))
  {
    tor1 = p; ord++;
  }
  if ((p = torspnt(e,pointell(e,w22,prec),2))
   || (!ord && (p = torspnt(e,pointell(e,gadd(w12,w22),prec),2))))
  {
    tor2 = p; ord += 2;
  }

  switch(ord)
  {
    case 0:
      for (i=9; i>1; i-=2)
      {
        w1j=gdivgs((GEN)e[15],i);
        p = (b%i==0)? torspnt(e,pointell(e,w1j,prec),i): NULL;
        if (p) { k = i; break; }
      }
      break;

    case 1:
      for (i=12; i>2; i-=2)
      {
        w1j=gdivgs((GEN)e[15],i);
        if (b%i==0)
        {
          p = torspnt(e,pointell(e,w1j,prec),i);
          if (i%4==0 && p==NULL)
            p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);
        }
        else p = NULL;
        if (p) { k = i; break; }
      }
      if (!p) { p = tor1; k = 2; }
      break;

    case 2:
      for (i=5; i>1; i-=2)
      {
        w1j = gdivgs((GEN)e[15],i);
        p = (b%i==0)? torspnt(e,pointell(e,gadd(w22,w1j),prec),i+i): NULL;
        if (p) { k = i+i; break; }
      }
      if (!p) { p = tor2; k = 2; }
      tor2 = NULL; break;

    case 3:
      for (i=8; i>2; i-=2)
      {
        w1j=gdivgs((GEN)e[15],i);
        p = (b%(2*i)==0)? torspnt(e,pointell(e,w1j,prec),i): NULL;
        if (p) { k = i; break; }
      }
      if (!p) { p = tor1; k = 2; }
      break;
  }
  return gerepileupto(av, tors(e,k,p,tor2, v));
}

GEN
elltors0(GEN e, long flag)
{
  switch(flag)
  {
    case 0: return torselldoud(e);
    case 1: return torsellnagelllutz(e);
    default: err(flagerr,"torsell");
  }
  return NULL; /* not reached */
}

/* par compatibilite */
GEN torsell(GEN e) {return torselldoud(e);}

/* LOCAL ROOT NUMBERS, D'APRES HALBERSTADT halberst@math.jussieu.fr */

/* ici p=2 ou 3 */
static long
neron(GEN e, GEN p, long* ptkod)
{
  long av=avma,kod,v4,v6,vd;
  GEN c4, c6, d, nv;

  nv=localreduction(e,p);
  kod=itos((GEN)nv[2]); *ptkod=kod;
  c4=(GEN)e[10]; c6=(GEN)e[11]; d=(GEN)e[12];
  v4=gcmp0(c4) ? 12 : ggval(c4,p);
  v6=gcmp0(c6) ? 12 : ggval(c6,p);
  vd=ggval(d,p);
  avma=av;
  switch(itos(p))
  {
    case 3:
      if (labs(kod)>4) return 1;
      else
      {
	switch(kod)
	{
	  case -1: case 1: return v4&1 ? 2 : 1;
	  case -3: case 3: return (2*v6>vd+3) ? 2 : 1;
	  case -4: case 2:
	    switch (vd%6)
	    {
	      case 4: return 3;
	      case 5: return 4;
	      default: return v6%3==1 ? 2 : 1;
	    }
	  default: /* kod = -2 et 4 */
	    switch (vd%6)
	    {
	      case 0: return 2;
	      case 1: return 3;
	      default: return 1;
	    }
	}
      }
    case 2:
      if (kod>4) return 1;
      else
      {
	switch(kod)
        {
	  case 1: return (v6>0) ? 2 : 1;
	  case 2:
	    if (vd==4) return 1;
	    else
	    {
	      if (vd==7) return 3;
	      else return v4==4 ? 2 : 4;
	    }
	  case 3:
	    switch(vd)
	    {
	      case 6: return 3;
	      case 8: return 4;
	      case 9: return 5;
	      default: return v4==5 ? 2 : 1;
	    }
	  case 4: return v4>4 ? 2 : 1;
	  case -1:
	    switch(vd)
	    {
	      case 9: return 2;
	      case 10: return 4;
	      default: return v4>4 ? 3 : 1;
	    }
	  case -2:
	    switch(vd)
	    {
	      case 12: return 2;
	      case 14: return 3;
	      default: return 1;
	    }
	  case -3:
	    switch(vd)
	    {
	      case 12: return 2;
	      case 14: return 3;
	      case 15: return 4;
	      default: return 1;
	    }
	  case -4: return v6==7 ? 2 : 1;
	  case -5: return (v6==7 || v4==6) ? 2 : 1;
	  case -6:
	    switch(vd)
	    {
	      case 12: return 2;
	      case 13: return 3;
	      default: return v4==6 ? 2 : 1;
	    }
	  case -7: return (vd==12 || v4==6) ? 2 : 1;
	  default: return v4==6 ? 2 : 1;
	}
      }
    default: return 0; /* should not occur */
  }
}

static long
ellrootno_2(GEN e)
{
  long n2,kod,u,v,x1,y1,d1,av=avma,v4,v6,w2;
  GEN p=gdeux,c4,c6,tmp,p6;

  n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p6=stoi(64);
  if (gcmp0(c4)) {v4=12; u=0;}
  else {v4=pvaluation(c4,p,&tmp); u=itos(modii(tmp,p6));}
  if (gcmp0(c6)) {v6=12; v=0;}
  else {v6=pvaluation(c6,p,&tmp); v=itos(modii(tmp,p6));}
  (void)pvaluation((GEN)e[12],p,&tmp); d1=itos(modii(tmp,p6));
  avma=av; x1=u+v+v;
  if (kod>=5)
    {w2=mpodd(addii((GEN)e[2],(GEN)e[3])) ? 1 : -1; avma=av; return w2;}
  if (kod<-9) return (n2==2) ? -kross(-1,v) : -1;
  switch(kod)
  {
    case 1: return 1;
    case 2:
      switch(n2)
      {
	case 1:
	  switch(v4)
	  {
	    case 4: return kross(-1,u);
	    case 5: return 1;
	    default: return -1;
	  }
	case 2: return (v6==7) ? 1 : -1;
	case 3: return (v%8==5 || (u*v)%8==5) ? 1 : -1;
	case 4: if (v4>5) return kross(-1,v);
	  return (v4==5) ? -kross(-1,u) : -1;
      }
    case 3:
      switch(n2)
      {
	case 1: return -kross(2,u*v);
	case 2: return -kross(2,v);
	case 3: y1=itos(modis(gsubsg(u,gmul2n(c6,-5)),16)); avma=av;
	  return (y1==7 || y1==11) ? 1 : -1;
	case 4: return (v%8==3 || (2*u+v)%8==7) ? 1 : -1;
	case 5: return v6==8 ? kross(2,x1) : kross(-2,u);
      }
    case -1:
      switch(n2)
      {
	case 1: return -kross(2,x1);
	case 2: return (v%8==7) || (x1%32==11) ? 1 : -1;
	case 3: return v4==6 ? 1 : -1;
	case 4: if (v4>6) return kross(-1,v);
	  return v4==6 ? -kross(-1,u*v) : -1;
      }
    case -2: return n2==1 ? kross(-2,v) : kross(-1,v);
    case -3:
      switch(n2)
      {
	case 1: y1=(u-2*v)%64; if (y1<0) y1+=64;
	  return (y1==3) || (y1==19) ? 1 : -1;
	case 2: return kross(2*kross(-1,u),v);
	case 3: return -kross(-1,u)*kross(-2*kross(-1,u),u*v);
	case 4: return v6==11 ? kross(-2,x1) : -kross(-2,u);
      }
    case -5:
      if (n2==1) return x1%32==23 ? 1 : -1;
      else return -kross(2,2*u+v);
    case -6:
      switch(n2)
      {
	case 1: return 1;
	case 2: return v6==10 ? 1 : -1;
	case 3: return (u%16==11) || ((u+4*v)%16==3) ? 1 : -1;
      }
    case -7:
      if (n2==1) return 1;
      else
      {
	y1=itos(modis(gaddsg(u,gmul2n(c6,-8)),16)); avma=av;
	if (v6==10) return (y1==9) || (y1==13) ? 1 : -1;
	else return (y1==9) || (y1==5) ? 1 : -1;
      }
    case -8: return n2==2 ? kross(-1,v*d1) : -1;
    case -9: return n2==2 ? -kross(-1,d1) : -1;
    default: return -1;
  }
}

static long
ellrootno_3(GEN e)
{
  long n2,kod,u,v,d1,av=avma,r6,k4,k6,v4;
  GEN p=stoi(3),c4,c6,tmp,p4;

  n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p4=stoi(81);
  if (gcmp0(c4)) { v4=12; u=0; }
  else { v4=pvaluation(c4,p,&tmp); u=itos(modii(tmp,p4)); }
  if (gcmp0(c6)) v=0;
  else {(void)pvaluation(c6,p,&tmp); v=itos(modii(tmp,p4));}
  (void)pvaluation((GEN)e[12],p,&tmp); d1=itos(modii(tmp,p4));
  avma=av;
  r6=v%9; k4=kross(u,3); k6=kross(v,3);
  if (kod>4) return k6;
  switch(kod)
  {
    case 1: case 3: case -3: return 1;
    case 2:
      switch(n2)
      {
	case 1: return (r6==4 || r6>6) ? 1 : -1;
	case 2: return -k4*k6;
	case 3: return 1;
	case 4: return -k6;
      }
    case 4:
      switch(n2)
      {
	case 1: return k6*kross(d1,3);
	case 2: return -k4;
	case 3: return -k6;
      }
    case -2: return n2==2 ? 1 : k6;
    case -4:
      switch(n2)
      {
	case 1:
	  if (v4==4) return (r6==4 || r6==8) ? 1 : -1;
	  else return (r6==1 || r6==2) ? 1 : -1;
	case 2: return -k6;
	case 3: return (r6==2 || r6==7) ? 1 : -1;
	case 4: return k6;
      }
    default: return -1;
  }
}

static long
ellrootno_not23(GEN e, GEN p, GEN ex)
{
  GEN j;
  long ep,z;

  if (gcmp1(ex)) return -kronecker(negi((GEN)e[11]),p);
  j=(GEN)e[13];
  if (!gcmp0(j) && ggval(j,p) < 0) return kronecker(negi(gun),p);
  ep=12/cgcd(12,ggval((GEN)e[12],p));
  if (ep==4) z=2;
  else z=(ep%2==0) ? 1 : 3;
  return kronecker(stoi(-z),p);
}

static long
ellrootno_intern(GEN e, GEN p, GEN ex)
{
  if (cmpis(p,3) > 0) return ellrootno_not23(e,p,ex);
  switch(itos(p))
  {
    case 3: return ellrootno_3(e);
    case 2: return ellrootno_2(e);
    default: err(talker,"incorrect prime in ellrootno_intern");
  }
  return 0; /* not reached */
}
	
/* local epsilon factor at p, including p=0 for the infinite place. Global
   if p==1. The equation can be non minimal, but must be over Q. Internal,
   no garbage collection. */
static long
ellrootno_all(GEN e, GEN p, GEN* ptcond)
{
  long s,exs,i;
  GEN fa,gr,cond,pr,ex;

  gr=globalreduction(e);
  e=coordch(e,(GEN)gr[2]);
  cond=(GEN)gr[1]; if(ptcond) *ptcond=cond;
  if (typ(e[12]) != t_INT)
    err(talker,"not an integral curve in ellrootno");
  if (typ(p) != t_INT || signe(p)<0)
    err(talker,"not a nonnegative integer second arg in ellrootno");
  if (cmpis(p,2)>=0)
  {
    exs=ggval(cond,p);
    if (!exs) return 1;
  }
  if (cmpis(p,3)>0) return ellrootno_not23(e,p,stoi(exs));
  switch(itos(p))
  {
    case 3: return ellrootno_3(e);
    case 2: return ellrootno_2(e);
    case 1: s=-1; fa=factor(cond); pr=(GEN)fa[1]; ex=(GEN)fa[2];
      for (i=1; i<lg(pr); i++) s*=ellrootno_intern(e,(GEN)pr[i],(GEN)ex[i]);
      return s;
    case 0: return -1; /* local factor at infinity = -1 */
    default: return 0; /* never reached */
  }
}

long
ellrootno(GEN e, GEN p)
{
  long av=avma,s;
  if (!p) p = gun;
  s=ellrootno_all(e, p, NULL);
  avma=av; return s;
}