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.