Re: ellap() bug?

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

On Wed, 23 Jan 2002, Iwao KIMURA wrote:
> I heard that ellap() returns wrong answer for large p.
>
> ? e=ellinit([0,0,0,0,24])
> ? p=557018866141
> ? ellap(e,p)
> time = 2mn, 24,587 ms.
> %4 = 47745597502
> ? 47745597502 < 2*sqrt(p)
> time = 0 ms.
> %5 = 0
>
> ellap(e,p) returns 47745597502, but this is not smaller than
> 2*sqrt(p), violating Hasse's theorem.

Stack corruption, easy to fix once the problem is spotted. Patch follows
[should apply to any of the 2.1.* and 2.2.* releases; I updated both CVS
archives (stable and unstable)].

(20:39) gp > e=ellinit([0,0,0,0,24]); p=557018866141;
(20:41) gp > ellap(e,p)
time = 10 ms.
%2 = 1228567

    Karim.

Index: src/modules/elliptic.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/modules/elliptic.c,v
retrieving revision 1.35
retrieving revision 1.37
diff -c -r1.35 -r1.37
*** src/modules/elliptic.c	2002/01/08 22:19:34	1.35
--- src/modules/elliptic.c	2002/01/23 19:56:48	1.37
***************
*** 1339,1352 ****
  static GEN
  addsell(GEN e, GEN z1, GEN z2, GEN p)
  {
!   GEN p1,p2,x,x1,x2,y,y1,y2;
!   long av = avma;

    if (!z1) return z2;
    if (!z2) return z1;

    x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
    x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
    p2 = subii(x2, x1);
    if (p2 == gzero)
    {
--- 1339,1353 ----
  static GEN
  addsell(GEN e, GEN z1, GEN z2, GEN p)
  {
!   GEN z,p1,p2,x,x1,x2,y,y1,y2;
!   ulong av;

    if (!z1) return z2;
    if (!z2) return z1;

    x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
    x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
+   z = cgetg(3, t_VEC); av = avma;
    p2 = subii(x2, x1);
    if (p2 == gzero)
    {
***************
*** 1358,1368 ****
    else p1 = subii(y2,y1);
    p1 = mulii(p1, mpinvmod(p2, p));
    p1 = resii(p1, p);
!   x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p);
    y = negi(addii(y1, mulii(p1,subii(x,x1))));
!   avma = av; p1 = cgetg(3,t_VEC);
!   p1[1] = licopy(x);
!   p1[2] = lmodii(y, p); return p1;
  }

  /* z1 <-- z1 + z2 */
--- 1359,1370 ----
    else p1 = subii(y2,y1);
    p1 = mulii(p1, mpinvmod(p2, p));
    p1 = resii(p1, p);
!   x = subii(sqri(p1), addii(x1,x2));
    y = negi(addii(y1, mulii(p1,subii(x,x1))));
!   x = modii(x,p);
!   y = modii(y,p); avma = av;
!   z[1] = licopy(x);
!   z[2] = licopy(y); return z;
  }

  /* z1 <-- z1 + z2 */

-- 
Karim Belabas                    Tel: (+33) (0)1 69 15 57 48
Dép. de Mathematiques, Bat. 425  Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud             Email: Karim.Belabas@math.u-psud.fr
F-91405 Orsay (France)           http://www.math.u-psud.fr/~belabas
--
PARI/GP Home Page: http://www.parigp-home.de/