Karim BELABAS on Tue, 9 May 2000 21:04:13 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: precision peculiarity |
[Ilya:] > > The bnfinit/quadclassunit computation don't guarantee anything about the > > actual precision of the output. They use an internal (much higher) precision > > using heuristics with the goal to get an exact class group, not an exact > > regulator. > > Precision policies of functions should be documented (with function > docs, but preferably also in a separate place). Is the answer's > precision based on argument precision? Is it based on default > precision? Are the internal calculations done with default precision, > and the answer has whatever precisioon it would have (:-()? Or is it > something from blue sky, as you describe above? /* Warning: long message ahead. I typed it as I was thinking/experimenting */ In this specific example, it's something from blue sky. The input has infinite precision (integer), and so has part of the output (class number and integral ideals generating the class group). The only guaranteed thing we know about the regulator at this point [assuming GRH and technical parameters were set correctly. This part is documented] is that it has the correct order of magnitude (not off by a multiplicative integer factor). Internal computations are done using an internal precision which depends on default precision, discriminant (proportional to log^2(D)), and whether a number of computations are successful or not (if not, double current precision and start over). It is extremely difficult to actually control what goes on there (and bnfinit is much much worse), but it is expected that in the final regulator at most 1 or 2 words will be off. If default precision is very small, this is a problem. On the other hand, it is trivial to check results a posteriori and give a regulator with guaranteed precision. For the example at hand: N = 139619637069004 R = 22944792.765428 /* crude approximation */ a = qfbred(Qfb(1,0,-N/4)) c = b = qfbred(a,1) (19:56) gp > n = round(R / component(c, 4)) %1 = 2753519 (19:56) gp > c = b ^ n %2 = Qfb(-4212135, 7476308, 4969241, 25010212.12717241812458238146) /* oops we overshoot */ (19:56) gp > n = round(R / component(c, 4)); /* ..., a few iterations later, we reach n = 2576292 */ (19:56) gp > c = b^n %11 = Qfb(1, 11816074, -8073882, 22944792.76542885320856372809) ^^^^^^^^^^^^^^^^ good approximation to R ! (19:56) gp > a == c %12 = 1 /* as expected we made a full turn on the whole principal cycle */ Now that we know n, we can compute R to any precision we want: (19:56) gp > \p2000 (19:56) gp > a = qfbred(Qfb(1,0,-N/4)); b = qfbred(a,1) (19:56) gp > component(b^n, 4) and we get R to 2000 digits in 2mn computing time, which should be the time needed to compute < 2 log_2(n) logarithms to this accuracy. A few benchmarks show that under GP it is at least 10 times slower than that, and 3 times SLOWER than to recompute the quadclassunit FROM SCRATCH starting from that precision !!! I wasn't aware that t_QFR was such an inefficient type [ should have known better: none of the PARI routines actually use this type, it's purely there for user interaction with GP ] After a glance at the code, it's obvious why: after each single reduction STEP a logarithm is computed, instead of accumulating a product and taking a single log at the end... ********************************************* Ok, this mail is long enough. Summary: (1) it's trivial to output regulator to any given accuracy, provided a crude initial value. It can be done very efficiently in library mode (without using the dedicated types!). (2) qfbred(t_QFR) [hence all t_QFR computations] are abysmally slow and can be just as easily improved. I won't do (1) before the stable release. I may do (2) if I have a little time (qfbred is unused: it won't affect stability...) Cheers, Karim. __ 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/