File: [local] / OpenXM_contrib / pari / src / modules / Attic / elliptic.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:33 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/********************************************************************/
/** **/
/** 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;
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to elliptic.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / modules