ellisoncurve bug & suggested patch

[Gerhard Niklasch <nikl@mathematik.tu-muenchen.de>]

Xavier Roblot recently pointed out the following weirdness to me:
> BTW, someone working on elliptic curves sent me the following
> bug report (which is still present in the new version):
> 
> ? E40 = [0,0,0,-7,-6];
> ? pol = x^3-7*x-7;
> // clearly, [a,1] is on the curve for any root a of pol;
> // but PARI doesn't agree
> ? r = polroots(pol);
> ? ellisoncurve(E40,[r[1],1])
> %1 = 1
> ? ellisoncurve(E40,[r[2],1])
> %2 = 0
> ? ellisoncurve(E40,[r[3],1])
> %3 = 0

There were in fact two bugs.  One was of the one-character-typo
kind and had been introduced waaaaaaaaay back in 2.0.2.alpha.
The other was the logic in the internal function oncurve().
It computes the lefthand side  y^2 + a_1*x*y + a_3*y  and the
righthand side  (cubic in x)  of the equation, and then their
difference.  When the difference is zero, it returns true (1).
When the difference is nonzero and both sides are exact, it
returns false (0).  So far, so good.  When one of the sides is
exact and the other isn't, as was the case above  (a_1 is zero
and x is the only inexact quantity in sight, so the lefthand
side is exact),  it used to return false, which accounts for
the strange output above.

With the patch below, results look more intuitive:
(01:28) gp > ellisoncurve(E40,[r[1],1])
%4 = 1
(01:28) gp > ellisoncurve(E40,[r[2],1])
%5 = 1
(01:28) gp > ellisoncurve(E40,[r[3],1])
%6 = 1
(01:29) gp > ellisoncurve(E40,[r[1],1.000000000000000000001])
%9 = 0
(01:29) gp > ellisoncurve(E40,[r[1],1.0000000000000000000000001])
%10 = 0
(01:29) gp > ellisoncurve(E40,[r[1],1.00000000000000000000000000001])
%11 = 1
(01:29) gp > ellisoncurve(E40,[r[3]+0.000000000000000001,1])
%12 = 0
(01:29) gp > ellisoncurve(E40,[r[3]+0.00000000000000000000000001,1])
%13 = 0
(01:29) gp > ellisoncurve(E40,[r[3]+0.00000000000000000000000000000001,1])
%14 = 1

Diffs are against the latest pre-2.0.17 snapshot, but it should go in
with mild protest if you have 2.0.16.beta  (line numbers are off by 2
in that case).  Apply by hand against older versions;  this particular
corner of code hasn't changed during the last several beta versions.

Please do let me know if you find this breaking anything else.

Enjoy, Gerhard


$ diff -c src/modules/elliptic.c.orig src/modules/elliptic.c
*** src/modules/elliptic.c.orig	Wed Sep  1 01:03:15 1999
--- src/modules/elliptic.c	Wed Sep  1 01:24:14 1999
***************
*** 449,455 ****
  static long
  ellexpo(GEN e)
  {
!   long i,k2, k = -gexpo((GEN)e[1]);
    for (i=2; i<6; i++)
    {
      k2 = gexpo((GEN)e[i]);
--- 449,455 ----
  static long
  ellexpo(GEN e)
  {
!   long i,k2, k = gexpo((GEN)e[1]);
    for (i=2; i<6; i++)
    {
      k2 = gexpo((GEN)e[i]);
***************
*** 458,463 ****
--- 458,468 ----
    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.-- Changed to get rid of the resulting
+    weirdness;  hope this doesn't break anything p-adically. --GN1999Sep01 */
  int
  oncurve(GEN e, GEN z)
  {
***************
*** 469,476 ****
    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;
!   if (!q) { avma=av; return 0; } /* one of p1, p2 is exact */
    /* the constant 0.93 is arbitrary */
    q = (gexpo(x)-ellexpo(e) < - bit_accuracy(q) * 0.93);
    avma = av; return q;
--- 474,482 ----
    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;