Karim Belabas on Fri, 13 May 2005 17:41:23 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: erfc() behavior change |
* 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]