Karim Belabas on Thu, 12 May 2005 11:02:21 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: erfc() behavior change


* Walter Neumann [2005-05-12 07:31]:
> On Wed, 11 May 2005, Igor Schein wrote:
>>\\ ver 2.2.9
>>? for(k=1,10,print(k" "precision(erfc(2^k))" "precision(erfc(-2^k))))
>>...
>>10 455407 455446
>>\\ ver 2.2.10
>>? for(k=1,10,print(k" "precision(erfc(2^k))" "precision(erfc(-2^k))))
>>...
>>10 38 455446
>
> The second (2.2.10) looks better to me:
>
> GP/PARI CALCULATOR Version 2.2.11 (development CHANGES-1.1205)
>
> ? erfc(2^10)
> %1 = 9.342620665669385261706140592 E-455395
> ? precision(%)
> %2 = 28
> ? erfc(-2^10)
> %3 = 2.000000000000000000000000000
> ? precision(%)
> %4 = 455427
> ? 2-%3
> %5 = 9.342620665669385261706140592 E-455395
> ? precision(%)
> %6 = 38

Indeed. As for the ridiculous accuracy of %3 above, we have conflicting 
"specifications":
1) PARI functions give as precise a result as is possible from the input,
2) floating point computations are meant to foster speed by truncating
operands.

Only 1) is specified in the documentation, 2) is only a general understanding. 
And a rather misleading one as far as PARI is concerned; it is a common
source of misapprehension to assume that

* 'realprecision' is  "the relative accuracy used to truncate operands in 2)".
Which it is not: it is used to convert exact objects to inexact ones.

* operands with n digits of accuracy will yield a result with at most the
same accuracy. Which is wrong: indeed 1 + 1e-50000 may be computed to
more than 50000 digits of accuracy.


In most cases, the second behaviour is a bug from the user's point of view
(what's the point of getting 455000 trailing zeroes ?). I believe it is
better to stick to strict specifications and let the user sort out
numerical problems from this point.

The 2.2.9 problem was quite different: the _apparent_ accuracy of erfc(huge)
was huge, but almost all printed decimals were wrong due to catastrophic
cancellation. I fixed it so that only meaningful digits are output
[ and so that complex arguments are accepted, but that's irrelevant to our
present topic ].

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]