Re: ellpointtoz() bug?

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

> 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);
!     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();
--- 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/