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]