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

[Aaron Denney]

On 2008-03-17, John Meacham [EMAIL PROTECTED] wrote:
 On Mon, Mar 17, 2008 at 12:59:09PM -0400, David Roundy wrote:
 foreign import ccall unsafe math.h log10 log10 :: Double - Double
 
 since in ghc CDouble and Double are identical.
 
 It's a bit sloppier, but shouldn't cause any trouble.  And I've no
 idea how realToFrac is implemented, but would worry about it behaving
 oddly... for instance when given NaNs.

 Yes. 'realToFrac' is inherently pretty broken and should be avoided
 whenever possible. It is not all all obvious to me what the correct
 primitive should be.. but we really should do something better for
 haskell'. relying on RULES firing as ghc currently does doesn't seem
 ideal..

 hmm.. maybe a 'FloatMax' type and have 'fromFloatMax' and 'toFloatMax'
 in 'Fractional' and 'Real' respectively? hmm.. hc has 'fromDouble' and
 'toDouble' there, but jhc also supports a 'Float128' type (when the
 underlying hardware does). so this still isn't quite right.

Well, the whole numeric hierarchy needs to be redone to support proper
mathematical structures like groups, rings, and fields.  Once that's
done, this might end up being clarified a bit.

-- 
Aaron Denney
--

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

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

[Wilhelm B. Kloke]

Jacob Schwartz [EMAIL PROTECTED] 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(n0)  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

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe