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

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

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

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);
!     r1=gsub(a,b);
!     if (gexpo(gsub(x1,x0)) < gexpo(x1) - bit_accuracy(prec) + 5)
{
if (fl) break;
fl = 1;
}
else fl = 0;
}
!   if (gexpo(t) <  5 - bit_accuracy(prec))
t = negi(gun);
else
!   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;

!     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 (DEBUGLEVEL)
{
if (DEBUGLEVEL>4)
{
fprintferr("  z  = %Z\n",z);
!         fprintferr("  u  = %Z\n",u);
!         fprintferr("  x1 = %Z\n",x1);
}
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);
!     r1 = gsub(a,b);
!     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
!   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;

!     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 (DEBUGLEVEL)
{
if (DEBUGLEVEL>4)
{
fprintferr("  z  = %Z\n",z);
!         fprintferr("  z1 = %Z\n",z1);
!         fprintferr("  z2 = %Z\n",z2);
}