Ian Lynagh igloo at earth.li
Thu Jul 17 11:42:41 EDT 2008

```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

```