Karim BELABAS on Tue, 19 Oct 1999 15:20:16 +0200 (MET DST)

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

 Re: ellpointtoz() bug?

```> gp> E=ellinit([0,0,1,-1,0]);
> gp> for(n=-6,6,print(n" "ellpow(E,[0,0],n)" "ellpointtoz(E,ellpow(E,[0,0],n))))
[...]
> -1 [0, -1] 2.063865930946563955426810634 + 1.225694690993395030427112415*I
> 0  0
> 1 [0, 0] 2.063865930946563955426810634 + 1.225694690993395030427112415*I
[...] which is clearly nonsense. What happens is that ellpointtoz (zell
internally) often "mistakes" a point for its opposite [notice that
[0,0] = - [0,-1]].

In fact the y coordinate is not used at all in the main algorithm [ it used
to be in 1.39, but mistakenly, since it endeavoured to give a wrong answer
from times to times. A comment acknowledged this in the code. ]

The chosen "solution" (implemented in the function as it stood in 2.0.17) was
to compute the inverse function (ellztopoint = pointell internally) at the
lowest possible precision [in fact just \wp'(z)] and set z := -z if it didn't
give back the original point.

There are two problems with this:

1) it'd be much nicer to have a nice theoretical test [if you look at the
code, the ambiguity is due to a (complex) square root taken at the end.
Changing determination only involves periods, hence is OK, but I have no idea
how to choose the sign (I'm no expert on elliptic curves...)]

2) it doesn't generalize well. E.g ellztopoint is not implemented for p-adic
fields; as it stands, ellztopoint over Q_p will always yield the same value
when taken at P or -P !

There was a third problem but it was a bug in the implementation of my
test... This is corrected by the following patch [ which formally corrects
the behaviour in your example above ].

Cheers,

Karim.

P.S: I'd be very glad to have an expert opinion on this [ and even happier
with an actual patch ! ]

P.S2: the patch is too big, but the CVS version evolved quite a bit since the
2.0.17 release (I already deleted many irrelevant chunks).

Index: src/modules/elliptic.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/modules/elliptic.c,v
retrieving revision 1.2
retrieving revision 1.10
diff -c -r1.2 -r1.10
***************
*** 778,784 ****
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;

checkbell(e);
if (!oncurve(e,z)) err(heller1);
--- 768,774 ----
zell(GEN e, GEN z, long prec)
{
long av=avma,ty,sw,fl;
!   GEN t,u,p1,p2,a,b,x1,u2,D = (GEN)e;

checkbell(e);
if (!oncurve(e,z)) err(heller1);
***************
*** 799,848 ****
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) >= gexpo((GEN)u)
!                         || gexpo((GEN)x1) >= gexpo((GEN)u));
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();
--- 789,839 ----
return gerepileupto(av,t);
}

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

!     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)));
!     r1 = gsub(x1,x0);
!     if (gcmp0(r1) || gexpo(r1) < gexpo(x1) - bit_accuracy(prec) + 5)
{
if (fl) break;
fl = 1;
}
else fl = 0;
}
!   u = gdiv(x1,a); t = gaddsg(1,u);
!   if (gcmp0(t) || 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 = gmul(u, glog(t,prec));

/* which square root? test the reciprocal function (pointell) */
if (!gcmp0(t))
{
!     GEN z1,z2;
!     int bad;

!     z1 = pointell(e,t,3); /* we don't need much precision */
!     /* Either z = z1 (ok: keep t), or z = z2 (bad: t <-- -t) */
!     z2 = invell(e, z1);
!     bad = (gexpo(gsub(z,z1)) > gexpo(gsub(z,z2)));
if (bad) t = gneg(t);
if (DEBUGLEVEL)
{
if (DEBUGLEVEL>4)
{
fprintferr("  z  = %Z\n",z);
!         fprintferr("  z1 = %Z\n",z1);
!         fprintferr("  z2 = %Z\n",z2);
}
fprintferr("ellpointtoz: %s square root\n", bad? "bad": "good");
flusherr();
__
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/
```