[Haskell-cafe] Complex numbers and (**)

[Lennart Augustsson]

A2: Yes, this seem unfortunate, so perhaps a different definition for
Complex is warranted.
Or maybe the default implementation for (**) should be changed so that
0**x is 0, except if x is 0 (in which case I think it should be
undefined).

  -- Lennart


On Sat, Aug 8, 2009 at 2:55 PM, Paul Sargent<psarge+haskell at gmail.com> wrote:
> Hi,
>
> First post to the cafe, so "Hello everybody!".
> Hope this is reasonable subject matter and not too long.
>
> I've been working on some algorithms that involved taking the n-th root of
> complex numbers. In my code I've implemented this as raising the complex
> number ('z') to 1/n using the (**) operator. Obviously, there are n roots,
> but I only need one of them so this is fine.
>
> Q1) Have I missed a method that's a little less general than 'raising to a
> Complex'? We have integer powers, but not integer roots?
>
> All seems to work fine, except I have a little wrapper function to prefer
> real roots of real numbers, until I started seeing NaNs appearing. This
> happened when I tried to take the root of 0+0i. In fact raising 0+0i to any
> power with (**) causes NaNs to appear. (^) and (^^) have no problem,
> assuming the calculation is one that can be represented with those
> operators. Neither is there a problem when the values being raised are not
> in complex form.
>
> Prelude Data.Complex> let xs = [0.0 :+ 0.0, 1.0 :+ 0.0, 2.0 :+ 0.0, 3.0 :+
> 0.0]
>
> Prelude Data.Complex> [x ^ 2 | x <- xs]
> [0.0 :+ 0.0,1.0 :+ 0.0,4.0 :+ 0.0,9.0 :+ 0.0]
>
> Prelude Data.Complex> [x ^^ 2 | x <- xs]
> [0.0 :+ 0.0,1.0 :+ 0.0,4.0 :+ 0.0,9.0 :+ 0.0]
>
> Prelude Data.Complex> [x ** 2 | x <- xs]
> [NaN :+ NaN,1.0 :+ 0.0,4.0 :+ 0.0,9.000000000000002 :+ 0.0]
>
> Prelude Data.Complex> let xs = [0.0,1.0,2.0,3.0]
> Prelude Data.Complex> [x ** 2 | x <- xs]
> [0.0,1.0,4.0,9.0]
>
> Digging deeper I've discovered this is because Complex inherits it's
> definition of (**) as "x ** y = exp (log x * y)". Well... the log of 0+0i is
> -Inf+0i. Multiply this by a real number in complex form and you end up with
> -Infinity * 0.0 as one of the terms. According to the IEEE floating point
> spec, this is NaN. That NaN propagates through exp, and you end up with NaN
> :+ NaN  as the result.
>
> Q2) Do people agree this is a bug in the definition of Data.Complex?
>
> Seems like the thing to do to fix this is have an instance of (**) for
> Data.Complex that special cases (0 :+ 0) ** _ to always return (0 :+ 0). An
> alternative would be to use the underlying non-complex (**) operator for
> arguments with no imaginary parts. One downside is that this would change
> the output of Complex (**) so that raising a real argument to a real power
> always produced a real result (which is actually what I want, but may not be
> what others expect / have got used to)
>
> Q3) Do people agree with these options? Any opinions? How would I submit a
> patch?
>
> I did send a mail to the glasgow-haskell-bugs list, but it doesn't appear to
> shown up in the archives, so I assume it didn't make it. It also didn't seem
> quite the right place as this is in the libraries. Apologies if anybody
> reading this is getting deja-vu.
>
> Paul
>
>
>
>
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>
>