[Haskell-cafe] Data.Complex.magnitude slow?

[Ian Lynagh]

On Thu, Jul 17, 2008 at 05:18:01PM +0200, Henning Thielemann wrote:
> 
> On Thu, 17 Jul 2008, stefan kersten wrote:
> 
> >i've attached an example program which seems to indicate that the 
> >magnitude function from Data.Complex is very slow compared to a more naive 
> >implementation (for Complex Float). on my machine (intel core2 duo, osx 
> >10.4) the CPU time using the library function is about 6-7 times as much 
> >as when using the other function. any ideas what might be going on? any 
> >flaws in my measurement code?
> 
> 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. Stefan, perhaps you can try timing
with the above code, and also with:

{-# SPECIALISE magnitude :: Complex Double -> Double #-}
magnitude :: (RealFloat a) => Complex a -> a
magnitude (x:+y) =  scaleFloat k
                     (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
                    where k  = max (exponent x) (exponent y)
                          mk = - k
                          sqr x = x * x

and let us know what the results are?


Thanks
Ian