[Haskell-cafe] Re: floating point operations and representation
Wilhelm B. Kloke
wb at arb-phys.uni-dortmund.de
Thu Mar 13 05:51:39 EDT 2008
Jacob Schwartz <quark at bluespec.com> schrieb:
> I have two questions about using the Double data type and the
> operations in the Floating typeclass on a computer that uses IEEE
> floating point numbers.
>
> I notice that the Floating class only provides "log" (presumably log
> base 'e') and "logBase" (which, in the latest source that I see for
> GHC is defined as "log y / log x"). However, in C, the "math.h"
> library provides specific "log2" and "log10" functions, for extra
> precision. A test on IEEE computers (x86 and x86-64), shows that for
> a range of 64-bit "double" values, the answers in C do differ (in the
> last bit) if you use "log2(x)" and "log10(x)" versus "log (x) /
> log(2)" and "log(x) / log(10)".
This is to be expected in about 25% of cases, as the GHC definition rounds
two times, whereas the implementation provided in the SUN math library
(file lib/msun/src/e_log10.c on FreeBSD) uses a representation in two
doubles for log(10) to avoid one of the rounding errors:
* Return the base 10 logarithm of x
*
* Method :
* Let log10_2hi = leading 40 bits of log10(2) and
* log10_2lo = log10(2) - log10_2hi,
* ivln10 = 1/log(10) rounded.
* Then
* n = ilogb(x),
* if(n<0) n = n+1;
* x = scalbn(x,-n);
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
--
Dipl.-Math. Wilhelm Bernhard Kloke
Institut fuer Arbeitsphysiologie an der Universitaet Dortmund
Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257
PGP: http://vestein.arb-phys.uni-dortmund.de/~wb/mypublic.key
More information about the Haskell-Cafe
mailing list