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

[Paul Sargent]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.haskell.org/pipermail/haskell-cafe/attachments/20090808/d65950ba/attachment.html