Re: ellinit() problem

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

> I noticed the problem about the following command:
> 
> ellinit([0,0,0,-21002000290904852499,37045810753211102630776943500])
> 
> It keeps doubling the memory.  I ran it to as high as 1GB RAM.
> Debugging showed that it stays forever inside for(;;) loop in do_agm()
> function.  I couldn't understand, however, what exactly eats the memory.

Precision loss due to the fact that two of the roots are really close. An
error term estimate happens to be an inexact 0 with a (relatively) high
exponent [the minimal one that its precision warrants but still too big...],
hence it's never recognized as being small enough.


The best solution would be to compute the roots to a higher accuracy on
startup, and iterate as long as the periods, etc. are not obtained to the
desired accuracy [using the fact that precision of the objects decreases when
it can't be guaranteed anymore]. That wouldn't be completely rigorous since
most PARI functions don't guarantee the precision (polroots is an exception,
not the rule), and that in any case it would not deal with the accumulation
of errors (one would need some form of interval arithmetic for that...). But
it would be much more satisfactory.

For the time being, I'm going with the easy solution:

  This patch considers any real 0 as satisfactory, and lets it go through. If
  that happens, functions using inexact components of the elliptic curve
  structure will return objects with lower precision than the current
  realprecision.

If the precision drops below than realprecision/2, one will get "precision
too low in initell" error messages. Hence meaningless results won't be
returned; only less precise ones.

Karim.

Index: src/modules/elliptic.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/modules/elliptic.c,v
retrieving revision 1.2
diff -c -r1.2 elliptic.c
*** src/modules/elliptic.c      1999/09/23 17:50:57     1.2
--- src/modules/elliptic.c      1999/09/30 11:38:16
***************
*** 183,189 ****
      r1=gsub(a1,b1);
      p1=gsqrt(gdiv(gadd(x,r1),x),prec);
      x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
!     if (gexpo(r1) <= G + gexpo(b1)) break;
    }
    if (gprecision(x1)*2 <= (prec+2))
      err(talker,"precision too low in initell");
--- 183,189 ----
      r1=gsub(a1,b1);
      p1=gsqrt(gdiv(gadd(x,r1),x),prec);
      x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
!     if (gcmp0(r1) || gexpo(r1) <= G + gexpo(b1)) break;
    }
    if (gprecision(x1)*2 <= (prec+2))
      err(talker,"precision too low in initell");

__
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/