[Haskell-cafe] Re: floating point operations and representation

[Wilhelm B. Kloke]

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