| Karim BELABAS on Thu, 3 Jun 1999 16:19:44 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| elltors |
[Optional patch to 2.0.15]
Following John Cremona's suggestion, the following replaces the old elltors
[still available with elltors(e, 1)] by an implementation of Darrin Doud's
algorithm to find torsion points on an elliptic curve defined over Q. It
should perform orders of magnitude faster than the old function, but will
now require a "big elliptic curve" (initell(..., 1) won't do). On the other
hand, it shouldn't require an integral model.
The patch is huge, and probably useless to you if you don't intend to test
the elltors function. It cleans up a lot of code, but it looks very useful,
and I think that it can be included at this stage without endangering the
release. Give it a rough treatment...
Karim.
*** src/modules/elliptic.c.orig Thu Jun 3 15:54:33 1999
--- src/modules/elliptic.c Mon May 31 21:11:13 1999
***************
*** 40,45 ****
--- 40,96 ----
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)
{
***************
*** 101,115 ****
}
static GEN
! do_padic_agm(GEN *ptx1, GEN p, GEN a1, GEN b1)
{
GEN p1,r1,a,b,x,bmod1, bmod = modii((GEN)b1[4],p), x1 = *ptx1;
for(;;)
{
a=a1; b=b1; x=x1;
b1=gsqrt(gmul(a,b),0); bmod1=modii((GEN)b1[4],p);
! if (!gegal(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);
--- 152,188 ----
}
static GEN
! do_agm(GEN *ptx1, GEN a1, GEN b1, long prec, long sw)
! {
! long G = 6 - bit_accuracy(prec);
! GEN p1,r1,a,b,x,x1;
!
! x1 = gmul2n(gsub(a1,b1),-2);
! if (gcmp0(x1)) err(talker,"precision too low in initell");
! 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;
! }
! *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);
***************
*** 120,133 ****
*ptx1 = x1; return ginv(gmul2n(a1,2));
}
! GEN
! initell(GEN x, long prec)
{
! GEN y,b2,b4,disc,c4,c6,p1,p2,p,u,w,pv;
! GEN aa1,bb1,r1,x1,u2,q,e0,e1;
! long ty,i,av,tetpil,e,alpha;
! av=avma; y=cgetg(20,t_VEC);
smallinitell0(x,y);
e=BIGINT;
--- 193,260 ----
*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,c4,c6,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;
***************
*** 140,203 ****
p = (GEN)q[2];
}
}
! if (e<BIGINT)
! {
! q=ggrandocp(p,e);
! for (i=1; i<=5; i++) y[i]=ladd(q,(GEN)x[i]);
! }
b2 = (GEN) y[6];
b4 = (GEN) y[7];
c4 = (GEN) y[10];
c6 = (GEN) y[11];
! disc = (GEN) y[12];
!
! ty=typ(disc);
! if (ty != t_PADIC)
! {
! if (prec && is_const_t(ty) && ty!=t_INTMOD)
! {
! GEN pi=mppi(prec), pi2=gmul2n(pi,1), tau, w2;
! long sw;
!
! p1=cgetg(6,t_POL); p1[1]=evalsigne(1) | evallgef(6);
! p1[2]=y[8]; /* b6 */
! p1[3]=lmul2n(b4,1);
! p1[4]=y[6]; /* b2 */
! p1[5]=lstoi(4);
! p1=roots(p1,prec);
! if (gsigne(disc) < 0) p1[1]=lreal((GEN)p1[1]);
! else
! { /* sort roots, real parts in decreasing order */
! GEN tmp;
! i=1; p1 = greal(p1); tmp = (GEN)p1[1];
! if (gcmp((GEN)p1[2],tmp)>0) { i=2; tmp=(GEN)p1[2]; }
! if (gcmp((GEN)p1[3],tmp)>0) { i=3; tmp=(GEN)p1[3]; }
! p1[i]=p1[1]; p1[1]=(long)tmp;
!
! if (gcmp((GEN)p1[2],(GEN)p1[3])<0)
! { tmp=(GEN)p1[3]; p1[3]=p1[2]; p1[2]=(long)tmp; }
! }
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);
! sw=signe(w);
!
! aa1=gmul2n(gsub(w,p2),-2);
! bb1=gmul2n(w,-1); r1=gsub(aa1,bb1); x1=gmul2n(r1,-2);
! if (gcmp0(r1)) err(talker,"precision too low in initell");
! do
! {
! GEN aa0=aa1, bb0=bb1, x0=x1;
! bb1=gsqrt(gmul(aa0,bb0),prec); setsigne(bb1,sw);
! aa1=gmul2n(gadd(gadd(aa0,bb0),gmul2n(bb1,1)),-2);
! r1=gsub(aa1,bb1);
! x1=gmul(x0,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(x0,r1),x0),prec)),-1)));
! }
! while (gexpo(r1) > gexpo(bb1) - bit_accuracy(prec) + 6);
! u2=ginv(gmul2n(aa1,2));
w = gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
q = gsqrt(gaddgs(gsqr(w),-1),prec);
--- 267,296 ----
p = (GEN)q[2];
}
}
! if (e<BIGINT) return padic_initell(y,p,e);
!
b2= (GEN)y[6];
b4= (GEN)y[7];
c4= (GEN)y[10];
c6= (GEN)y[11];
! 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);
***************
*** 210,281 ****
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;
else
{
y[15]=lmul2n(gabs((GEN)w2[1],prec),1);
! q=gexp(gmul(gmul(pi2,gi),gdiv(w2,(GEN)y[15])),prec);
}
y[16]=(long)w2;
! p1=gdiv(gsqr(pi),gmulsg(6,(GEN)y[15])); q=gsqrt(q,prec);
y[17]=lmul(p1,gdiv(thetanullk(q,3,prec),thetanullk(q,1,prec)));
y[18]=ldiv(gsub(gmul((GEN)y[17],(GEN)y[16]),gmul(gi,pi)),(GEN)y[15]);
! }
! else { y[14]=y[15]=y[16]=y[17]=y[18]=y[19]=zero; }
! tetpil=avma; return gerepile(av,tetpil,gcopy(y));
! }
!
! if (valp(y[13]) >= 0) /* p | j */
! err(talker,"valuation of j must be negative in p-adic ellinit");
!
! p=(GEN)disc[2];
! if (cmpis(p,2)) pv=p;
! else
! {
! pv=stoi(4); err(impl,"initell for 2-adic numbers");
! }
! 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;
!
! aa1=gmul2n(gsub(w,gadd(gmulsg(3,e1),gmul2n(b2,-2))),-2);
! bb1=gmul2n(w,-1); r1=gsub(aa1,bb1); x1=gmul2n(r1,-2);
! u2 = do_padic_agm(&x1,pv,aa1,bb1);
!
! 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);
!
! y[15]=(long)u2;
! y[16] = (kronecker((GEN)u2[4],p) > 0 && (valp(u2)&1) == 0)?
! lsqrt(u2,0): zero;
! y[17]=(long)q; p1=cgetg(2,t_VEC);
! y[14]=(long)p1; p1[1]=(long)e1; y[19]=zero;
! tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}
! static GEN
! ellpol_eval(GEN e, GEN r)
{
! GEN p1;
! p1 = gadd((GEN)e[2],r);
! p1 = gadd((GEN)e[4], gmul(r,p1));
! p1 = gadd((GEN)e[5], gmul(r,p1));
! return p1;
}
GEN
--- 303,330 ----
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]));
y[17]=lmul(p1,gdiv(thetanullk(q,3,prec),thetanullk(q,1,prec)));
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
***************
*** 290,302 ****
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))));
! p1 = gadd(gadd(gmul(r,(GEN)e[1]),gmul2n(t,1)),(GEN)e[3]);
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)))));
! p1 = ellpol_eval(e,r);
! p2 = gmul(t,gadd(gadd(gmul(r,(GEN)e[1]),t),(GEN)e[3]));
! y[5] = lmul(v6,gsub(p1,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)))))));
--- 339,351 ----
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)))))));
***************
*** 321,327 ****
y[14]=(long)p2;
y[15]=lmul(gsqr(u),(GEN)e[15]);
y[16]=lmul(u,(GEN)e[16]);
! /* A MODIFIER : comment changent q et w ??? */
y[17]=e[17];
y[18]=e[18];
y[19]=zero;
--- 370,376 ----
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;
***************
*** 391,403 ****
int
oncurve(GEN e, GEN z)
{
! GEN p1,p2,x,y;
long av=avma,p,q;
checksell(e); checkpt(z); if (lg(z)<3) return 1; /* oo */
! x=(GEN)z[1]; y=(GEN)z[2];
! p1 = gmul(y,gadd(y,gadd((GEN)e[3],gmul((GEN)e[1],x))));
! p2 = ellpol_eval(e,x); x = gsub(p1,p2);
if (gcmp0(x)) { avma=av; return 1; }
p = precision(p1);
q = precision(p2); if (p < q) q = p;
--- 440,451 ----
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) q = p;
***************
*** 419,436 ****
x1=(GEN)z1[1]; y1=(GEN)z1[2];
x2=(GEN)z2[1]; y2=(GEN)z2[2];
! if (x1 == x2 || gegal(x1,x2)) /* implies y1 = y2 or -y2 */
! {
if (y1 != y2)
{
! long equal;
if (precision(y1) || precision(y2))
! equal = (gexpo(gadd(y1,y2)) >= gexpo(y1));
else
! equal = gegal(y1,y2);
! if (!equal) { avma=av; y=cgetg(2,t_VEC); y[1]=zero; return y; }
}
! p2=gadd((GEN)e[3],gadd(gmul((GEN)e[1],x1),gmul2n(y1,1)));
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))));
--- 467,484 ----
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))));
***************
*** 438,444 ****
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(gadd((GEN)e[3],y1),gmul(x,(GEN)e[1])),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);
}
--- 486,492 ----
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);
}
***************
*** 450,456 ****
if (lg(z)<3) return z;
p1=cgetg(3,t_VEC); p1[1]=z[1];
! p1[2]=(long)gneg_i(gadd(gadd((GEN)z[2],(GEN)e[3]),gmul((GEN)e[1],(GEN)z[1])));
return p1;
}
--- 498,504 ----
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;
}
***************
*** 460,502 ****
long av=avma,tetpil;
checksell(e); checkpt(z2);
-
z2=invell(e,z2); tetpil=avma;
return gerepile(av,tetpil,addell(e,z1,z2));
}
- static GEN
- quad_sol(long av, GEN d, GEN b, GEN sqrtd)
- {
- GEN y, p1 = gsub(sqrtd,b);
- long tetpil = avma;
-
- if (! signe(d))
- {
- y=cgetg(2,t_VEC);
- y[1]=lmul2n(p1,-1);
- }
- else
- {
- y=cgetg(3,t_VEC);
- y[1]=lmul2n(p1,-1);
- y[2]=lsub((GEN)y[1],sqrtd);
- }
- return gerepile(av,tetpil,y);
- }
-
- static GEN
- cgetimod(GEN x, GEN y)
- {
- GEN p1 = cgetg(3,t_INTMOD);
- p1[1]=(long)x; p1[2]=(long)y; return p1;
- }
-
GEN
ordell(GEN e, GEN x, long prec)
{
! GEN p1,p2,p3,p4,p5,d,y,pd;
! long av=avma,tetpil,td,i,lx,tx=typ(x);
checksell(e);
if (is_matvec_t(tx))
--- 508,522 ----
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))
***************
*** 506,586 ****
return y;
}
! p1 = ellpol_eval(e,x);
! p2 = gadd((GEN)e[3],gmul((GEN)e[1],x));
! d = gadd(gsqr(p2),gmul2n(p1,2));
! td=typ(d);
! if (td==t_INT)
! {
! if (!carrecomplet(d,&p3))
! { avma=av; return cgetg(1,t_VEC); }
! return quad_sol(av,d,p2, p3);
! }
! if (is_frac_t(td))
! {
! pd=(GEN)d[2]; d=(GEN)d[1];
! if (!carrecomplet(mulii(d,pd),&p3))
! { avma=av; return cgetg(1,t_VEC); }
! return quad_sol(av,d,p2, gdiv(p3,pd));
! }
!
! if (td==t_INTMOD)
! {
! if (! cmpis((GEN)d[1],2)) /* modulo 2 */
! {
! avma=av;
! if (gcmp0(p2))
{
! y=cgetg(2,t_VEC); p1 = gcmp0(p1)? gzero: gun;
! y[1] = (long) cgetimod(gdeux,p1);
! return y;
! }
! if (gcmp0(p1))
! {
! y=cgetg(3,t_VEC);
! y[1] = (long) cgetimod(gdeux,gzero);
! y[2] = (long) cgetimod(gdeux,gun);
! return y;
! }
! return cgetg(1,t_VEC);
! }
! if (kronecker((GEN)d[2],(GEN)d[1]) == -1)
! {
! avma=av; y=cgetg(1,t_VEC); return y;
! }
}
! p3=gsqrt(d,prec); p5=gsub(p3,p2); tetpil=avma;
! if (gcmp0(d))
{
! y=cgetg(2,t_VEC);
! y[1]=lmul2n(p5,-1);
}
else
{
! if (td==t_REAL || td==t_COMPLEX)
{
! p1=gneg_i(p1); p4=gadd(p3,p2);
! if (gcmp(gnorm(p5),gnorm(p4)) < 0) p5 = gneg_i(p4);
! tetpil=avma;
! y=cgetg(3,t_VEC);
! y[1]=lmul2n(p5,-1);
! y[2]=ldiv(p1,(GEN)y[1]);
! }
! else
{
y=cgetg(3,t_VEC);
! y[1]=lmul2n(p5,-1);
! y[2]=lsub((GEN)y[1],p3);
}
}
! return gerepile(av,tetpil,y);
}
static GEN
CM_powell(GEN e, GEN z, GEN n)
{
! GEN y,p0,p1,q0,q1,p2,q2,z1,z2,pol,grdx,resx,resy;
long av=avma,tetpil,ln,ep,vn;
if (lg(z)<3) return gcopy(z);
--- 526,579 ----
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);
***************
*** 597,603 ****
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;
--- 590,597 ----
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;
***************
*** 615,625 ****
while (lgef(p1)-3 < vn);
if (lgef(p1)-3 > vn || signe(z2))
err(talker,"not a complex multiplication in powell");
! resx=gdiv(p1,q1); resy=gdiv(deriv(resx,0),n);
! resx=gsub(poleval(resx,grdx), gdivgs((GEN)e[6],12));
! resy=gmul2n(gsub(gmul(gadd(gadd(gmul2n((GEN)z[2],1),gmul((GEN)e[1],(GEN)z[1])),(GEN)e[3]),poleval(resy,grdx)),gadd(gmul((GEN)e[1],resx),(GEN)e[3])),-1);
! tetpil=avma; y=cgetg(3,t_VEC); y[1]=lcopy(resx); y[2]=lcopy(resy);
! return gerepile(av,tetpil,y);
}
GEN
--- 609,619 ----
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
***************
*** 758,777 ****
zell(GEN e, GEN z, long prec)
{
long av=avma,ty,sw,fl;
! GEN t,u,p1,p2,r1,a,b,x1,p,u2;
checkbell(e);
if (!oncurve(e,z)) err(heller1);
! ty=typ(e[12]);
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)
{
! p=gmael(e,12,2);
! u2 = do_padic_agm(&x1,p,a,b);
!
if (!gcmp0((GEN)e[16]))
{
t=gsqrt(gaddsg(1,gdiv(x1,a)),prec);
--- 752,769 ----
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);
***************
*** 807,813 ****
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;
--- 799,805 ----
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;
***************
*** 826,832 ****
fprintferr(" z = ");outerr(z);
fprintferr(" u = ");outerr(u);
fprintferr(" x1 = ");outerr(x1);
- fprintferr(" bad= %ld\n",bad);
}
fprintferr("ellpointtoz: `sqrt' gave the %s square root\n",s);
flusherr();
--- 818,823 ----
***************
*** 1269,1275 ****
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(gadd((GEN)e[3],gmul((GEN)e[1],y)),-1));
tetpil=avma; v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lcopy(yp);
return gerepile(av,tetpil,v);
}
--- 1260,1266 ----
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);
}
***************
*** 1914,1927 ****
hell(GEN e, GEN a, long prec)
{
long av=avma,tetpil,n;
! GEN p1,p2,y,z,q,psi2,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));
- psi2=gadd((GEN)e[3],gadd(gmul((GEN)e[1],(GEN)a[1]),gmul2n((GEN)a[2],1)));
y=gsin(z,prec); n=0; qn=gun; ps=gneg_i(q);
do
{
--- 1905,1917 ----
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
{
***************
*** 1929,1935 ****
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),psi2)),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);
--- 1919,1925 ----
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);
***************
*** 2036,2045 ****
if (lg(a)<3) return gzero;
if (!oncurve(e,a)) err(heller1);
! x=(GEN)a[1]; y=(GEN)a[2];
! psi2=numer(gadd((GEN)e[3],gadd(gmul((GEN)e[1],x),gmul2n(y,1))));
if (!signe(psi2)) { avma=av; return gzero; }
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; }
--- 2026,2035 ----
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; }
***************
*** 2231,2237 ****
case 1: r = divis(subii(p2k, a2prime), 3); break;
case 2: r = negi(divis(addii(a2prime, p2k), 3)); break;
}
! a3prime = addii((GEN)e[3], mulii(r, (GEN)e[1]));
if (mpodd(a3prime))
t = shifti(subii(mulii(pk, p2k), a3prime), -1);
else
--- 2221,2227 ----
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
***************
*** 2370,2376 ****
{
r = negi(modis((GEN)e[8], 3));
s = modis((GEN)e[1], 3);
! t = modis(addii((GEN)e[3], mulii((GEN)e[1], r)), 3);
}
cumule(&v, &e, gun, r, s, t); /* p | a1, a2, a3, a4 et a6 */
p2 = stoi(p*p);
--- 2360,2366 ----
{
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);
***************
*** 2561,2567 ****
}
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(gadd((GEN)e[3], gmul(r,(GEN)e[1])), -2);
cumule(&v, &e, gun, r, s, t);
tetpil = avma;
result = cgetg(4, t_VEC); result[1] = lcopy(N); result[2] = lcopy(v);
--- 2551,2557 ----
}
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);
***************
*** 2674,2689 ****
v[n+6]=lneg(gdiv((GEN)s2[2],(GEN)s2[3]));
}
}
! w=gsub(gmul(polx[0],gmul(d,deriv(v,0))),gcmp0((GEN)e[1])?(GEN)e[3]:gadd((GEN)e[3],gmul((GEN)e[1],v)));
tetpil=avma; s1=cgetg(3,t_VEC); s1[1]=lcopy(v); s1[2]=lmul2n(w,-1);
return gerepile(av,tetpil,s1);
}
/********************************************************************/
/** **/
! /** Points de torsion **/
/** **/
/********************************************************************/
/* p is a polynomial of degree exactly 3 with integral coefficients
* and leading term 4. Outputs the vector of rational roots of p
--- 2664,2700 ----
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;
+ }
+
+ /* 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
***************
*** 2700,2741 ****
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=0;
! 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(gsubst(p,0,a))) v[++t]=(long)a;
a = gneg_i(a);
! if (!gsigne(gsubst(p,0,a))) v[++t]=(long)a;
! }
! setlg(v,t+1); return v;
! }
!
! /* we assume e is defined over Q */
! 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;
}
static int
! is_new_torsion(GEN e, GEN v, GEN p, long t2)
! {
GEN pk = p, pkprec = NULL;
long k,l;
--- 2711,2731 ----
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;
***************
*** 2747,2828 ****
for (l=2; l<=t2; l++)
if (gegal((GEN)pk[1],gmael(v,l,1))) return 1;
! if (pkprec && lg(pkprec)>2)
! if (k<=5 && gegal((GEN)pk[1],(GEN)pkprec[1])) return 1;
pkprec=pk;
}
return 0;
}
GEN
! torsell(GEN e)
{
! GEN n,d,ld,pol,p1,lr,v,w;
! long i,j,nlr,t,t2,k,k2,av=avma,tetpil,av1;
! checkell(e); t=1;
! v=cgetg(17,t_VEC); p1=cgetg(2,t_VEC); p1[1]=zero; v[1]=(long)p1;
! pol=gadd((GEN)e[8],gmul(polx[0],gadd(gmul2n((GEN)e[7],1),gmul(polx[0],gadd((GEN)e[6],gmul2n(polx[0],2))))));
lr=ratroot(pol); nlr=lg(lr)-1;
! for (i=1; i<=nlr; i++)
{
! p1=cgetg(3,t_VEC); p1[1]=lr[i];
! p1[2]=ldivgs(gadd(gmul((GEN)e[1],(GEN)lr[i]),(GEN)e[3]),-2);
v[++t]=(long)p1;
}
! t2=t;
! ld=factor(gabs((GEN)e[12],0)); n=stoi(4);
! for (i=1; i<lg(ld[1]); i++)
! n = gmul(n,gpui(gcoeff(ld,i,1),shifti(gcoeff(ld,i,2),-1),0));
! ld=divisors(n);
! for (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,gadd(gmul((GEN)e[1],(GEN)lr[i]),(GEN)e[3])),-1);
if (is_new_torsion(e,v,p1,t2))
{
GEN p2 = cgetg(3,t_VEC);
! p2[1]=p1[1]; p2[2]=lsub((GEN)p1[2],d);
! v[++t]=(long)p1; v[++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;
}
- tetpil=avma; w=cgetg(4,t_VEC); w[1]=lstoi(t);
if (nlr<3)
{
! p1=cgetg(2,t_VEC); p1[1]=lstoi(t); w[2]=(long)p1;
! av1=avma;
! k=2; while (k<=t && itos(orderell(e,(GEN)v[k])) != t) k++;
! avma=av1;
if (k>t) err(bugparier,"torsell (bug1)");
! p1=cgetg(2,t_VEC); p1[1]=lcopy((GEN)v[k]); w[3]=(long)p1;
! return gerepile(av,tetpil,w);
}
if (t&3) err(bugparier,"torsell (bug2)");
!
! p1=cgetg(3,t_VEC); t2=t>>1; p1[1]=lstoi(t2); p1[2]=(long)gdeux;
! w[2]=(long)p1; av1=avma;
! k=2; while (k<=t && itos(orderell(e,(GEN)v[k])) != t2) k++;
if (k>t) err(bugparier,"torsell (bug3)");
! k2=2; p1=powell(e,(GEN)v[k],stoi(t>>2));
! if (lg(p1)==3 && gegal((GEN)v[2],p1)) k2++;
! avma=av1;
! p1=cgetg(3,t_VEC); p1[1]=lcopy((GEN)v[k]); p1[2]=lcopy((GEN)v[k2]);
! w[3]=(long)p1; return gerepile(av,tetpil,w);
}
/* LOCAL ROOT NUMBERS, D'APRES HALBERSTADT halberst@math.jussieu.fr */
--- 2737,3031 ----
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,v,w,w2,w3;
! long i,j,nlr,t,t2,k,k2,av=avma;
! checkell(e);
! for (i=1; i<=5; i++)
! if (typ(e[i]) != t_INT) err(talker, "not an integral model in torsell");
! pol = RHSpol(e);
lr=ratroot(pol); nlr=lg(lr)-1;
! v=cgetg(17,t_VEC); p1=cgetg(2,t_VEC); p1[1]=zero; v[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);
v[++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,v,p1,t2))
{
GEN p2 = cgetg(3,t_VEC);
! p2[1] = p1[1];
! p2[2] = lsub((GEN)p1[2],d);
! v[++t]=(long)p1;
! v[++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)v[k])) == t) break;
if (k>t) err(bugparier,"torsell (bug1)");
!
! w3=cgetg(2,t_VEC); w3[1]=v[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)v[k])) == t2) break;
if (k>t) err(bugparier,"torsell (bug3)");
! p1 = powell(e,(GEN)v[k],stoi(t>>2));
! k2 = (lg(p1)==3 && gegal((GEN)v[2],p1))? 3: 2;
! w3=cgetg(3,t_VEC); w3[1]=v[k]; w3[2]=v[k2];
}
+ 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;
+ }
+
+ /* Input the curve, a point, and an integer n, returns a point of order n
+ on the curve, or 0 if the point p is not a proper point. */
+ static GEN
+ torspnt(GEN e, GEN q, long n)
+ {
+ GEN p = cgetg(3,t_VEC);
+ p[1] = lmul2n(ground(gmul2n((GEN)q[1],2)),-2);
+ p[2] = lmul2n(ground(gmul2n((GEN)q[2],3)),-3);
+ return (oncurve(e,p) && gcmp0(gimag(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 p1, r = cgetg(4,t_VEC);
+ 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);
+
+ 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);
+ 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[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 w1,w22,w1j,w12,p,tor1,tor2;
+
+ checkbell(e);
+
+ prec=max(MEDDEFAULTPREC,gprecision(e));
+ b = torsbound(e,3);
+ if (b==1) { avma=av; return tors(e,1,NULL,NULL); }
+ 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));
+ }
+
+ 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));
+ }
+
+ 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 */
*** src/language/init.c.orig Thu Jun 3 16:10:55 1999
--- src/language/init.c Wed May 26 21:11:40 1999
***************
*** 1374,1380 ****
{"ellsigma",99,(void*)ellsigma,5,"GGD0,L,p"},
{"ellsub",99,(void*)subell,5,"GGGp"},
{"elltaniyama",1,(void*)taniyama,5,"Gp"},
! {"elltors",1,(void*)torsell,5,"Gp"},
{"ellwp",99,(void*)ellwp0,5,"GDGD0,L,pP"},
{"ellzeta",99,(void*)ellzeta,5,"GGp"},
{"ellztopoint",29,(void*)pointell,5,"GGp"},
--- 1374,1380 ----
{"ellsigma",99,(void*)ellsigma,5,"GGD0,L,p"},
{"ellsub",99,(void*)subell,5,"GGGp"},
{"elltaniyama",1,(void*)taniyama,5,"Gp"},
! {"elltors",99,(void*)elltors0,5,"GD0,L,p"},
{"ellwp",99,(void*)ellwp0,5,"GDGD0,L,pP"},
{"ellzeta",99,(void*)ellzeta,5,"GGp"},
{"ellztopoint",29,(void*)pointell,5,"GGp"},
*** src/language/helpmsg.c.orig Thu Jun 3 16:11:34 1999
--- src/language/helpmsg.c Wed May 26 21:13:39 1999
***************
*** 121,127 ****
"ellsigma(om,z,{flag=0}): om=[om1,om2], value of the Weierstrass sigma function of the lattice generated by om at z if flag = 0 (default). If flag = 1, arbitrary determination of the logarithm of sigma. If flag = 2 or 3, same but using the product expansion instead of theta series",
"ellsub(e,z1,z2): difference of the points z1 and z2 on elliptic curve e",
"elltaniyama(e): modular parametrization of elliptic curve e",
! "elltors(e): torsion subgroup of elliptic curve e: order, structure, generators",
"ellwp(e,{z=x},{flag=0}): Complex value of Weierstrass P function at z on the lattice generated over Z by e=[om1,om2] (e as given by ellinit is also accepted). Optional flag means 0 (default), compute only P(z), 1 compute [P(z),P'(z)], 2 consider om as an elliptic curve and compute P(z) for that curve (identical to ellztopoint in that case). If z is omitted or is a simple variable, return formal expansion in z",
"ellzeta(om,z): om=[om1,om2], value of the Weierstrass zeta function of the lattice generated by om at z",
"ellztopoint(e,z): coordinates of point P on the curve e corresponding to the complex number z",
--- 121,127 ----
"ellsigma(om,z,{flag=0}): om=[om1,om2], value of the Weierstrass sigma function of the lattice generated by om at z if flag = 0 (default). If flag = 1, arbitrary determination of the logarithm of sigma. If flag = 2 or 3, same but using the product expansion instead of theta series",
"ellsub(e,z1,z2): difference of the points z1 and z2 on elliptic curve e",
"elltaniyama(e): modular parametrization of elliptic curve e",
! "elltors(e,{flag=0}): torsion subgroup of elliptic curve e: order, structure, generators. If flag = 0, use Doud's algorithm; if flag = 1, use Lutz-Nagell",
"ellwp(e,{z=x},{flag=0}): Complex value of Weierstrass P function at z on the lattice generated over Z by e=[om1,om2] (e as given by ellinit is also accepted). Optional flag means 0 (default), compute only P(z), 1 compute [P(z),P'(z)], 2 consider om as an elliptic curve and compute P(z) for that curve (identical to ellztopoint in that case). If z is omitted or is a simple variable, return formal expansion in z",
"ellzeta(om,z): om=[om1,om2], value of the Weierstrass zeta function of the lattice generated by om at z",
"ellztopoint(e,z): coordinates of point P on the curve e corresponding to the complex number z",
__
Karim Belabas email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France) Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://hasse.mathematik.tu-muenchen.de/ntsw/pari/