Re: precision peculiarity

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

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