Re: PARI:Failure of ellheight on elliptic curve [0,-1,0,-784,8704]

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

On Fri, 6 Jul 2001, Bill Allombert wrote:
> On Thu, Jul 05, 2001 at 06:14:52PM +0100, Martin Prickett wrote:
>> I am Martin Prickett, a phD student of John Cremona at Nottingham. I am
>> using PARI and have discovered that if I take curve
>> e=ellinit([0,-1,0,-784,8704]) and point P=[-8,120] then
>> ellheight(e,P) crashes the computer by causing it to run out of memory on
>> Linux machines. However on Dec alpha machines (which are 64 bit) it gives
>> the correct answer.
>
> First I suppose your Linux machine is a x86 ?
>
> I try this on a x86 Linux machine and on a sparc solaris machine
> and I got.
>   ***   the PARI stack overflows !
>   current stack size: 33554432 (32.000 Mbytes)
>
> which is not a computer crash, just a stack overflow condition.
>
> I also try it on a alpha Linux machine and it work
> ? ellheight(e,P)
> %3 = 0.74629920869451370890943760363155317690
>
> So you are right is a 32bit-only issue.

Numerical instability in ellpointtoz (pre-AGM phase): the AGM entered an
infinite loop since requested accuracy couldn't possibly be attained [was
still possible on 64bit machines]. Eventually, the GP stack would overflow.

Please apply the patch below [CVS repository is also up to date]

    Karim.

Index: src/modules/elliptic.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/modules/elliptic.c,v
retrieving revision 1.21
diff -c -r1.21 elliptic.c
*** src/modules/elliptic.c	2001/03/20 00:04:00	1.21
--- src/modules/elliptic.c	2001/07/06 09:17:39
***************
*** 766,778 ****
    b = gmul2n(w,-1);
    r1 = gsub(a,b);
    p1 = gadd(x, gmul2n(gadd(e1,r0),-1));
!   if (gcmp0(p1))
!     p1 = gsqrt(gneg_i(gmul(a,r1)),prec);
!   else
!   {
!     GEN delta = gdiv(gmul(a,r1),gsqr(p1));
!     p1 = gmul2n(gmul(p1,gaddsg(1,gsqrt(gsubsg(1,gmul2n(delta,2)),prec))),-1);
!   }
    *pta = a; *ptb = b;
    return gmul(p1,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(p1,r1),p1),prec)),-1)));
  }
--- 766,773 ----
    b = gmul2n(w,-1);
    r1 = gsub(a,b);
    p1 = gadd(x, gmul2n(gadd(e1,r0),-1));
!   p1 = gmul2n(p1,-1);
!   p1 = gadd(p1, gsqrt(gsub(gsqr(p1), gmul(a,r1)), prec));
    *pta = a; *ptb = b;
    return gmul(p1,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(p1,r1),p1),prec)),-1)));
  }

-- 
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://www.parigp-home.de/