On Thu, 17 Jul 2008, Ian Lynagh wrote:
On Thu, Jul 17, 2008 at 05:18:01PM +0200, Henning Thielemann wrote:
Complex.magnitude must prevent overflows, that is, if you just square
1e200::Double you get an overflow, although the end result may be also
around 1e200. I guess, that to this end Complex.magnitude will separate
mantissa and exponent, but this is done via Integers, I'm afraid.
Here's the code:
{-# SPECIALISE magnitude :: Complex Double -> Double #-}
magnitude :: (RealFloat a) => Complex a -> a
magnitude (x:+y) = scaleFloat k
(sqrt ((scaleFloat mk x)^(2::Int) + (scaleFloat mk
y)^(2::Int)))
where k = max (exponent x) (exponent y)
mk = - k
So the slowdown may be due to the scaling, presumably to prevent
overflow as you say. However, the e^(2 :: Int) may also be causing a
slowdown, as (^) is lazy in its first argument; I'm not sure if there is
a rule that will rewrite that to e*e.
I guess the rule should be INLINE.
Indeed, here you can see
http://darcs.haskell.org/packages/base/GHC/Float.lhs
that scaleFloat calls decodeFloat and encodeFloat. Both of them use
Integer. I expect that most FPUs are able to divide a floating point
number into exponent and mantissa, but GHC seems not to have a primitive
for it? As a quick work-around, Complex.magnitude could check whether the
arguments are too big, then use scaleFloat and otherwise it should use the
naive algorithm.
In the math library of C there is the function 'frexp'
http://www.codecogs.com/reference/c/math.h/frexp.php?alias=
but calling an external functions will be slower than a comparison that
can be performed by a single FPU instruction.
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
