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

[Henning Thielemann]

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.