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
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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe