Re: erfc() behavior change

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

* Walter Neumann [2005-05-13 16:19]:
> Thanks for the clarification. I think you are saying that for a zero real 
> only the precision is relevant.

Only the exponent is, yes.

> But the following is then still confusing:
> 
> ? 0.e-10
> %1 = 0.E-38
> ? precision(%)
> %2 = 48
> ? %1*10^40
> %3 = 0.E1
> ? precision(%)
> %4 = 0
> 
> I would expect the answers to be 0.E-48, 48, 0.E-10, 10 respectively.

%1: this is a bug. The answer should be 0.E-10 as input ( and just as 0E-10
  would produce ). Fixed in CVS.

%2 & %3: the precision is actually a number of significant _words_ (64 bits
  in your case). Then it is translated to a number of significant    
  _decimal_ digits, which is the documented user interface.

  This means that the precision is (roughly again !) a multiple of 
  19 = floor( log_10 2^64 ) on 64-bit machines, and a multiple of
   9 = floor( log_10 2^32 ) on 32-bit machines.

Here are the full rules for a real 0: what is internally stored as "any
real number x such that |x| < 2^-n"

1) is approximately printed as 0.E -N, with N = floor(-n * log(2)/log(10))
[ approximately since log(2)/log(10) is itself approximated to at most 16
significant digits... ]

2) has precision
  -- 0 if n <= 0
  -- floor( floor(n/2^sizeof(long)) * log(2)/log(10) ) otherwise (!).

That's why I said my initial explanation attempt was only approximately true...

One way to justify the above formula is that the only thing that makes
sense (given the internal representation) is the number of words.
precision() is just a wrapper, restricted to the GP user interface, never
used anywhere in the library.

It would be immensely simpler and clearer if PARI had a notion of
"bit accuracy" instead of the present (machine-dependent) "word accuracy". 
But this is a deep design bug, and I can't change that in the near future.

----------------------

To sum up: as a GP user, I do the following

1) Use #x (= length(x)) to dermine how many significant words there
truly are in a t_REAL. This is the only sensible data, everything else
is approximate.

2) A related side note: sizedigit() is likewise only approximate, use 
#Str(x) to get the exact number of decimal digits in an integer x
[ granted, sizedigit is slightly faster... ]

3) To estimate sizes, use only binary exponents
    ? install(gexpo, lG, expo)
    ? expo(0e-20)
    %1 = -67
[ that internal 'binary exponent' function should be exported in GP ]
     
Cheers,

    Karim.
-- 
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dep. de Mathematiques, Bat. 425   Fax: (+33) (0)1 69 15 60 19
Universite Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]