Karim BELABAS on Thu, 30 Sep 1999 14:59:18 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: ellinit() problem |
> 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/