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
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@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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
On Sat, 8 Aug 2009, Paul Sargent wrote:
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.
I have written something on that topic: http://haskell.org/haskellwiki/Power_function I think the problem cannot be fully solved, especially not within the Haskell 98 numeric type classes. There is no satisfying implementation of (**) even for Float and Double.
On Sat, Aug 8, 2009 at 18:17, Henning Thielemann < lemming@henning-thielemann.de> wrote:
On Sat, 8 Aug 2009, Paul Sargent wrote:
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.
I have written something on that topic: http://haskell.org/haskellwiki/Power_function
I think the problem cannot be fully solved, especially not within the Haskell 98 numeric type classes. There is no satisfying implementation of (**) even for Float and Double.
Well, I'm an engineer rather than a mathematician, so I'm not 100% sure I followed everything on that page, but you seem to be saying that nth-root can't be done because we're using floating point, and floating point isn't accurate enough to represent 1/n. I basically agree with this, which was why I asked about nth-root functions, but for my uses (**) is good enough as long as complex zero works. (To be honest I'd expect any complex nth-root function to work with the polar form, there-by getting around some of the issues I think you're referring to.)
On Sun, 9 Aug 2009, Paul Sargent wrote:
Well, I'm an engineer rather than a mathematician, so I'm not 100% sure I followed everything on that page, but you seem to be saying that nth-root can't be done because we're using floating point, and floating point isn't accurate enough to represent 1/n.
Yes that's one problem.
I basically agree with this, which was why I asked about nth-root functions, but for my uses (**) is good enough as long as complex zero works.
In NumericPrelude we have the 'root' function for that purpose.
participants (3)
-
Henning Thielemann -
Lennart Augustsson -
Paul Sargent