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<[email protected]> 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 > [email protected] > http://www.haskell.org/mailman/listinfo/haskell-cafe > > _______________________________________________ Haskell-Cafe mailing list [email protected] http://www.haskell.org/mailman/listinfo/haskell-cafe
