[Haskell-cafe] Complex numbers and (**)
Paul Sargent
psarge+haskell at gmail.com
Sat Aug 8 08:55:49 EDT 2009
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
More information about the Haskell-Cafe
mailing list