OpenXM_contrib/pari-2.2/src/modules/elliptic.c - diff

Return to elliptic.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / modules

Diff for /OpenXM_contrib/pari-2.2/src/modules/Attic/elliptic.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:11

version 1.2, 2002/09/11 07:27:05

Line 26 checkpt(GEN z)

if (typ(z)!=t_VEC) err(elliper1);

}

long

void

checkell(GEN e)

{

long lx=lg(e);

if (typ(e)!=t_VEC || lx<14) err(elliper1);

return lx;

}

void

Line 142 smallinitell0(GEN x, GEN y)

Line 141 smallinitell0(GEN x, GEN y)

GEN

smallinitell(GEN x)

{

ulong av = avma;

gpmem_t av = avma;

GEN y = cgetg(14,t_VEC);

smallinitell0(x,y); return gerepilecopy(av,y);

}

Line 162 ellinit0(GEN x, long flag,long prec)

Line 161 ellinit0(GEN x, long flag,long prec)

void

ellprint(GEN e)

{

long av = avma;

gpmem_t av = avma;

long vx = fetch_var();

long vy = fetch_var();

GEN z = cgetg(3,t_VEC);

Line 170 ellprint(GEN e)

Line 169 ellprint(GEN e)

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]));

fprintferr("%Z - (%Z)\n", ellLHS(e, z), ellRHS(e, polx[vx]));

(void)delete_var();

(void)delete_var(); avma = av;

}

Line 204 static GEN

Line 203 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;

long mi;

if (!x1) x1 = gmul2n(gsub(a1,b1),-2);

mi = min(precp(a1),precp(b1));

for(;;)

{

a=a1; b=b1; x=x1;

b1=gsqrt(gmul(a,b),0); bmod1=modii((GEN)b1[4],p);

b1=gprec(gsqrt(gmul(a,b),0),mi); bmod1=modii((GEN)b1[4],p);

if (!egalii(bmod1,bmod)) b1 = gneg_i(b1);

a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2);

a1=gprec(gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2),mi);

r1=gsub(a1,b1);

if (gcmp0(r1)) {x1=x; break;}

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));

}

Line 227 padic_initell(GEN y, GEN p, long prec)

Line 228 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 */

q=gadd(gun,ggrandocp(p,prec));

for (i=1; i<=13; i++) y[i]=lmul(q,(GEN)y[i]);

if (gcmp0((GEN)y[13]) || valp((GEN)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=stoi(4);

err(impl,"initell for 2-adic numbers");

}

else pv=p;

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];

Line 275 padic_initell(GEN y, GEN p, long prec)

Line 280 padic_initell(GEN y, GEN p, long prec)

y[19]=zero; return y;

}

/* mis pour debugger do_padic_agm. On peut enlever quand on veut */

GEN

dopad(GEN a, GEN b, GEN pv)

{

GEN x1=NULL;

return ginv(gmul2n(do_padic_agm(&x1,a,b,pv),2));

}

static int

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

Line 355 initell0(GEN x, long prec)

Line 370 initell0(GEN x, long prec)

GEN

initell(GEN x, long prec)

{

ulong av = avma;

gpmem_t av = avma;

return gerepilecopy(av, initell0(x,prec));

}

GEN

coordch(GEN e, GEN ch)

_coordch(GEN e, GEN u, GEN r, GEN s, GEN t)

{

GEN y,p1,p2,v,v2,v3,v4,v6,r,s,t,u;

GEN y, p1, p2, v, v2, v3, v4, v6;

long i,lx = checkell(e);

long i, lx = lg(e);

ulong av = avma;

gpmem_t av = avma;

checkch(ch);

y = cgetg(lx,t_VEC);

u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];

v = ginv(u); v2 = gsqr(v); v3 = gmul(v,v2);v4 = gsqr(v2); v6 = gsqr(v3);

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);

Line 388 coordch(GEN e, GEN ch)

Line 401 coordch(GEN e, GEN ch)

y[11] = lmul(v6,(GEN)e[11]);

y[12] = lmul(gsqr(v6),(GEN)e[12]);

y[13] = e[13];

if (lx>14)

if (lx > 14)

{

p1=(GEN)e[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)

{

y[14] = (long)_vec( gmul(v2, gsub((GEN)p1[1],r)) );

y[15] = lmul((GEN)e[15], gsqr(u));

y[16] = lmul((GEN)e[16], u);

y[17] = e[17];

y[18] = e[18]; /* FIXME: how do q and w change ? */

y[19] = zero;

}

else

{

if (typ(e[1])==t_PADIC)

p2 = cgetg(4,t_COL);

{

for (i=1; i<=3; i++) p2[i] = lmul(v2, gsub((GEN)p1[i],r));

p2=cgetg(2,t_VEC); p2[1]=lmul(v2,gsub((GEN)p1[1],r));

y[14] = (long)p2;

y[14]=(long)p2;

y[15] = lmul((GEN)e[15], u);

y[15]=lmul(gsqr(u),(GEN)e[15]);

y[16] = lmul((GEN)e[16], u);

y[16]=lmul(u,(GEN)e[16]);

y[17] = ldiv((GEN)e[17], u);

/* FIXME: how do q and w change ??? */

y[18] = ldiv((GEN)e[18], u);

y[17]=e[17];

y[19] = lmul((GEN)e[19], gsqr(u));

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]);

}

return gerepilecopy(av,y);

}

GEN

coordch(GEN e, GEN ch)

{

checkch(ch); checkell(e);

return _coordch(e, (GEN)ch[1], (GEN)ch[2], (GEN)ch[3], (GEN)ch[4]);

}

static GEN

pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t)

{

Line 431 pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t)

Line 446 pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t)

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

z = cgetg(3,t_VEC); p1=gadd((GEN)x[1],mor);

z = cgetg(3,t_VEC); p1 = gadd((GEN)x[1],mor);

z[1] = lmul(v2,p1);

z[1] = lmul(v2, p1);

z[2] = lmul(v3,gsub((GEN)x[2],gadd(gmul(s,p1),t)));

z[2] = lmul(v3, gsub((GEN)x[2], gadd(gmul(s,p1),t)));

return z;

}

Line 441 GEN

Line 456 GEN

pointch(GEN x, GEN ch)

{

GEN y,v,v2,v3,mor,r,s,t,u;

long tx,lx=lg(x),i;

long tx, i, lx = lg(x);

ulong av = avma;

gpmem_t av = avma;

checkpt(x); checkch(ch);

if (lx < 2) return gcopy(x);

u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];

u = (GEN)ch[1];

tx=typ(x[1]); v=ginv(u); v2=gsqr(v); v3=gmul(v,v2); mor=gneg_i(r);

r = (GEN)ch[2];

s = (GEN)ch[3];

t = (GEN)ch[4];

v = ginv(u); v2 = gsqr(v); v3 = gmul(v,v2);

mor = gneg_i(r);

tx = typ(x[1]);

if (is_matvec_t(tx))

{

y=cgetg(lx,tx);

y = cgetg(lx,tx);

for (i=1; i<lx; i++)

y[i]=(long) pointch0((GEN)x[i],v2,v3,mor,s,t);

y[i] = (long)pointch0((GEN)x[i],v2,v3,mor,s,t);

}

else

y = pointch0(x,v2,v3,mor,s,t);

Line 467 int

Line 487 int

oncurve(GEN e, GEN z)

{

GEN p1,p2,x;

long av=avma,p,q;

long p, q;

gpmem_t av=avma;

checksell(e); checkpt(z); if (lg(z)<3) return 1; /* oo */

p1 = ellLHS(e,z);

Line 485 GEN

Line 506 GEN

addell(GEN e, GEN z1, GEN z2)

{

GEN p1,p2,x,y,x1,x2,y1,y2;

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

checksell(e); checkpt(z1); checkpt(z2);

if (lg(z1)<3) return gcopy(z2);

Line 531 invell(GEN e, GEN z)

Line 552 invell(GEN e, GEN z)

GEN

subell(GEN e, GEN z1, GEN z2)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

checksell(e); checkpt(z2);

z2=invell(e,z2); tetpil=avma;

Line 541 subell(GEN e, GEN z1, GEN z2)

Line 562 subell(GEN e, GEN z1, GEN z2)

GEN

ordell(GEN e, GEN x, long prec)

{

long av=avma,td,i,lx,tx=typ(x);

long td, i, lx, tx=typ(x);

gpmem_t av=avma;

GEN D,a,b,d,pd,p1,y;

checksell(e);

Line 600 static GEN

Line 622 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;

long ln, ep, vn;

gpmem_t av=avma, tetpil;

if (lg(z)<3) return gcopy(z);

pol=(GEN)n[1];

Line 623 CM_powell(GEN e, GEN z, GEN n)

Line 646 CM_powell(GEN e, GEN z, GEN n)

GEN ss=gzero;

{

ep=(-valp(z2))>>1; ss=gadd(ss,gmul((GEN)z2[2],gpuigs(polx[0],ep)));

ep=(-valp(z2))>>1; ss=gadd(ss,gmul((GEN)z2[2],gpowgs(polx[0],ep)));

z2=gsub(z2,gmul((GEN)z2[2],gpuigs(z1,ep)));

z2=gsub(z2,gmul((GEN)z2[2],gpowgs(z1,ep)));

}

while (valp(z2)<=0);

p2=gadd(p0,gmul(ss,p1)); p0=p1; p1=p2;

Line 646 GEN

Line 669 GEN

powell(GEN e, GEN z, GEN n)

{

GEN y;

long av=avma,i,j,tetpil,s;

long i, j, s;

gpmem_t av=avma, tetpil;

ulong m;

checksell(e); checkpt(z);

Line 679 mathell(GEN e, GEN x, long prec)

Line 703 mathell(GEN e, GEN x, long prec)

{

GEN y,p1,p2, *pdiag;

long lx=lg(x),i,j,tx=typ(x);

ulong av = avma;

gpmem_t av = avma;

if (!is_vec_t(tx)) err(elliper1);

lx=lg(x); y=cgetg(lx,t_MAT); pdiag=(GEN*) new_chunk(lx);

Line 704 mathell(GEN e, GEN x, long prec)

Line 728 mathell(GEN e, GEN x, long prec)

static GEN

bilhells(GEN e, GEN z1, GEN z2, GEN h2, long prec)

{

long lz1=lg(z1),tx,av=avma,tetpil,i;

long lz1=lg(z1), tx, i;

gpmem_t av=avma, tetpil;

GEN y,p1,p2;

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

Line 726 GEN

Line 751 GEN

bilhell(GEN e, GEN z1, GEN z2, long prec)

{

GEN p1,h2;

long av=avma,tetpil,tz1=typ(z1),tz2=typ(z2);

long tz1=typ(z1), tz2=typ(z2);

gpmem_t av=avma, tetpil;

if (!is_matvec_t(tz1) || !is_matvec_t(tz2)) err(elliper1);

if (lg(z1)==1) return cgetg(1,tz1);

Line 775 new_coords(GEN e, GEN x, GEN *pta, GEN *ptb, long prec

Line 801 new_coords(GEN e, GEN x, GEN *pta, GEN *ptb, long prec

GEN

zell(GEN e, GEN z, long prec)

{

long av=avma,ty,sw,fl;

long ty, sw, fl;

gpmem_t av=avma;

GEN t,u,p1,p2,a,b,x1,u2,D = (GEN)e[12];

checkbell(e);

Line 857 zell(GEN e, GEN z, long prec)

Line 884 zell(GEN e, GEN z, long prec)

return gerepileupto(av,t);

}

typedef struct {

GEN w1,w2; /* original basis for L = <w1,w2> */

GEN W1,W2,tau; /* new basis for L = <W1,W2> = W2 <1,tau> */

GEN a,b,c,d; /* tau in F = h/Sl2, tau = g.t, g=[a,b;c,d] in SL(2,Z) */

GEN x,y; /* z/w2 defined mod <1,tau> --> z + x tau + y reduced mod <1,tau> */

} SL2_red;

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

static void

getgamma(GEN *ptt)

set_gamma(SL2_red *T)

{

GEN t = *ptt,a,b,c,d,n,m,p1,p2,run;

GEN t = T->tau,a,b,c,d,n,m,p1,run;

run = gsub(realun(DEFAULTPREC), gpuigs(stoi(10),-8));

run = gsub(realun(DEFAULTPREC), gpowgs(stoi(10),-8));

a=d=gun; b=c=gzero;

a = d = gun;

b = c = gzero;

for(;;)

{

n = ground(greal(t));

Line 880 getgamma(GEN *ptt)

Line 915 getgamma(GEN *ptt)

p1=negi(c); c=a; a=p1;

p1=negi(d); d=b; b=p1;

}

m=cgetg(3,t_MAT); *ptt = t;

T->a = a; T->b = b;

p1=cgetg(3,t_COL); m[1]=(long)p1;

T->c = c; T->d = d;

p2=cgetg(3,t_COL); m[2]=(long)p2;

p1[1]=(long)a; p2[1]=(long)b;

p1[2]=(long)c; p2[2]=(long)d; return m;

}

static GEN

#define swap(x,y) { GEN _t=x; x=y; y=_t; }

get_tau(GEN *ptom1, GEN *ptom2, GEN *ptga)

/* swap w1, w2 so that Im(t := w1/w2) > 0. Set tau = representative of t in

* the standard fundamental domain, and g in Sl_2, such that tau = g.t */

static void

red_modSL2(SL2_red *T)

{

GEN om1 = *ptom1, om2 = *ptom2, tau = gdiv(om1,om2);

long s;

long s = gsigne(gimag(tau));

T->tau = gdiv(T->w1,T->w2);

if (!s)

s = gsigne(gimag(T->tau));

err(talker,"omega1 and omega2 R-linearly dependent in elliptic function");

if (!s) err(talker,"w1 and w2 R-linearly dependent in elliptic function");

if (s < 0) { *ptom1=om2; *ptom2=om1; tau=ginv(tau); }

if (s < 0) { swap(T->w1, T->w2); T->tau=ginv(T->tau); }

*ptga = getgamma(&tau); return tau;

set_gamma(T);

/* update lattice */

T->W1 = gadd(gmul(T->a,T->w1), gmul(T->b,T->w2));

T->W2 = gadd(gmul(T->c,T->w1), gmul(T->d,T->w2));

T->tau = gdiv(T->W1, T->W2);

}

static int

get_periods(GEN e, GEN *om1, GEN *om2)

get_periods(GEN e, SL2_red *T)

{

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 3: T->w1=(GEN)e[1]; T->w2=(GEN)e[2]; red_modSL2(T); return 1;

case 20: *om1=(GEN)e[16]; *om2=(GEN)e[15]; return 1;

case 20: T->w1=(GEN)e[16]; T->w2=(GEN)e[15];red_modSL2(T); return 1;

}

return 0;

}

extern GEN PiI2(long prec);

static GEN

check_real(GEN q)

{

return (typ(q) == t_COMPLEX && gcmp0((GEN)q[2]))? (GEN)q[1]: q;

}

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

_elleisnum(SL2_red *T, long k, long prec)

{

long av=avma,lim,av1;

gpmem_t lim, av;

GEN om1,om2,p1,pii2,tau,q,y,qn,ga,court,asub = NULL; /* gcc -Wall */

GEN p1,pii2,q,y,qn,n;

if (k%2 || k<=0) err(talker,"k not a positive even integer in elleisnum");

if (!get_periods(om, &om1, &om2)) err(typeer,"elleisnum");

pii2 = PiI2(prec);

tau = get_tau(&om1,&om2, &ga);

q = gexp(gmul(pii2,T->tau), prec);

if (k==2) asub=gdiv(gmul(pii2,mulsi(12,gcoeff(ga,2,1))),om2);

q = check_real(q);

om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));

y = gzero; n = utoi(1);

if (k==2) asub=gdiv(asub,om2);

av = avma; lim = stack_lim(av,1); qn = gun; n[2] = 0;

q=gexp(gmul(pii2,tau),prec);

y=gzero; court=stoi(3);

av1=avma; lim=stack_lim(av1,1); qn=gun; court[2]=0;

for(;;)

{

court[2]++; qn=gmul(q,qn);

n[2]++; qn = gmul(q,qn);

p1=gdiv(gmul(gpuigs(court,k-1),qn),gsub(gun,qn));

p1 = gdiv(gmul(gpowgs(n,k-1),qn), gsub(gun,qn));

y=gadd(y,p1);

if (gcmp0(p1) || gexpo(p1) <= - bit_accuracy(prec) - 5) break;

if (low_stack(lim, stack_lim(av1,1)))

y = gadd(y, p1);

if (low_stack(lim, stack_lim(av,1)))

{

GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;

if(DEBUGMEM>1) err(warnmem,"elleisnum");

gerepilemany(av1,gptr,2);

gerepilemany(av,gptr,2);

}

y = gadd(gun, gmul(y, gdiv(gdeux, gzeta(stoi(1-k),prec)))); /* = E_k(tau) */

y=gadd(gun,gmul(gdiv(gdeux,gzeta(stoi(1-k),prec)),y));

y = gmul(y, gpowgs(gdiv(pii2,T->W2),k));

p1=gpuigs(gdiv(pii2,om2),k);

return check_real(y);

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

/* Return (2iPi)^k E_k(L) = (2iPi/w2)^k E_k(tau), with L = <w1,w2>, k > 0 even

* E_k(tau) = 1 + 2/zeta(1-k) * sum(n>=1, n^(k-1) q^n/(1-q^n))

* If flag is != 0 and k=4 or 6, compute g2 = E4/12 or g3 = E6/216 resp. */

GEN

elleta(GEN om, long prec)

elleisnum(GEN om, long k, long flag, long prec)

{

long av=avma;

gpmem_t av = avma;

GEN p1, y;

SL2_red T;

if (k&1 || k<=0) err(talker,"k not a positive even integer in elleisnum");

if (!get_periods(om, &T)) err(typeer,"elleisnum");

y = _elleisnum(&T, k, prec);

if (k==2 && signe(T.c))

{

p1 = gmul(PiI2(prec), mulsi(12, T.c));

y = gsub(y, gdiv(p1, gmul(T.w2, T.W2)));

}

else if (k==4 && flag) y = gdivgs(y, 12);

else if (k==6 && flag) y = gdivgs(y,216);

return gerepileupto(av,y);

}

static GEN

_elleta(SL2_red *T, long prec)

{

GEN e2,y1,y2,y;

e2 = gdivgs(elleisnum(om,2,0,prec),12);

e2 = gdivgs(_elleisnum(T,2,prec),12);

y2 = gmul((GEN)om[2],e2);

y2 = gmul(T->W2, e2);

y1 = gadd(gdiv(PiI2(prec),(GEN)om[2]), gmul((GEN)om[1],e2));

y1 = gadd(gdiv(PiI2(prec), T->W2), gmul(T->W1,e2));

y = cgetg(3,t_VEC);

y[1] = lneg(y1);

y[2] = lneg(y2); return gerepileupto(av, y);

y[2] = lneg(y2); return y;

}

/* computes the numerical value of wp(z | om1 Z + om2 Z),

/* compute eta1, eta2 */

If flall=1, compute also wp'. Reduce to the fundamental domain first. */

GEN

elleta(GEN om, long prec)

{

gpmem_t av = avma;

SL2_red T;

if (!get_periods(om, &T)) err(typeer,"elleta");

return gerepileupto(av, _elleta(&T,prec));

}

static GEN

weipellnumall(GEN om1, GEN om2, GEN z, long flall, long prec)

reduce_z(GEN z, long prec, SL2_red *T)

{

long av=avma,tetpil,lim,av1,toadd;

GEN Z = gdiv(z, T->W2);

GEN p1,pii2,pii4,a,tau,q,u,y,yp,u1,u2,qn,v,ga;

long t = typ(z);

pii2 = PiI2(prec);

if (!is_scalar_t(t) || t == t_INTMOD || t == t_PADIC || t == t_POLMOD)

tau = get_tau(&om1,&om2, &ga);

err(typeer,"reduction mod SL2 (reduce_z)");

om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));

T->x = ground(gdiv(gimag(Z), gimag(T->tau)));

z=gdiv(z,om2);

Z = gsub(Z, gmul(T->x,T->tau));

a=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(a,tau));

T->y = ground(greal(Z));

a=ground(greal(z)); z=gsub(z,a);

Z = gsub(Z, T->y);

if (gcmp0(z) || gexpo(z) < 5 - bit_accuracy(prec))

if (gcmp0(Z) || gexpo(Z) < 5 - bit_accuracy(prec)) Z = NULL; /* z in L */

{

return Z;

avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v;

}

q=gexp(gmul(pii2,tau),prec);

/* computes the numerical value of wp(z | L), L = om1 Z + om2 Z

u=gexp(gmul(pii2,z),prec);

* return NULL if z in L. If flall=1, compute also wp' */

u1=gsub(gun,u); u2=gsqr(u1);

static GEN

y=gadd(gdivgs(gun,12),gdiv(u,u2));

weipellnumall(SL2_red *T, GEN z, long flall, long prec)

if (flall) yp=gdiv(gadd(gun,u),gmul(u1,u2));

{

toadd=(long)ceil(9.065*gtodouble(gimag(z)));

long toadd;

gpmem_t av=avma, lim, av1;

GEN p1,pii2,q,u,y,yp,u1,u2,qn,v;

av1=avma; lim=stack_lim(av1,1); qn=q;

z = reduce_z(z,prec, T);

if (!z) return NULL;

/* Now L,z normalized to <1,tau>. z in fund. domain of <1, tau> */

pii2 = PiI2(prec);

q = gexp(gmul(pii2,T->tau),prec);

u = gexp(gmul(pii2,z),prec);

u1= gsub(gun,u); u2=gsqr(u1);

y = gadd(ginv(utoi(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;

for(;;)

{

GEN p2,qnu,qnu1,qnu2,qnu3,qnu4;

GEN qnu,qnu1,qnu2,qnu3,qnu4;

qnu=gmul(qn,u); qnu1=gsub(gun,qnu); qnu2=gsqr(qnu1);

qnu = gmul(qn,u); /* q^n u */

qnu3=gsub(qn,u); qnu4=gsqr(qnu3);

qnu1 = gsub(gun,qnu); /* 1 - q^n u */

p1=gsub(gmul(u,gadd(ginv(qnu2),ginv(qnu4))),

qnu2 = gsqr(qnu1); /* (1 - q^n u)^2 */

gmul2n(ginv(gsqr(gsub(gun,qn))),1));

qnu3 = gsub(qn,u); /* q^n - u */

p1=gmul(qn,p1);

qnu4 = gsqr(qnu3); /* (q^n - u)^2 */

y=gadd(y,p1);

p1 = gsub(gmul(u, gadd(ginv(qnu2),ginv(qnu4))),

gmul2n(ginv(gsqr(gsub(gun,qn))), 1));

y = gadd(y, gmul(qn,p1));

if (flall)

{

p2=gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)),

p1 = gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)),

gdiv(gadd(qn,u),gmul(qnu3,qnu4)));

p2=gmul(qn,p2);

yp=gadd(yp,p2);

yp = gadd(yp, gmul(qn,p1));

}

qn=gmul(q,qn);

qn = gmul(q,qn);

if (gexpo(qn) <= - bit_accuracy(prec) - 5 - toadd) break;

if (low_stack(lim, stack_lim(av1,1)))

{

Line 1025 weipellnumall(GEN om1, GEN om2, GEN z, long flall, lon

Line 1100 weipellnumall(GEN om1, GEN om2, GEN z, long flall, lon

}

pii2=gdiv(pii2,om2);

u1 = gdiv(pii2, T->W2);

pii4=gsqr(pii2);

u2 = gsqr(u1);

y = gmul(pii4,y);

y = gmul(u2,y); /* y *= (2i pi / w2)^2 */

if (flall) yp=gmul(u,gmul(gmul(pii4,pii2),yp));

if (flall)

tetpil=avma;

{

if (flall) { v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lmul2n(yp,-1); }

yp = gmul(u, gmul(gmul(u1,u2),yp));/* yp *= u (2i pi / w2)^3 */

else v=gcopy(y);

v = cgetg(3,t_VEC);

return gerepile(av,tetpil,v);

v[1] = (long)y;

v[2] = lmul2n(yp,-1);

}

else v = y;

return gerepilecopy(av, v);

}

GEN

ellzeta(GEN om, GEN z, long prec)

{

long av=avma,tetpil,lim,av1,toadd;

long toadd;

GEN zinit,om1,om2,p1,pii2,tau,q,u,y,u1,qn,ga,x1,x2,et;

gpmem_t av=avma, lim, av1;

GEN Z,p1,pii2,q,u,y,u1,qn,et;

SL2_red T;

if (!get_periods(om, &om1, &om2)) err(typeer,"ellzeta");

if (!get_periods(om, &T)) err(typeer,"ellzeta");

pii2 = PiI2(prec);

Z = reduce_z(z, prec, &T);

tau = get_tau(&om1,&om2, &ga);

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));

et = _elleta(&T,prec);

x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);

et = gadd(gmul(T.x,(GEN)et[1]), gmul(T.y,(GEN)et[2]));

et=elleta(om,prec);

if (!Z) return gerepileupto(av, gadd(ginv(z),et)); /* FIXME ??? */

et=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));

if (gcmp0(z) || gexpo(z) < 5 - bit_accuracy(prec))

pii2 = PiI2(prec);

{

q=gexp(gmul(pii2,T.tau),prec);

p1=ginv(zinit); tetpil=avma; return gerepile(av,tetpil,gadd(p1,et));

u=gexp(gmul(pii2,Z),prec);

}

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=gdiv(gmul(gsqr(T.W2),_elleisnum(&T,2,prec)),pii2);

y=gadd(ghalf,gdivgs(gmul(z,y),-12));

y=gadd(ghalf,gdivgs(gmul(Z,y),-12));

y=gadd(y,ginv(u1));

toadd=(long)ceil(9.065*gtodouble(gimag(z)));

toadd=(long)ceil(9.065*gtodouble(gimag(Z)));

av1=avma; lim=stack_lim(av1,1); qn=q;

for(;;)

{

Line 1079 ellzeta(GEN om, GEN z, long prec)

Line 1153 ellzeta(GEN om, GEN z, long prec)

}

y=gmul(gdiv(pii2,om2),y);

y=gmul(gdiv(pii2,T.W2),y);

tetpil=avma;

return gerepileupto(av, gadd(y,et));

return gerepile(av,tetpil,gadd(y,et));

}

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

GEN

ellsigma(GEN om, GEN z, long flag, long prec)

ellsigma(GEN w, GEN z, long flag, long prec)

{

long av=avma,lim,av1,toadd;

long toadd;

GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,u1,qn,ga,negu,uinv,x1,x2,et,etnew,uhalf;

gpmem_t av=avma, lim, av1;

GEN Z,zinit,p1,pii2,q,u,y,y1,u1,qn,negu,uinv,et,etnew,uhalf;

int doprod = (flag >= 2);

int dolog = (flag & 1);

SL2_red T;

if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmaprod");

if (!get_periods(w, &T)) err(typeer,"ellsigma");

pii2 = PiI2(prec);

Z = reduce_z(z, prec, &T);

tau = get_tau(&om1,&om2, &ga);

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));

et = _elleta(&T, prec);

x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);

etnew = gadd(gmul(T.x,(GEN)et[1]), gmul(T.y,(GEN)et[2]));

et=elleta(om,prec);

zinit = Z? gmul(Z,T.W2): gzero;

etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));

p1 = gadd(zinit, gmul2n(gadd(gmul(T.x,T.W1),gmul(T.y,T.W2)),-1));

etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit));

etnew = gmul(etnew, p1);

if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1));

pii2 = PiI2(prec);

if (gexpo(z) < 5 - bit_accuracy(prec))

if (mpodd(T.x) || mpodd(T.y)) etnew = gadd(etnew, gmul2n(pii2,-1));

{

if (!Z)

if (dolog)

{ /* FIXME ??? */

return gerepileupto(av, gadd(etnew,glog(zinit,prec)));

if (dolog) Z = gadd(etnew, glog(z,prec));

else

else Z = gmul(gexp(etnew,prec), z);

return gerepileupto(av, gmul(gexp(etnew,prec),zinit));

return gerepileupto(av, Z);

}

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

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

/* 9.065 = 2*Pi/log(2) */

toadd = (long)ceil(9.065*fabs(gtodouble(gimag(z))));

toadd = (long)ceil(9.065*fabs(gtodouble(gimag(Z))));

uhalf = gexp(gmul(gmul2n(pii2,-1),z),prec);

uhalf = gexp(gmul(gmul2n(pii2,-1),Z),prec);

u = gsqr(uhalf);

if (doprod)

{ /* use product */

q=gexp(gmul(pii2,tau),prec);

q=gexp(gmul(pii2,T.tau),prec);

uinv=ginv(u);

u1=gsub(uhalf,ginv(uhalf));

y=gdiv(gmul(om2,u1),pii2);

y=gdiv(gmul(T.W2,u1),pii2);

av1=avma; lim=stack_lim(av1,1); qn=q;

negu=stoi(-1);

for(;;)

Line 1147 ellsigma(GEN om, GEN z, long flag, long prec)

Line 1218 ellsigma(GEN om, GEN z, long flag, long prec)

{ /* use sum */

GEN q8,qn2,urn,urninv;

long n;

q8=gexp(gmul2n(gmul(pii2,tau),-3),prec);

q8=gexp(gmul2n(gmul(pii2,T.tau),-3),prec);

q=gpuigs(q8,8);

q=gpowgs(q8,8);

u=gneg_i(u); uinv=ginv(u);

y=gzero;

av1=avma; lim=stack_lim(av1,1); qn=q; qn2=gun;

Line 1169 ellsigma(GEN om, GEN z, long flag, long prec)

Line 1240 ellsigma(GEN om, GEN z, long flag, long prec)

}

p1=gmul(q8,gmul(gdiv(gdiv((GEN)om[2],pii2),gpuigs(trueeta(tau,prec),3)),y));

p1=gmul(q8,gmul(gdiv(gdiv(T.W2,pii2),gpowgs(trueeta(T.tau,prec),3)),y));

}

if (dolog)

Line 1181 ellsigma(GEN om, GEN z, long flag, long prec)

Line 1252 ellsigma(GEN om, GEN z, long flag, long prec)

GEN

pointell(GEN e, GEN z, long prec)

{

long av=avma,tetpil;

gpmem_t av = avma;

GEN y,yp,v,p1;

GEN v;

SL2_red T;

checkbell(e);

checkbell(e); (void)get_periods(e, &T);

p1=weipellnumall((GEN)e[16],(GEN)e[15],z,1,prec);

v = weipellnumall(&T,z,1,prec);

if (lg(p1)==2) { avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v; }

if (!v) { avma = av; return _vec(gzero); }

y = gsub((GEN)p1[1], gdivgs((GEN)e[6],12));

v[1] = lsub((GEN)v[1], gdivgs((GEN)e[6],12));

yp = gsub((GEN)p1[2], gmul2n(ellLHS0(e,y),-1));

v[2] = lsub((GEN)v[2], gmul2n(ellLHS0(e,(GEN)v[1]),-1));

tetpil=avma; v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lcopy(yp);

return gerepilecopy(av, v);

return gerepile(av,tetpil,v);

}

GEN

static GEN

weipell(GEN e, long prec)

_weipell(GEN c4, GEN c6, long PREC)

{

long av1,tetpil,precres,i,k,l;

long i, k, l, precres = 2*PREC;

GEN res,p1,s,t;

gpmem_t av1, tetpil;

GEN P,p1,s,t, res = cgetg(precres+2,t_SER);

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; }

if (!PREC) { setsigne(res,0); return res; }

for (i=3; i<precres; i+=2) res[i]=zero;

switch(prec)

P = res + 2;

for (i=1; i<precres; i+=2) P[i]=zero;

switch(PREC)

{

default: res[8]=ldivgs((GEN)e[11],6048);

default:P[6] = ldivgs(c6,6048);

case 3: res[6]=ldivgs((GEN)e[10],240);

case 3: P[4] = ldivgs(c4, 240);

case 2: res[4]=zero;

case 2: P[2] = zero;

case 1: res[2]=un;

case 1: P[0] = un;

case 0: break;

}

for (k=4; k<prec; k++)

for (k=4; k<PREC; k++)

{

av1 = avma;

s = k&1? gzero: gsqr((GEN)res[k+2]);

s = k&1? gzero: gsqr((GEN)P[k]);

t = gzero;

for (l=2; l+l<k; l++)

t = gadd(t, gmul((GEN)res[(l+1)<<1],(GEN)res[(k-l+1)<<1]));

t = gadd(t, gmul((GEN)P[l<<1],(GEN)P[(k-l)<<1]));

p1=gmulsg(3,gadd(s,gmul2n(t,1)));

p1 = gmulsg(3,gadd(s,gmul2n(t,1)));

tetpil=avma;

tetpil = avma;

p1=gdivgs(p1,(k-3)*(2*k+1));

p1 = gdivgs(p1, (k-3)*(2*k+1));

res[(k+1)<<1] = lpile(av1,tetpil,p1);

P[k<<1] = lpile(av1,tetpil, p1);

}

return res;

}

GEN

ellwp0(GEN om, GEN z, long flag, long prec, long PREC)

weipell(GEN e, long PREC)

{

GEN v,om1,om2;

GEN c4 = (GEN)e[10];

long av = avma;

GEN c6 = (GEN)e[11];

checkell(e); return _weipell(c4,c6,PREC);

}

if (z==NULL) return weipell(om,PREC);

GEN

weipell0(GEN e, long prec, long PREC)

{

GEN c4,c6;

if (lg(e) > 3)

{

checkell(e);

c4 = (GEN)e[10];

c6 = (GEN)e[11];

}

else

{

c4 = elleisnum(e, 4, 0, prec);

c6 = elleisnum(e, 6, 0, prec); c6 = gneg(c6);

}

return _weipell(c4,c6,PREC);

}

/* assume x a t_POL */

int

is_simple_var(GEN x)

{

return (degpol(x) == 1 && gcmp0((GEN)x[2]) && gcmp1((GEN)x[3]));

}

GEN

ellwp0(GEN w, GEN z, long flag, long prec, long PREC)

{

GEN v;

gpmem_t av = avma;

SL2_red T;

if (!z) return weipell0(w,prec,PREC);

if (typ(z)==t_POL)

{

if (lgef(z) != 4 || !gcmp0((GEN)z[2]) || !gcmp1((GEN)z[3]))

if (!is_simple_var(z)) err(talker,"expecting a simple variable in ellwp");

err(talker,"expecting a simple variable in ellwp");

v = weipell0(w,prec,PREC); setvarn(v, varn(z));

v = weipell(om,PREC); setvarn(v, varn(z));

return v;

}

if (!get_periods(om, &om1, &om2)) err(typeer,"ellwp");

if (!get_periods(w, &T)) err(typeer,"ellwp");

switch(flag)

{

case 0: v=weipellnumall(om1,om2,z,0,prec);

case 0: v = weipellnumall(&T,z,0,prec);

if (typ(v)==t_VEC && lg(v)==2) { avma=av; v=gpuigs(z,-2); }

if (!v) { avma = av; v = gpowgs(z,-2); }

return v;

case 1: v=weipellnumall(om1,om2,z,1,prec);

case 1: v = weipellnumall(&T,z,1,prec);

if (typ(v)==t_VEC && lg(v)==2)

if (!v)

{

GEN p1 = gmul2n(gpuigs(z,3),1);

GEN p1 = gmul2n(gpowgs(z,3),1);

long tetpil=avma;

gpmem_t tetpil = avma;

v=cgetg(3,t_VEC);

v = cgetg(3,t_VEC);

v[1]=lpuigs(z,-2);

v[1] = lpuigs(z,-2);

v[2]=lneg(p1); return gerepile(av,tetpil,v);

v[2] = lneg(p1); return gerepile(av,tetpil,v);

}

return v;

case 2: return pointell(om,z,prec);

case 2: return pointell(w,z,prec);

default: err(flagerr,"ellwp"); return NULL;

}

Line 1266 ellwp0(GEN om, GEN z, long flag, long prec, long PREC)

Line 1372 ellwp0(GEN om, GEN z, long flag, long prec, long PREC)

static GEN

_a_2(GEN e)

{

long av = avma;

gpmem_t 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]);

Line 1282 apell2_intern(GEN e, ulong p)

Line 1388 apell2_intern(GEN e, ulong p)

if (p == 2) return _a_2(e);

else

{

ulong av = avma, i;

ulong i;

gpmem_t av = avma;

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]);

Line 1337 multi_invmod(GEN x, GEN p)

Line 1444 multi_invmod(GEN x, GEN p)

static GEN

addsell(GEN e, GEN z1, GEN z2, GEN p)

{

GEN p1,p2,x,x1,x2,y,y1,y2;

GEN z,p1,p2,x,x1,x2,y,y1,y2;

long av = avma;

gpmem_t av;

if (!z1) return z2;

if (!z2) return z1;

x1 = (GEN)z1[1]; y1 = (GEN)z1[2];

x2 = (GEN)z2[1]; y2 = (GEN)z2[2];

z = cgetg(3, t_VEC); av = avma;

p2 = subii(x2, x1);

if (p2 == gzero)

{

Line 1356 addsell(GEN e, GEN z1, GEN z2, GEN p)

Line 1464 addsell(GEN e, GEN z1, GEN z2, GEN 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);

x = subii(sqri(p1), addii(x1,x2));

y = negi(addii(y1, mulii(p1,subii(x,x1))));

avma = av; p1 = cgetg(3,t_VEC);

x = modii(x,p);

p1[1] = licopy(x);

y = modii(y,p); avma = av;

p1[2] = lmodii(y, p); return p1;

z[1] = licopy(x);

z[2] = licopy(y); return z;

}

/* z1 <-- z1 + z2 */

Line 1387 addsell_part2(GEN e, GEN z1, GEN z2, GEN p, GEN p2inv)

Line 1496 addsell_part2(GEN e, GEN z1, GEN z2, GEN p, GEN p2inv)

}

static GEN

negsell(GEN f)

{

GEN g = cgetg(3, t_VEC);

g[1] = f[1];

g[2] = lnegi((GEN)f[2]);

return g;

}

static GEN

powsell(GEN e, GEN z, GEN n, GEN p)

{

GEN y;

Line 1396 powsell(GEN e, GEN z, GEN n, GEN p)

Line 1514 powsell(GEN e, GEN z, GEN n, GEN p)

if (!s || !z) return NULL;

if (s < 0)

{

n = negi(n); y = cgetg(3,t_VEC);

z = negsell(z);

y[2] = lnegi((GEN)z[2]);

n = negi(n);

y[1] = z[1]; z = y;

}

if (is_pm1(n)) return z;

Line 1417 powsell(GEN e, GEN z, GEN n, GEN p)

Line 1534 powsell(GEN e, GEN z, GEN n, GEN p)

return addsell(e,y,z,p);

}

/* assume H.f = 0, return exact order of f */

static GEN

exact_order(GEN H, GEN f, GEN cp4, GEN p)

{

GEN P, e, g, h = H, fa = decomp(H);

long i, j, l;

P = (GEN)fa[1]; l = lg(P);

e = (GEN)fa[2];

for (i=1; i<l; i++)

for (j=itos((GEN)e[i]); j; j--)

{

g = diviiexact(h,(GEN)P[i]);

if (powsell(cp4,f,g,p)) break;

h = g;

}

return h;

}

/* make sure *x has lgefint >= k */

static void

_fix(GEN x, long k)

Line 1425 _fix(GEN x, long k)

Line 1561 _fix(GEN x, long k)

if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }

}

/* low word of integer x */

/* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */

#define _low(x) (__x=(GEN)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;

long *tx, *ty, *ti, pfinal, i, j, s, k, k1, x, nb;

GEN p1,p2,p3,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,acon,bcon,c4,c6,cp4,pts;

gpmem_t av = avma, av2;

GEN __x;

GEN p1,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,A,B,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;

tx = ty = ti = NULL; /* gcc -Wall */

if (DEBUGLEVEL) (void)timer2();

c4 = gmod(gdivgs((GEN)e[10], -48), p);

c6 = gmod(gdivgs((GEN)e[11], -864), p);

pordmin = gceil(gmul2n(gsqrt(p,DEFAULTPREC),2)); /* ceil( sqrt(4p) ) */

p1p = addsi(1,p);

p2p = shifti(p1p,1);

x = 0; u = c6; k1 = 0; k = kronecker(c6, p);

A = gzero; B = gun; h = p1p;

for(;;)

{

while (flc==flcc || !flc)

while (k == k1 || !k)

{

x++;

u = gadd(c6, gmulsg(x, gaddgs(c4,x*x)));

u = modii(addii(c6, mulsi(x, addii(c4, mulss(x,x)))), p);

flc = kronecker((GEN)u[2],p);

k = kronecker(u, p);

}

flcc = flc;

k1 = k;

f = cgetg(3,t_VEC);

f[1] = (long)lift_intern(gmulsg(x,u));

f[1] = lmodii(mulsi(x,u), p);

f[2] = (long)lift_intern(gsqr(u));

f[2] = lmodii(sqri(u), p);

cp4 = lift_intern(gmul(c4, (GEN)f[2]));

cp4 = modii(mulii(c4, (GEN)f[2]), p);

fh = powsell(cp4,f,h,p);

s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1;

s = itos(gceil(gsqrt(gdiv(pordmin,B),DEFAULTPREC))) >> 1;

nb = min(128, s >> 1);

/* look for h s.t f^h = 0 */

if (bcon == gun)

if (B == gun)

{ /* first time: initialize */

tx = newbloc(s+1);

ty = newbloc(s+1);

ti = newbloc(s+1);

}

else f = powsell(cp4,f,bcon,p); /* F */

else f = powsell(cp4,f,B,p); /* F = B.f */

*tx = evaltyp(t_VECSMALL) | evallg(s+1);

if (!fh) goto FOUND;

p1 = gcopy(fh);

pts = new_chunk(nb+1);

if (s < 3)

{ /* we're nearly done: naive search */

GEN q1 = p1, mf = negsell(f); /* -F */

for (i=1;; i++)

{

p1 = addsell(cp4,p1, f,p); /* h.f + i.F */

if (!p1) { h = addii(h, mulsi( i,B)); goto FOUND; }

q1 = addsell(cp4,q1,mf,p); /* h.f - i.F */

if (!q1) { h = addii(h, mulsi(-i,B)); goto FOUND; }

}

/* Baby Step/Giant Step */

nb = min(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */

pts = cgetg(nb+1, t_VEC);

j = lgefint(p);

for (i=1; i<=nb; i++)

{ /* baby steps */

pts[i] = (long)p1; /* h.f + (i-1).F */

_fix(p1+1, j); tx[i] = _low((GEN)p1[1]);

_fix(p1+1, j); tx[i] = modBIL((GEN)p1[1]);

_fix(p1+2, j); ty[i] = _low((GEN)p1[2]);

_fix(p1+2, j); ty[i] = modBIL((GEN)p1[2]);

p1 = addsell(cp4,p1,f,p); /* h.f + i.F */

if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; }

if (!p1) { h = addii(h, mulsi(i,B)); goto FOUND; }

}

mfh = dummycopy(fh);

mfh = negsell(fh);

mfh[2] = lnegi((GEN)mfh[2]);

fg = addsell(cp4,p1,mfh,p); /* h.f + nb.F - h.f = nb.F */

fg = addsell(cp4,p1,mfh,p); /* nb.F */

if (!fg) { h = mulsi(nb,B); goto FOUND; }

if (!fg) { 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)

Line 1499 apell1(GEN e, GEN p)

Line 1644 apell1(GEN e, GEN p)

if (u[j] == zero) /* sum = 0 or doubling */

{

long k = i+j-2;

if (egalii((GEN)p1[2],(GEN)fg[2])) k -= 2*nb; /* fg = p1 */

if (egalii((GEN)p1[2],(GEN)fg[2])) k -= 2*nb; /* fg == p1 */

h = addii(h, mulsi(k,bcon));

h = addii(h, mulsi(k,B)); goto FOUND;

goto FOUND;

}

v = multi_invmod(u, p);

Line 1510 apell1(GEN e, GEN p)

Line 1654 apell1(GEN e, GEN p)

{

p1 = (GEN)pts[j];

addsell_part2(cp4,p1,fg,p, (GEN)v[j]);

tx[i] = _low((GEN)p1[1]);

tx[i] = modBIL((GEN)p1[1]);

ty[i] = _low((GEN)p1[2]);

ty[i] = modBIL((GEN)p1[2]);

}

avma = av2;

}

p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = f^(s-1) */

p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = (s-1).F */

if (!p1) { h = mulsi(s-1,B); goto FOUND; }

if (DEBUGLEVEL) msgtimer("[apell1] baby steps, s = %ld",s);

/* giant steps: fg = f^s */

/* giant steps: fg = s.F */

fg = addsell(cp4,p1,f,p);

if (!fg) { h = addii(h, mulsi(s,bcon)); goto FOUND; }

if (!fg) { h = mulsi(s,B); goto FOUND; }

pfinal = _low(p); av2 = avma;

pfinal = modBIL(p); av2 = avma;

p1 = gen_sort(tx, cmp_IND | cmp_C, NULL);

for (i=1; i<=s; i++) ti[i] = tx[p1[i]];

Line 1531 apell1(GEN e, GEN p)

Line 1676 apell1(GEN e, GEN p)

avma = av2;

gaffect(fg, (GEN)pts[1]);

for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */

for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */

gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]);

{

p1 = addsell(cp4,(GEN)pts[j-1],fg,p);

if (!p1) { h = mulii(mulss(s,j), B); goto FOUND; }

gaffect(p1, (GEN)pts[j]);

}

/* replace fg by nb.fg since we do nb points at a time */

avma = av2;

fg = gcopy((GEN)pts[nb]);

Line 1542 apell1(GEN e, GEN p)

Line 1691 apell1(GEN e, GEN p)

{

GEN ftest = (GEN)pts[j];

ulong m, l = 1, r = s+1;

long k, k2;

long k, k2, j2;

avma = av2;

k = _low((GEN)ftest[1]);

k = modBIL((GEN)ftest[1]);

while (l<r)

{

m = (l+r) >> 1;

Line 1554 apell1(GEN e, GEN p)

Line 1703 apell1(GEN e, GEN p)

if (r <= (ulong)s && tx[r] == k)

{

while (tx[r] == k && r) r--;

k2 = _low((GEN)ftest[2]);

k2 = modBIL((GEN)ftest[2]);

for (r++; tx[r] == k && r <= (ulong)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;

if (DEBUGLEVEL) msgtimer("[apell1] giant steps, i = %ld",i);

p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p);

if (egalii((GEN)p1[1], (GEN)ftest[1]))

{

if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i;

h = addii(h, mulii(addis(mulss(s,i), j2), bcon));

h = addii(h, mulii(addis(mulss(s,i), j2), B));

goto FOUND;

}

if (++j > nb)

{ /* compute next nb points */

long save = 0; /* gcc -Wall */

long save = 0; /* gcc -Wall */;

for (j=1; j<=nb; j++)

{

p1 = (GEN)pts[j];

Line 1590 apell1(GEN e, GEN p)

Line 1739 apell1(GEN e, GEN p)

}

FOUND: /* success, found a point of exponent h */

FOUND: /* found a point of exponent h */

p2 = decomp(h); p1=(GEN)p2[1]; p2=(GEN)p2[2];

h = exact_order(h, f,cp4,p);

for (i=1; i<lg(p1); i++)

/* now h is the exact order, and divides E(Fp) = A (mod B) */

for (j=itos((GEN)p2[i]); j; j--)

if (B == gun) B = h;

{

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));

p1 = chinois(gmodulcp(A,B), gmodulsg(0,h));

acon = (GEN)p1[2];

A = (GEN)p1[2];

bcon = (GEN)p1[1];

B = (GEN)p1[1];

}

i = (cmpii(bcon, pordmin) < 0);

i = (cmpii(B, pordmin) < 0);

if (i) acon = centermod(subii(p2p,acon), bcon);

if (i) A = centermod(subii(p2p,A), B);

p1 = ground(gdiv(gsub(p1p,acon),bcon));

p1 = diviiround(gsub(p1p,A), B);

h = addii(acon, mulii(p1,bcon));

h = addii(A, mulii(p1,B));

/* h = A mod B, closest lift to p+1 */

if (!i) break;

}

gunclone(tx);

gunclone(ty);

gunclone(ti);

p1 = (flc==1)? subii(p1p,h): subii(h,p1p);

p1 = (k==1)? subii(p1p,h): subii(h,p1p);

return gerepileupto(av,p1);

return gerepileuptoint(av,p1);

}

typedef struct

Line 1667 powssell1(long e, long p, long n, sellpt *p1, sellpt *

Line 1810 powssell1(long e, long p, long n, sellpt *p1, sellpt *

}

/* assume H.f = 0, return exact order of f, cf. exact_order */

static long

sexact_order(long H, sellpt *f, long cp4, long p)

{

GEN P, e, fa = decomp(stoi(H));

long h = H, g, pp;

long i, j, l;

sellpt fh;

P = (GEN)fa[1]; l = lg(P);

e = (GEN)fa[2];

for (i=1; i<l; i++)

{

pp = itos((GEN)P[i]);

for (j=itos((GEN)e[i]); j; j--)

{

g = h / pp;

powssell1(cp4, p, g, f, &fh);

if (!fh.isnull) break;

h = g;

}

return h;

}

typedef struct

{

long x,y,i;

Line 1682 compare_multiples(multiple *a, multiple *b)

Line 1850 compare_multiples(multiple *a, multiple *b)

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 pordmin,u,p1p,p2p,A,B,c4,c6,cp4;

long av,i,j,s,flc,flcc,x,q,h,p3,l,r,m;

long i, s, k, k1, x, q, h, l, r, m;

gpmem_t av;

multiple *table;

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

table = NULL; /* gcc -Wall */

av = avma; p1 = gmodulss(1,p);

av = avma;

c4 = itos((GEN)gdivgs(gmul(p1,(GEN)e[10]), -48)[2]);

c4 = itos( gmodgs(gdivgs((GEN)e[10], -48), p) );

c6 = itos((GEN)gdivgs(gmul(p1,(GEN)e[11]), -864)[2]);

c6 = itos( gmodgs(gdivgs((GEN)e[11], -864), p) );

pordmin = (long)(1 + 4*sqrt((float)p));

p1p = p+1; p2p = p1p << 1;

p1p = p+1;

x=0; flcc=0; flc = kross(c6, p);

p2p = p1p << 1;

u=c6; acon=0; bcon=1; h=p1p;

x = 0; u = c6; k1 = 0; k = kross(c6, p);

table = NULL; /* gcc -Wall */

A = 0; B = 1; h = p1p;

for(;;)

{

while (flc==flcc || !flc)

while (k == k1 || !k)

{

x++;

u = addssmod(c6, mulssmod(x, c4+mulssmod(x,x,p), p), p);

flc = kross(u,p);

k = kross(u,p);

}

flcc = flc;

k1 = k;

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);

s = (long) (sqrt(((float)pordmin)/B) / 2);

if (!s) s=1;

if (bcon==1)

if (B==1)

{

table = (multiple *) gpmalloc((s+1)*sizeof(multiple));

f2 = f;

}

else powssell1(cp4, p, bcon, &f, &f2);

else powssell1(cp4, p, B, &f, &f2);

for (i=0; i < s; i++)

{

if (fh.isnull) { h += bcon*i; goto FOUND; }

if (fh.isnull) { h += B*i; goto FOUND; }

table[i].x = fh.x;

table[i].y = fh.y;

table[i].i = i;

Line 1742 apell0(GEN e, long p)

Line 1911 apell0(GEN e, long p)

if (r < s && table[r].x == ftest.x) break;

addsell1(cp4, p, &ftest, &fg);

}

h += table[r].i * bcon;

h += table[r].i * B;

if (table[r].y == ftest.y) i = -i;

h += s * i * bcon;

h += s * i * B;

FOUND:

p2=decomp(stoi(h)); p1=(GEN)p2[1]; p2=(GEN)p2[2];

h = sexact_order(h, &f, cp4, p);

for (i=1; i < lg(p1); i++)

if (B == 1) B = h;

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));

GEN p1 = chinois(gmodulss(A,B), gmodulss(0,h));

acon = itos((GEN)p1[2]);

A = itos((GEN)p1[2]);

if (is_bigint(p1[1])) { h = acon; break; }

if (is_bigint(p1[1])) { h = A; break; }

bcon = itos((GEN)p1[1]);

B = itos((GEN)p1[1]);

}

i = (bcon < pordmin);

i = (B < pordmin);

if (i)

{

acon = (p2p - acon) % bcon;

A = (p2p - A) % B;

if ((acon << 1) > bcon) acon -= bcon;

if ((A << 1) > B) A -= B;

}

q = ((ulong)(p2p + bcon - (acon << 1))) / (bcon << 1);

q = ((ulong)(p2p + B - (A << 1))) / (B << 1);

h = acon + q*bcon;

h = A + q*B;

avma = av; if (!i) break;

}

free(table); return stoi((flc==1)? p1p-h: h-p1p);

free(table); return stoi(k==1? p1p-h: h-p1p);

}

GEN

Line 1785 apell(GEN e, GEN p)

Line 1946 apell(GEN e, GEN p)

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;

long s;

gpmem_t av = avma;

GEN p0 = egalii(p,gdeux)? stoi(8): p;

GEN c6 = gmul((GEN)e[11],gmodulsg(1,p0));

s = kronecker(lift_intern(c6),p); avma=av;

Line 1809 GEN

Line 1971 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;

long p, i, m;

gpmem_t av, tetpil;

GEN p1,p2,an;

checkell(e);

Line 1873 anell(GEN e, long n)

Line 2036 anell(GEN e, long n)

GEN

akell(GEN e, GEN n)

{

long i,j,ex,av=avma;

long i, j, ex;

gpmem_t av=avma;

GEN p1,p2,ap,u,v,w,y,pl;

checkell(e);

Line 1908 akell(GEN e, GEN n)

Line 2072 akell(GEN e, GEN n)

GEN

hell(GEN e, GEN a, long prec)

{

long av=avma,tetpil,n;

long n;

gpmem_t av=avma, tetpil;

GEN p1,p2,y,z,q,pi2surw,pi2isurw,qn,ps;

checkbell(e);

Line 1951 hells(GEN e, GEN x, long prec)

Line 2116 hells(GEN e, GEN x, long prec)

return mu;

}

/* [1,0,0,0] */

static GEN

init_ch() {

GEN v = cgetg(5,t_VEC);

v[1] = un; v[2] = v[3] = v[4] = zero; return v;

}

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;

long lx, lc, i, j, tx;

gpmem_t av=avma, tetpil;

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 = init_ch(); p1[2] = laddgs(gfloor(e3),-1);

p1[3]=p1[4]=zero; ep=coordch(e,p1); xp=pointch(x,p1);

ep = coordch(e,p1);

xp = pointch(x,p1);

tx=typ(x[1]); lx=lg(x);

if (!is_matvec_t(tx))

{

Line 2016 hell0(GEN e, GEN z, long prec)

Line 2190 hell0(GEN e, GEN z, long prec)

static GEN

ghell0(GEN e, GEN a, long flag, long prec)

{

long av=avma,lx,i,n,n2,grandn,tx=typ(a);

long lx, i, n, n2, grandn, tx=typ(a);

gpmem_t av=avma;

GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;

checkbell(e); if (!is_matvec_t(tx)) err(elliper1);

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

Line 2281 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;

long l, n, eps, flun;

gpmem_t av=avma, av1, tetpil, lim;

GEN z,p1,p2,cg,cg1,v,cga,cgb,s2,ns,gs,N;

if (!A) A = gun;

Line 2118 lseriesell(GEN e, GEN s, GEN A, long prec)

Line 2294 lseriesell(GEN e, GEN s, GEN A, long prec)

}

flun = gcmp1(A) && gcmp1(s);

eps = ellrootno_all(e,gun,&N);

if (flun && eps<0) { z=cgetr(prec); affsr(0,z); return z; }

if (flun && eps < 0) return realzero(prec);

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((ulong)l,TEMPMAX));

s2 = ns = NULL; /* gcc -Wall */

if (!flun) { s2=gsubsg(2,s); ns=gpui(cg,gsubgs(gmul2n(s,1),2),prec); }

if (!flun) { s2=gsubsg(2,s); ns=gpow(cg,gsubgs(gmul2n(s,1),2),prec); }

z=gzero;

if (typ(s)==t_INT)

{

Line 2136 lseriesell(GEN e, GEN s, GEN A, long prec)

Line 2313 lseriesell(GEN e, GEN s, GEN A, long 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));

p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpow(stoi(n),s,prec));

p2=flun? p1: gdiv(gmul(ns,incgam(s2,gmulsg(n,cgb),prec)),

gpui(stoi(n),s2,prec));

gpow(stoi(n),s2,prec));

if (eps<0) p2=gneg_i(p2);

z = gadd(z, gmul(gadd(p1,p2),

((ulong)n<=TEMPMAX)? (GEN)v[n]: akell(e,stoi(n))));

Line 2177 lseriesell(GEN e, GEN s, GEN A, long prec)

Line 2354 lseriesell(GEN e, GEN s, GEN A, long prec)

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 void cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t);

static GEN

localreduction_result(long av, long f, long kod, long c, GEN v)

localred_result(long f, long kod, long c, GEN v)

{

long tetpil = avma;

GEN z = cgetg(5, t_VEC);

GEN result = cgetg(5, t_VEC);

z[1] = lstoi(f);

result[1] = lstoi(f); result[2] = lstoi(kod);

z[2] = lstoi(kod);

result[3] = lcopy(v); result[4] = lstoi(c);

z[3] = lcopy(v);

return gerepile(av,tetpil, result);

z[4] = lstoi(c); return z;

}

/* ici, p != 2 et p != 3 */

/* Here p > 3. e assumed integral */

static GEN

localreduction_carac_not23(GEN e, GEN p)

localred_carac_p(GEN e, GEN p, int minim)

{

long av = avma, k, f, kod, c, nuj, nudelta;

long k, f, kod, c, nuj, nudelta;

GEN pk, p2k, a2prime, a3prime;

GEN p2, v = init_ch();

GEN p2, r = gzero, s = gzero, t = gzero, v;

GEN c4, c6, delta, tri;

GEN c4, c6, delta, unmodp, xun, tri, var, p4k, p6k;

nudelta = ggval((GEN)e[12], p);

c4 = (GEN)e[10];

v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;

c6 = (GEN)e[11];

nuj = gcmp0((GEN)e[13])? 0: - ggval((GEN)e[13], p);

delta = (GEN)e[12];

nuj = gcmp0((GEN)e[13])? 0: - ggval((GEN)e[13], p);

nudelta = ggval(delta, p);

k = (nuj > 0 ? nudelta - nuj : nudelta) / 12;

c4 = (GEN)e[10]; c6 = (GEN)e[11]; delta = (GEN)e[12];

if (k <= 0)

if (k > 0) /* modele non minimal */

{

pk = gpuigs(p, k);

if (minim) return v;

if (mpodd((GEN)e[1]))

}

s = shifti(subii(pk, (GEN)e[1]), -1);

else

{ /* model not minimal */

s = negi(shifti((GEN)e[1], -1));

GEN pk = gpowgs(p,k), p2k = sqri(pk), p4k = sqri(p2k), p6k = mulii(p4k,p2k);

p2k = sqri(pk);

GEN r, s, t;

p4k = sqri(p2k);

p6k = mulii(p4k, p2k);

a2prime = subii((GEN)e[2], mulii(s, addii((GEN)e[1], s)));

s = negi((GEN)e[1]);

switch(smodis(a2prime, 3))

if (mpodd(s)) s = addii(s, pk);

s = shifti(s, -1);

r = subii(mulii(s, addii((GEN)e[1], s)), (GEN)e[2]); /* - a_2' */

switch(smodis(r, 3))

{

case 0: r = negi(divis(a2prime, 3)); break;

default: break; /* 0 */

case 1: r = divis(subii(p2k, a2prime), 3); break;

case 1: r = addii(r, p2k); break;

case 2: r = negi(divis(addii(a2prime, p2k), 3)); break;

case 2: r = subii(r, p2k); break;

}

a3prime = ellLHS0_i(e,r);

r = divis(r, 3);

if (mpodd(a3prime))

t = shifti(subii(mulii(pk, p2k), a3prime), -1);

t = negi(ellLHS0_i(e,r)); /* - a_3' */

else

if (mpodd(t)) t = addii(t, mulii(pk, p2k));

t = negi(shifti(a3prime, -1));

t = shifti(t, -1);

v[1] = (long)pk; v[2] = (long)r; v[3] = (long)s; v[4] = (long)t;

v[1] = (long)pk;

v[2] = (long)r;

v[3] = (long)s;

v[4] = (long)t;

if (minim) return v;

nudelta -= 12 * k;

c4 = divii(c4, p4k); c6 = divii(c6, p6k);

c4 = diviiexact(c4, p4k);

delta = divii(delta, sqri(p6k));

c6 = diviiexact(c6, p6k);

delta = diviiexact(delta, sqri(p6k));

}

if (nuj > 0) switch(nudelta - nuj)

{

case 0: f = 1; kod = 4+nuj; /* Inu */

Line 2244 localreduction_carac_not23(GEN e, GEN p)

Line 2431 localreduction_carac_not23(GEN e, GEN p)

break;

case 6: f = 2; kod = -4-nuj; /* Inu* */

if (nuj & 1)

c = 3 + kronecker(divii(mulii(c6, delta),gpuigs(p, 9+nuj)), p);

c = 3 + kronecker(divii(mulii(c6, delta),gpowgs(p, 9+nuj)), p);

else

c = 3 + kronecker(divii(delta, gpuigs(p, 6+nuj)), p);

c = 3 + kronecker(divii(delta, gpowgs(p, 6+nuj)), p);

break;

default: err(bugparier,"localred (nu_delta - nu_j != 0,6)");

return NULL; /* not reached */

}

else switch(nudelta)

{

case 0: f = 0; kod = 1; c = 1; break; /* I0, regulier */

case 0: f = 0; kod = 1; c = 1; break; /* I0, regular */

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(p)), p);

c = 2 + kronecker(divii(mulis(c6, -6), sqri(p)), p);

break;

case 6: f = 2; kod = -1; /* I0* */

p2 = sqri(p);

unmodp = gmodulsg(1,p);

/* x^3 - 3c4/p^2 x - 2c6/p^3 */

var = gmul(unmodp,polx[0]);

tri = coefs_to_pol(4, gun, gzero,

tri = gsub(gsqr(var),gmul(divii(gmulsg(3, c4), p2),unmodp));

negi(divii(gmul2n(c6,1), mulii(p2,p))),

tri = gsub(gmul(tri, var),

negi(divii(gmulsg(3, c4), p2)));

gmul(divii(gmul2n(c6,1), mulii(p2,p)),unmodp));

c = 1 + FpX_nbroots(tri, p);

xun = gmodulcp(var,tri);

c = lgef(ggcd((GEN)(gsub(gpui(xun,p,0),xun))[2], tri)) - 2;

break;

case 8: f = 2; kod = -4; /* IV* */

c = 2 + kronecker(gdiv(mulis(c6,-6), gpuigs(p,4)), p);

c = 2 + kronecker(gdiv(mulis(c6,-6), gpowgs(p,4)), p);

break;

case 9: f = 2; kod = -3; c = 2; break; /* III* */

case 10: f = 2; kod = -2; c = 1; break; /* II* */

default: err(bugparier,"localred");

return NULL; /* not reached */

}

return localreduction_result(av,f,kod,c,v);

return localred_result(f, kod, c, v);

}

/* renvoie a_{ k,l } avec les notations de Tate */

/* return a_{ k,l } in Tate's notation */

static int

aux(GEN ak, int p, int l)

{

long av = avma, pl = p, res;

gpmem_t av = avma;

while (--l) pl *= p;

long res = smodis(divis(ak, u_pow(p, l)), p);

res = smodis(divis(ak, pl), p);

avma = av; return res;

}

static int

aux2(GEN ak, int p, GEN pl)

{

long av = avma, res;

gpmem_t av = avma;

res = smodis(divii(ak, pl), p);

long res = smodis(divii(ak, pl), p);

avma = av;

avma = av; return res;

return res;

}

/* renvoie le nombre de racines distinctes du polynome XXX + aXX + bX + c

/* number of distinct roots of X^3 + aX^2 + bX + c modulo p

* modulo p s'il y a une racine multiple, elle est renvoyee dans *mult

* if there's a multiple root, put it un *mult */

static int

numroots3(int a, int b, int c, int p, int *mult)

{

Line 2317 numroots3(int a, int b, int c, int p, int *mult)

Line 2499 numroots3(int a, int b, int c, int p, int *mult)

}

/* idem pour aXX +bX + c */

/* same for aX^2 +bX + c */

static int

numroots2(int a, int b, int c, int p, int *mult)

{

Line 2325 numroots2(int a, int b, int c, int p, int *mult)

Line 2507 numroots2(int a, int b, int c, int p, int *mult)

else { *mult = a * b; return (b * b - a * c) % 3 ? 2 : 1; }

}

/* ici, p1 = 2 ou p1 = 3 */

/* p = 2 or 3 */

static GEN

localreduction_carac_23(GEN e, GEN p1)

localred_carac_23(GEN e, long p, int minim)

{

long av = avma, p, c, nu, nudelta;

long c, nu, nudelta;

int a21, a42, a63, a32, a64, theroot, al, be, ga;

int a21, a42, a63, a32, a64, theroot, al, be, ga, p2, p3, p4;

GEN pk, p2k, pk1, p4, p6;

GEN pk, p2k, pk1;

GEN p2, p3, r = gzero, s = gzero, t = gzero, v;

GEN r, s, t, v;

nudelta = ggval((GEN)e[12], p1);

nudelta = ggval((GEN)e[12], stoi(p));

v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;

v = init_ch();

for(;;)

for (;;)

{

if (!nudelta)

return localreduction_result(av, 0, 1, 1, v);

return localred_result(0, 1, 1, v);

/* I0 */

p = itos(p1);

if (smodis((GEN)e[6], p)) /* p \nmid b2 */

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);

return localred_result(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);

if (signe(r))

else t = modis((GEN)e[5], 2);

t = modis(addii(addii((GEN)e[4], (GEN)e[5]), s), 2);

else

t = modis((GEN)e[5], 2);

}

else /* p == 3 */

{

Line 2365 localreduction_carac_23(GEN e, GEN p1)

Line 2548 localreduction_carac_23(GEN e, GEN p1)

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

cumule(&v, &e, gun, r, s, t); /* p | (a1, a2, a3, a4, a6) */

p2 = stoi(p*p);

p2 = p * p;

if (!divise((GEN)e[5], p2))

if (smodis((GEN)e[5], p2))

return localreduction_result(av, nudelta, 2, 1, v);

return localred_result(nudelta, 2, 1, v);

/* II */

p3 = stoi(p*p*p);

p3 = p2 * p;

if (!divise((GEN)e[9], p3))

if (smodis((GEN)e[9], p3))

return localreduction_result(av, nudelta - 1, 3, 2, v);

return localred_result(nudelta - 1, 3, 2, v);

/* III */

if (!divise((GEN)e[8], p3))

if (smodis((GEN)e[8], p3))

{

if (smodis((GEN)e[8], (p==2)? 32: 27) == p*p)

if (smodis((GEN)e[8], (p==2)? 32: 27) == p2)

c = 3;

else

c = 1;

return localreduction_result(av, nudelta - 2, 4, c, v);

return localred_result(nudelta - 2, 4, c, v);

}

/* IV */

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

if (smodis((GEN)e[5], p3))

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);

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);

return localred_result(nudelta - 4, -1, c, v);

/* I0* */

case 2: /* calcul de nu */

case 2: /* compute 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;

pk = stoi(p2);

p2k = stoi(p * p * p * p);

p2k = stoi(p2 * p2);

for(;;)

{

if (numroots2(al = 1, be = aux2((GEN)e[3], p, pk),

be = aux2((GEN)e[3], p, pk);

ga = -aux2((GEN)e[5], p, p2k), p, &theroot) == 2)

ga = -aux2((GEN)e[5], p, p2k);

break;

al = 1;

if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk));

if (numroots2(al, be, ga, p, &theroot) == 2) break;

pk1 = pk; pk = mulsi(p, pk); p2k = mulsi(p, p2k);

if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk));

nu++;

pk1 = pk;

if (numroots2(al = a21, be = aux2((GEN)e[4], p, pk),

pk = mulsi(p, pk);

ga = aux2((GEN)e[5], p, p2k), p, &theroot) == 2)

p2k = mulsi(p, p2k); nu++;

break;

if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero);

al = a21;

p2k = mulsi(p, p2k);

be = aux2((GEN)e[4], p, pk);

nu++;

ga = aux2((GEN)e[5], p, p2k);

}

if (numroots2(al, be, ga, p, &theroot) == 2) break;

if (p == 2)

if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero);

c = 4 - 2 * (ga & 1);

p2k = mulsi(p, p2k); nu++;

else

}

c = 3 + kross(be * be - al * ga, 3);

if (p == 2)

return localreduction_result(av, nudelta - 4 - nu, -4 - nu, c, v);

c = 4 - 2 * (ga & 1);

/* Inu* */

else

c = 3 + kross(be * be - al * ga, 3);

return localred_result(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);

a32 = aux((GEN)e[3], p, 2);

if (numroots2(1, a32, -a64, p, &theroot) == 2)

a64 = aux((GEN)e[5], p, 4);

{

if (numroots2(1, a32, -a64, p, &theroot) == 2)

if (p == 2)

{

c = 3 - 2 * a64;

if (p == 2)

else

c = 3 - 2 * a64;

c = 2 + kross(a32 * a32 + a64, 3);

else

return localreduction_result(av, nudelta - 6, -4, c, v);

c = 2 + kross(a32 * a32 + a64, 3);

}

return localred_result(nudelta - 6, -4, c, v);

/* IV* */

}

if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p));

/* IV* */

/* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */

if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p));

p4 = sqri(p2);

/* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */

if (!divise((GEN)e[4], p4))

p4 = p2 * p2;

return localreduction_result(av, nudelta - 7, -3, 2, v);

if (smodis((GEN)e[4], p4))

/* III* */

return localred_result(nudelta - 7, -3, 2, v);

p6 = mulii(p4, p2);

/* III* */

if (!divise((GEN)e[5], p6))

return localreduction_result(av, nudelta - 8, -2, 1, v);

if (smodis((GEN)e[5], p4 * p2))

/* II* */

return localred_result(nudelta - 8, -2, 1, v);

cumule(&v, &e, p1, gzero, gzero, gzero); /* non minimal, on repart

/* II* */

pour un tour */

cumule(&v, &e, stoi(p), gzero, gzero, gzero); /* not minimal */

nudelta -= 12;

}

/* Not reached */

}

static GEN

localred(GEN e, GEN p, int minim)

{

if (gcmpgs(p, 3) > 0) /* p != 2,3 */

return localred_carac_p(e,p, minim);

else

{

GEN z = localred_carac_23(e,itos(p), minim);

return minim? (GEN)z[3]: z;

}

GEN

localreduction(GEN e, GEN p1)

localreduction(GEN e, GEN p)

{

gpmem_t av = avma;

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 gerepileupto(av, localred(e, p, 0));

return localreduction_carac_not23(e,p1);

else

return localreduction_carac_23(e,p1);

}

#if 0

static GEN

/* Calcul de toutes les fibres non elliptiques d'une courbe sur Z.

ell_to_small(GEN E)

* 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;

long i;

GEN p1, p2, p3, p4, prims, result;

GEN e;

if (lg(E) <= 14) return E;

e = cgetg(14,t_VEC);

for (i = 1; i < 14; i++) e[i] = E[i];

return e;

}

static GEN

ellintegralmodel(GEN e)

{

GEN a = cgetg(6,t_VEC), v, prims, d, u;

long i, l, k;

checkell(e);

prims = decomp((GEN)e[12]);

for (i = 1; i < 6; i++)

l = lg(p1 = (GEN)prims[1]);

{

p2 = (GEN)prims[2];

a[i] = e[i];

if ((long)prims == avma) cgiv(prims);

switch(typ(a[i]))

result = cgetg(5, t_MAT);

{

result[1] = (long)p1;

case t_INT: case t_FRAC: case t_FRACN: break;

result[2] = (long)p2;

default: err(talker, "not a rational curve in ellintegralmodel");

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));

/* a = [a1, a2, a3, a4, a6] */

for (k = 1; k < l; k++) p4[k] = lgeti(3);

d = denom(a);

av = avma;

prims = (GEN)auxdecomp(d, 0)[1]; /* partial factorization */

if (is_pm1(d)) return NULL;

l = lg(prims); u = gun;

for (k = 1; k < l; k++)

{

GEN q = localreduction(e, (GEN)p1[k]);

GEN p = (GEN)prims[k];

affii((GEN)q[1],(GEN)p2[k]);

int n = 0, m;

affii((GEN)q[2],(GEN)p3[k]);

for (i = 1; i < 6; i++)

affii((GEN)q[4],(GEN)p4[k]);

if (!gcmp0((GEN)a[i]))

avma = av;

{

int r = (i == 5)? 6: i; /* a5 is missing */

m = r * n + ggval((GEN)a[i], p);

while (m < 0) { n++; m += r; }

}

u = mulii(u, gpowgs(p, n));

}

return result;

v = init_ch(); v[1] = linv(u); return v;

}

#endif

/* Algorithme de reduction d'une courbe sur Q a sa forme standard. Etant

static void

* donne une courbe elliptique sous forme longue e, dont les coefficients

standard_model(GEN e, GEN *pv)

* 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

GEN r, s, t;

* e a sa forme minimale globale dans laquelle a1 et a3 valent 0 ou 1, et a2

s = gdiventgs((GEN)e[1], -2);

* vaut -1, 0 ou 1 et tel que u est un rationnel positif. Enfin c est le

r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3);

* produit des nombres de Tamagawa locaux cp.

t = gdiventgs(ellLHS0(e,r), -2);

cumulev(pv, gun, r, s, t);

}

GEN

globalreduction(GEN e1)

ellminimalmodel(GEN E, GEN *ptv)

{

long i, k, l, m, tetpil, av = avma;

gpmem_t av = avma;

GEN p1, c = gun, prims, result, N = gun, u = gun, r, s, t;

GEN c4, c6, e, v, prims;

GEN v = cgetg(5, t_VEC), a = cgetg(7, t_VEC), e = cgetg(20, t_VEC);

long l, k;

checkell(e1);

v = ellintegralmodel(E);

for (i = 1; i < 5; i++) a[i] = e1[i]; a[5] = zero; a[6] = e1[5];

e = ell_to_small(E);

prims = decomp(denom(a));

if (v) e = coordch(e, v); else v = init_ch();

p1 = (GEN)prims[1]; l = lg(p1);

c4 = (GEN)e[10];

c6 = (GEN)e[11];

prims = (GEN)decomp(mppgcd(c4,c6))[1];

l = lg(prims);

for (k = 1; k < l; k++)

{

int n = 0;

GEN w = localred(e, (GEN)prims[k], 1);

for (i = 1; i < 7; i++)

if (!gcmp1((GEN)w[1]))

if (!gcmp0((GEN)a[i]))

cumule(&v, &e, (GEN)w[1], (GEN)w[2], (GEN)w[3], (GEN)w[4]);

{

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;

standard_model(e, &v);

for (i = 1; i < 14; i++) e[i] = e1[i];

e = coordch(E, v);

for (; i < 20; i++) e[i] = zero;

if (ptv) { gerepileall(av, 2, &e, &v); *ptv = v; }

if (!gcmp1(u)) e = coordch(e, v);

else e = gerepileupto(av, e);

prims = decomp((GEN)e[12]);

return e;

l = lg(p1 = (GEN)prims[1]);

}

for (k = (signe(e[12]) < 0) + 1; k < l; k++)

/* Reduction of a rational curve E to its standard minimal model

* (a1,a3 = 0 or 1, a2 = -1, 0 or 1).

* Return [N, [u,r,s,t], c], where

* N = arithmetic conductor of E

* c = product of the local Tamagawa numbers cp

* [u, r, s, t] = the change of variable reducing E to its minimal model,

* with u > 0 */

GEN

globalreduction(GEN E)

{

long k, l;

gpmem_t av = avma;

GEN c, prims, result, N, v, e;

v = ellintegralmodel(E);

e = ell_to_small(E);

if (v) e = coordch(e, v); else v = init_ch();

prims = (GEN)decomp(absi((GEN)e[12]))[1];

l = lg(prims); c = N = gun;

for (k = 1; k < l; k++)

{

GEN q = localreduction(e, (GEN)p1[k]);

GEN p = (GEN)prims[k];

GEN v1 = (GEN)q[3];

GEN q = localreduction(e, p);

N = mulii(N, gpui((GEN)p1[k],(GEN)q[1],0));

GEN w = (GEN)q[3];

N = mulii(N, powgi(p, (GEN)q[1]));

c = mulii(c, (GEN)q[4]);

if (!gcmp1((GEN)v1[1])) cumule1(&v, &e, v1);

if (!gcmp1((GEN)w[1]))

cumule(&v, &e, (GEN)w[1], (GEN)w[2], (GEN)w[3], (GEN)w[4]);

}

s = gdiventgs((GEN)e[1], -2);

standard_model(e, &v);

r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3);

result = cgetg(4, t_VEC);

t = gdiventgs(ellLHS0(e,r), -2);

result[1] = (long)N;

cumule(&v, &e, gun, r, s, t);

result[2] = (long)v;

tetpil = avma;

result[3] = (long)c; return gerepilecopy(av, result);

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

/* accumulate the effects of variable changes [u,r,s,t] and [U,R,S,T].

* particuliers frequents, tels que soit u = 1, soit r' = s' = t' = 0

* Frequent special cases: (U = 1) or (r = s = t = 0) */

static void

cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t)

{

long av = avma, tetpil;

GEN U,R,S,T,UU,RR,SS,TT, v = *vtotal, w = cgetg(5, t_VEC);

GEN temp, v = *vtotal, v3 = cgetg(5, t_VEC);

gpmem_t av = avma;

if (gcmp1((GEN)v[1]))

U = (GEN)v[1];

R = (GEN)v[2];

S = (GEN)v[3];

T = (GEN)v[4];

if (gcmp1(U))

{

v3[1] = lcopy(u);

UU = gcopy(u);

v3[2] = ladd((GEN)v[2], r);

RR = gadd(R, r);

v3[3] = ladd((GEN)v[3], s);

SS = gadd(S, s); av = avma;

av = avma;

TT = gerepileupto(av, gadd(T, gadd(t, gmul(S, r))));

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);

UU = gmul(U, u);

v3[2] = lcopy((GEN)v[2]);

RR = gcopy(R);

v3[3] = lcopy((GEN)v[3]);

SS = gcopy(S);

v3[4] = lcopy((GEN)v[4]);

TT = gcopy(T);

}

else /* cas general */

else /* general case */

{

v3[1] = lmul((GEN)v[1], u);

GEN U2 = gsqr(U);

temp = gsqr((GEN)v[1]);

UU = gmul(U, u);

v3[2] = ladd(gmul(temp, r), (GEN)v[2]);

RR = gadd(R, gmul(U2, r));

v3[3] = ladd(gmul((GEN)v[1], s), (GEN)v[3]);

SS = gadd(S, gmul(U, s));

v3[4] = ladd((GEN)v[4], gmul(temp, gadd(gmul((GEN)v[1], t), gmul((GEN)v[3], r))));

TT = gadd(T, gmul(U2, gadd(gmul(U, t), gmul(S, r))));

gerepileall(av, 4, &UU,&RR,&SS,&TT);

v3 = gerepilecopy(av, v3);

}

*vtotal = v3;

w[1] = (long)UU;

w[2] = (long)RR;

w[3] = (long)SS;

w[4] = (long)TT; *vtotal = w;

}

static void

cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t)

{

long av = avma, tetpil;

*e = _coordch(*e, u, r, s, t);

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

Line 2628 GEN

Line 2848 GEN

taniyama(GEN e)

{

GEN v,w,c,d,s1,s2,s3;

long n,m,av=avma,tetpil;

long n, m;

gpmem_t 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++)

for (n=-3; n<=(long)precdl-4; n++)

{

if (n!=2)

{

Line 2681 GEN

Line 2902 GEN

orderell(GEN e, GEN p)

{

GEN p1;

long av=avma,k;

long k;

gpmem_t av=avma;

checkell(e); checkpt(p);

k=typ(e[13]);

Line 2696 orderell(GEN e, GEN p)

Line 2918 orderell(GEN e, GEN 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];

switch(typ(a[i]))

{

case t_INT: case t_FRAC: case t_FRACN: break;

default: err(talker, "not a rational curve in ellintegralmodel");

}

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

Line 2771 GEN

Line 2972 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;

long i, j, nlr, t, t2, k, k2;

gpmem_t av=avma;

checkell(e);

v = ellintegralmodel(e);

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

pol = RHSpol(e);

Line 2854 torsellnagelllutz(GEN e)

Line 3055 torsellnagelllutz(GEN e)

/* Using Doud's algorithm */

/* Input e and n, finds a bound for #Tor */

/* finds a bound for #E_tor */

static long

torsbound(GEN e, long n)

torsbound(GEN e)

{

long av = avma, m, b, c, d, prime = 2;

long m, b, c, prime = 2;

gpmem_t av = avma;

byteptr p = diffptr;

GEN D = (GEN)e[12];

long n = bit_accuracy(lgefint(D)) >> 3;

/* n = number of primes to try ~ 1 prime every 8 bits in D */

b = c = m = 0; p++;

while (m<n)

{

d = *p++; if (!d) err(primer1);

NEXT_PRIME_VIADIFF_CHECK(prime,p);

prime += d;

if (smodis(D, prime))

if (ggval(D,stoi(prime)) == 0)

{

b = cgcd(b, prime+1 - itos(apell0(e,prime)));

if (b==c) m++; else {c = b; m = 0;}

Line 2886 _round(GEN x, long *e)

Line 3088 _round(GEN x, long *e)

return y;

}

/* Input the curve, a point, and an integer n, returns a point of order n

/* E the curve, w in C/Lambda ~ E of order n, returns q = pointell(w) as a

on the curve, or NULL if q is not rational. */

* rational point on the curve, or NULL if q is not rational. */

static GEN

torspnt(GEN E, GEN q, long n)

torspnt(GEN E, GEN w, long n, long prec)

{

GEN p = cgetg(3,t_VEC);

GEN p = cgetg(3,t_VEC), q = pointell(E, w, prec);

long e;

p[1] = lmul2n(_round(gmul2n((GEN)q[1],2), &e),-2);

if (e > -5) return NULL;

Line 2924 best_in_cycle(GEN e, GEN p, long k)

Line 3126 best_in_cycle(GEN e, GEN p, long k)

return (gsigne(d_ellLHS(e,p)) < 0)? invell(e,p): p;

}

/* <p,q> = E_tors, possibly NULL (= oo), p,q independant unless NULL

* order p = k, order q = 2 unless NULL */

static GEN

tors(GEN e, long k, GEN p, GEN q, GEN v)

{

Line 2939 tors(GEN e, long k, GEN p, GEN q, GEN v)

Line 3143 tors(GEN e, long k, GEN p, GEN q, GEN v)

p = best_in_cycle(e,p,k);

if (v)

{

v[1] = linv((GEN)v[1]);

p = pointch(p,v);

q = pointch(q,v);

}

Line 2954 tors(GEN e, long k, GEN p, GEN q, GEN v)

Line 3157 tors(GEN e, long k, GEN p, GEN q, GEN v)

{

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);

Line 2977 tors(GEN e, long k, GEN p, GEN q, GEN v)

Line 3177 tors(GEN e, long k, GEN p, GEN q, GEN v)

GEN

torselldoud(GEN e)

{

long b,i,ord,av=avma,prec, k = 1;

long b, i, ord, prec, k = 1;

GEN v,w1,w22,w1j,w12,p,tor1,tor2;

gpmem_t av=avma;

GEN v,w,w1,w22,w1j,w12,p,tor1,tor2;

checkbell(e);

v = ellintegralmodel(e);

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

b = lgefint((GEN)e[12]) >> 1; /* b = size of sqrt(D) */

b = DEFAULTPREC + ((lgefint((GEN)e[12])-2) >> 1); /* b >= size of sqrt(D) */

prec = precision((GEN)e[15]);

w1 = (GEN)e[15];

prec = precision(w1);

if (prec < b) err(precer, "torselldoud");

b = max(b, DEFAULTPREC);

if (b < prec) { prec = b; e = gprec_w(e, b); w1 = (GEN)e[15]; }

if (b < prec) { prec = b; e = gprec_w(e, b); }

b = torsbound(e);

b = torsbound(e,3);

if (b==1) { avma=av; return tors(e,1,NULL,NULL, v); }

if (v) v[1] = linv((GEN)v[1]);

w22 = gmul2n((GEN)e[16],-1);

w1 = (GEN)e[15];

if (b % 4)

{

{ /* cyclic of order 1, p, 2p, p <= 5 */

p = NULL;

for (i=10; i>1; i--)

{

if (b%i != 0) continue;

w1j = gdivgs(w1,i);

p = torspnt(e,pointell(e,w1j,prec),i);

p = torspnt(e,w1j,i,prec);

if (!p && i%2==0)

{

p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);

p = torspnt(e,gadd(w22,w1j),i,prec);

if (!p) p = torspnt(e,pointell(e,gadd(w22,gmul2n(w1j,1)),prec),i);

if (!p) p = torspnt(e,gadd(w22,gmul2n(w1j,1)),i,prec);

}

if (p) { k = i; break; }

}

Line 3012 torselldoud(GEN e)

Line 3213 torselldoud(GEN e)

}

ord = 0; tor1 = tor2 = NULL;

w12 = gmul2n((GEN)e[15],-1);

w12 = gmul2n(w1,-1);

if ((p = torspnt(e,pointell(e,w12,prec),2)))

if ((p = torspnt(e,w12,2,prec)))

{

tor1 = p; ord++;

}

if ((p = torspnt(e,pointell(e,w22,prec),2))

w = w22;

|| (!ord && (p = torspnt(e,pointell(e,gadd(w12,w22),prec),2))))

if ((p = torspnt(e,w,2,prec)))

{

tor2 = p; ord += 2;

}

if (!ord)

{

w = gadd(w12,w22);

if ((p = torspnt(e,w,2,prec)))

{

tor2 = p; ord += 2;

}

p = NULL;

switch(ord)

{

case 0:

case 0: /* no point of order 2 */

for (i=9; i>1; i-=2)

{

if (b%i!=0) continue;

w1j=gdivgs((GEN)e[15],i);

w1j = gdivgs(w1,i);

p = torspnt(e,pointell(e,w1j,prec),i);

p = torspnt(e,w1j,i,prec);

if (p) { k = i; break; }

}

break;

case 1:

case 1: /* 1 point of order 2: w1 / 2 */

p = NULL;

for (i=12; i>2; i-=2)

{

if (b%i!=0) continue;

w1j=gdivgs((GEN)e[15],i);

w1j = gdivgs(w1,i);

p = torspnt(e,pointell(e,w1j,prec),i);

p = torspnt(e,w1j,i,prec);

if (!p && i%4==0)

p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);

p = torspnt(e,gadd(w22,w1j),i,prec);

if (p) { k = i; break; }

}

if (!p) { p = tor1; k = 2; }

break;

case 2:

case 2: /* 1 point of order 2: w = w2/2 or (w1+w2)/2 */

for (i=5; i>1; i-=2)

{

if (b%i!=0) continue;

w1j = gdivgs((GEN)e[15],i);

w1j = gdivgs(w1,i);

p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i+i);

p = torspnt(e,gadd(w,w1j),2*i,prec);

if (p) { k = 2*i; break; }

}

if (!p) { p = tor2; k = 2; }

tor2 = NULL; break;

case 3:

case 3: /* 2 points of order 2: w1/2 and w2/2 */

for (i=8; i>2; i-=2)

{

if (b%(2*i)!=0) continue;

w1j=gdivgs((GEN)e[15],i);

w1j = gdivgs(w1,i);

p = torspnt(e,pointell(e,w1j,prec),i);

p = torspnt(e,w1j,i,prec);

if (p) { k = i; break; }

}

if (!p) { p = tor1; k = 2; }

Line 3089 elltors0(GEN e, long flag)

Line 3297 elltors0(GEN e, long flag)

/* par compatibilite */

GEN torsell(GEN e) {return torselldoud(e);}

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

/* LOCAL ROOT NUMBERS, after HALBERSTADT */

/* ici p=2 ou 3 */

/* p = 2 or 3 */

static long

neron(GEN e, GEN p, long* ptkod)

{

long av=avma,kod,v4,v6,vd;

long kod, v4, v6, vd;

gpmem_t av=avma;

GEN c4, c6, d, nv;

nv=localreduction(e,p);

Line 3196 neron(GEN e, GEN p, long* ptkod)

Line 3405 neron(GEN e, GEN p, long* ptkod)

static long

ellrootno_2(GEN e)

{

long n2,kod,u,v,x1,y1,d1,av=avma,v4,v6,w2;

long n2, kod, u, v, x1, y1, d1, v4, v6, w2;

gpmem_t av=avma;

GEN p=gdeux,c4,c6,tmp,p6;

n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p6=stoi(64);

Line 3283 ellrootno_2(GEN e)

Line 3493 ellrootno_2(GEN e)

static long

ellrootno_3(GEN e)

{

long n2,kod,u,v,d1,av=avma,r6,K4,K6,v4;

long n2, kod, u, v, d1, r6, K4, K6, v4;

gpmem_t av=avma;

GEN p=stoi(3),c4,c6,tmp,p4;

n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p4=stoi(81);

Line 3328 ellrootno_3(GEN e)

Line 3539 ellrootno_3(GEN e)

}

/* assume p > 3 */

static long

ellrootno_not23(GEN e, GEN p, GEN ex)

ellrootno_p(GEN e, GEN p, GEN ex)

{

GEN j;

long ep,z;

Line 3346 ellrootno_not23(GEN e, GEN p, GEN ex)

Line 3558 ellrootno_not23(GEN e, GEN p, GEN ex)

static long

ellrootno_intern(GEN e, GEN p, GEN ex)

{

if (cmpis(p,3) > 0) return ellrootno_not23(e,p,ex);

if (cmpis(p,3) > 0) return ellrootno_p(e,p,ex);

switch(itos(p))

{

case 3: return ellrootno_3(e);

Line 3378 ellrootno_all(GEN e, GEN p, GEN* ptcond)

Line 3590 ellrootno_all(GEN e, GEN p, GEN* ptcond)

exs=ggval(cond,p);

if (!exs) return 1;

}

if (cmpis(p,3)>0) return ellrootno_not23(e,p,stoi(exs));

if (cmpis(p,3)>0) return ellrootno_p(e,p,stoi(exs));

switch(itos(p))

{

case 3: return ellrootno_3(e);

Line 3393 ellrootno_all(GEN e, GEN p, GEN* ptcond)

Line 3605 ellrootno_all(GEN e, GEN p, GEN* ptcond)

long

ellrootno(GEN e, GEN p)

{

long av=avma,s;

long s;

gpmem_t av=avma;

if (!p) p = gun;

s=ellrootno_all(e, p, NULL);

avma=av; return s;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>