#8539 closed bug (fixed)
Data.Complex shouldn't use default implementation of (**)
Reported by: | jjaredsimpson | Owned by: | |
---|---|---|---|
Priority: | low | Milestone: | 8.0.1 |
Component: | Prelude | Version: | 7.6.3 |
Keywords: | Cc: | ekmett | |
Operating System: | Unknown/Multiple | Architecture: | Unknown/Multiple |
Type of failure: | Incorrect result at runtime | Test Case: | |
Blocked By: | Blocking: | ||
Related Tickets: | Differential Rev(s): | ||
Wiki Page: |
Description
Prelude Data.Complex> 0**2 0.0 Prelude Data.Complex> 0**2 :: Complex Double NaN :+ NaN Prelude Data.Complex> exp $ 2 * log 0 :: Complex Double NaN :+ NaN Prelude Data.Complex> log 0 :: Complex Double (-Infinity) :+ 0.0 Prelude Data.Complex> 2 * it :: Complex Double (-Infinity) :+ NaN Prelude Data.Complex> exp it :: Complex Double NaN :+ NaN
So Complex uses the default implementation of (**). Then when 2*(-inf :+ 0) is evaluated. We do (2 * -inf - 0*0) :+ (2*0 + -inf*0). Which because of -inf*0 sets the imaginary part to NaN.
Then exp (-inf :+ NaN) = exp x cos y :+ exp x sin y which becomes 0 * cos NaN :+ 0 * sin NaN. So we end up with NaN :+ NaN.
Ross Paterson suggested:
I would say the default implementation of (**) is wrong: to match the Float/Double instances it should be
x ** y = if x == 0 then 0 else exp (log x * y)
Attachments (5)
Change History (51)
comment:1 Changed 6 years ago by
comment:2 Changed 6 years ago by
Agree with isaacdupree that
0**
0 == 1.0
The above
x ** y = if x == 0 then 0 else exp (log x * y)
is incorrect for negative y. Consider
$ ghci GHCi, version 7.6.3: http://www.haskell.org/ghc/ :? for help Prelude> 0.0**(-2.0) Infinity
There are other interesting cases. If the real part of y is positive, then the result should be 0. If the real part is negative, then the result should be Infinity. If y is imaginary, the result should be NaN.
For the most part,
Infinity ** y = 0 ** (negate y)
I won't address -0 and 0:+(-0) because they deserve more controversy than I can handle.
comment:3 Changed 5 years ago by
I think this will be enough:
(x:+0)**(y:+0) = (x ** y) :+ 0 x ** y = exp (log x * y)
Testing:
Prelude Data.Complex> 0**0 :: Complex Double 1.0 :+ 0.0 Prelude Data.Complex> 0**2 :: Complex Double 0.0 :+ 0.0 Prelude Data.Complex> 0**(-2) :: Complex Double Infinity :+ 0.0
Changed 5 years ago by
Attachment: | complex_power.path added |
---|
comment:4 Changed 5 years ago by
Status: | new → patch |
---|
comment:5 Changed 5 years ago by
The patch fails for negative x.
Prelude Data.Complex> (-1)**(0.5) :: Complex Double NaN :+ 0.0
The result needs to be near either 0:+1 or 0:+(-1).
The existing exp (log x * y) is basically good as long as x is not 0 or infinity.
comment:6 Changed 5 years ago by
What about just adding check x>=0?
(x:+0) ** (y:+0) | x >= 0 = (x ** y) :+ 0 x ** y = exp (log x * y)
It will be more accurate when raising real numbers to real power:
Prelude Data.Complex> exp (log (2) * 4) :: Complex Double 15.999999999999998 :+ 0.0 Prelude Data.Complex> 2**4 :: Double 16.0
comment:7 follow-up: 8 Changed 5 years ago by
Adding the x>=0 is an improvement.
If it's going to handle one infinity (i.e. +Infinity:+0) it's best to handle them all. There are lots of infinite complex numbers.
(x:+0) ** (y:+0) | x >= 0 = (x ** y) :+ 0 (re:+im) ** 0 | isInfinite re || isInfinite im = 1:+0 x ** y = exp (log x * y)
comment:8 Changed 5 years ago by
I think x raised to 0 must always be equal 1. Even for NaN:
Prelude> (1/0*0) :: Double NaN Prelude> (1/0*0) ** 0 :: Double 1.0
New complex power implementation:
_ ** (0:+0) = 1 :+ 0 (x:+0) ** (y:+0) | x >= 0 = (x ** y) :+ 0 x ** y = exp (log x * y)
But what about raising zero to complex power? Current implementation gets NaN:
Prelude Data.Complex> exp $ log 0 * (1:+1) :: Complex Double NaN :+ NaN
Python throws an exception:
>>> 0**(1+1j) Traceback (most recent call last): File "<stdin>", line 1, in <module> ZeroDivisionError: 0.0 to a negative or complex power
But wolphramalpha gets 0, infinity or indeterminate.
comment:9 Changed 5 years ago by
We can also add special case for raising to 1 to increase accuracy:
Prelude Data.Complex> exp $ log (1:+2) * (1:+0) :: Complex Double 0.9999999999999998 :+ 2.0000000000000004
By adding this line:
x ** (1:+0) = x
But I'm not sure. Is it necessary? This function already has 3 definitions.
comment:11 Changed 5 years ago by
And many other functions can't handle infinity:
Prelude Data.Complex> sqrt (1/0) :: Double Infinity Prelude Data.Complex> sqrt (1/0) :: Complex Double NaN :+ NaN Prelude Data.Complex> log (1/0) :: Double Infinity Prelude Data.Complex> log (1/0) :: Complex Double NaN :+ NaN Prelude Data.Complex> exp (1/0) :: Double Infinity Prelude Data.Complex> exp (1/0) :: Complex Double NaN :+ NaN Prelude Data.Complex> (1/0) * (1/0) :: Double Infinity Prelude Data.Complex> (1/0) * (1/0) :: Complex Double NaN :+ NaN Prelude Data.Complex> (1/0) / 2 :: Double Infinity Prelude Data.Complex> (1/0) / 2 :: Complex Double NaN :+ NaN
Is it ok? Do we need tickets for each of them? In my opinion complex numbers with zero imaginary part should behave like simple real numbers in all such functions.
Changed 5 years ago by
Attachment: | complex_functions.patch added |
---|
comment:12 Changed 5 years ago by
Testing complex_functions.patch:
Double | Complex Double (patched) | Complex Double (current) | |
---|---|---|---|
1/0 | Infinity | Infinity :+ 0.0 | Infinity :+ 0.0 |
1/0*0 | NaN | NaN :+ 0.0 | NaN :+ NaN |
1/0/0 | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
1/0/2 | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
(1/0)*(1/0) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
sqrt (1/0) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
sqrt (-1/0) | NaN | 0.0 :+ Infinity | NaN :+ NaN |
sqrt (-4) | NaN | 0.0 :+ 2.0 | (-0.0) :+ 2.0 |
exp (1/0) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
exp (-1/0) | 0.0 | 0.0 :+ 0.0 | NaN :+ NaN |
log (1/0) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
log (-1/0) | NaN | Infinity :+ 3.141592653589793 | NaN :+ NaN |
0**0 | 1.0 | 1.0 :+ 0.0 | NaN :+ NaN |
0**2 | 0.0 | 0.0 :+ 0.0 | NaN :+ NaN |
0**(-2) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
2**4 | 16.0 | 16.0 :+ 0.0 | 15.999999999999998 :+ 0.0 |
3**1 | 3.0 | 3.0 :+ 0 | 3.0000000000000004 :+ 0.0 |
This even partially fixes old raise to power implementation:
Double | Complex Double (patched) | Complex Double (current) | |
---|---|---|---|
exp (log 0 * 2) | 0.0 | 0.0 :+ 0.0 | NaN :+ NaN |
exp (log 0 * 0) | 1.0 | NaN :+ 0.0 | NaN :+ NaN |
exp (log 0 * (-2)) | Infinity | Infinity :+ 0.0 | NaN :+ NaN |
What do you think?
comment:13 Changed 5 years ago by
In my opinion complex numbers with zero imaginary part should behave like simple real numbers in all such functions.
This is a naive understanding of the complex numbers. They are not a plain extension of the real numbers; they are missing the linear ordering. So you can't say that log(0)= -infinity as you do with the real numbers, because there are additional ways in which to approach 0.
comment:14 Changed 5 years ago by
But wolphramalpha gets 0, infinity or indeterminate.
which agrees with what I said earlier,
If the real part of y is positive, then the result should be 0. If the real part is negative, then the result should be Infinity. If y is imaginary, the result should be NaN.
comment:15 Changed 5 years ago by
Wolfram Alpha is a good comparison. They're pros. I'm just a math major who remembers his complex analysis pretty well.
comment:16 follow-up: 17 Changed 5 years ago by
We can also add special case for raising to 1 to increase accuracy But I'm not so sure.
Special cases are not just messy code; they produce messy values. I've seen an explanation from a numerical analyst that they need a function that's mathematically monotonic, such as log, to be numerically monotonic when computed. It's hard to guarantee smooth computational behavior if special cases are thrown in, in order to gain perfect accuracy.
In my opinion it would be better to leave
x ** y = exp (log x * y)
in effect, except when x is 0 or infinite.
comment:17 Changed 5 years ago by
Replying to Scott Turner:
We can also add special case for raising to 1 to increase accuracy But I'm not so sure.
Special cases are not just messy code; they produce messy values. I've seen an explanation from a numerical analyst that they need a function that's mathematically monotonic, such as log, to be numerically monotonic when computed. It's hard to guarantee smooth computational behavior if special cases are thrown in, in order to gain perfect accuracy.
In my opinion it would be better to leave
x ** y = exp (log x * y)in effect, except when x is 0 or infinite.
Ok. Thank you for explanation. But why we want to handle infinity**0 and don't want to handle infinity**y or x**infinity? And what about handling infinity in other functions such as sqrt, exp, (/) or (*)? Why is (**) so special?
comment:18 follow-ups: 19 20 Changed 5 years ago by
I don't know that (**) is more special. Maybe we are talking about it because it has the rule
x ** 0 == 1
which is peculiarly effective. It does battle with another rule
0 ** x == 0
and in the judgement of most numerical analysis, it wins the fight.
I think the reason people aren't diving in to address how all the other operators handle complex infinities is that it's very complex. Your patch handles a few cases well, but it doesn't come near to handling them all. For example, it doesn't take care of
(infinity:+infinity)*(infinity:+infinity) == (NaN:+infinity)
which is poor because the result should be a plain infinite value. The same kind of thing can happen, even if one of the operands is finite. Then there are cases in which there's an intermediate overflow but the best result is not infinite. That's because
(x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x')
subtracts after multiplying. The intermediate x*x' can be unrepresentable, while x*x'-y*y' is representable.
It's all a lot of fun. It would be nice work if someone would pay us for doing it.
comment:19 Changed 5 years ago by
Replying to Scott Turner:
It's all a lot of fun. It would be nice work if someone would pay us for doing it.
I can not disagree with you. It's very interesting even for free, at least for now because it's my first try to contribute to open source project.
comment:20 Changed 5 years ago by
Replying to Scott Turner:
Maybe we are talking about it because it has the rule
x ** 0 == 1
I think it is completely safe to define:
_ ** (0:+0) = 1 :+ 0
Because exp $ any_finite_number * (0:+0)
always equals 1:+0
.
The only difference is for NaN, but NaN**0 :: Double
equals 1
too.
Changed 5 years ago by
Attachment: | complex_power_without_inf.patch added |
---|
comment:21 follow-ups: 22 23 Changed 5 years ago by
NaN**0 :: Double equals 1 too.
See http://wikipedia.org/wiki/NaN#Signaling_NaN for an exception to that. The same Wikipedia article has interesting things to say about 0**0 and infinity**0 also.
comment:22 Changed 5 years ago by
Replying to Scott Turner:
NaN**0 :: Double equals 1 too.
See http://wikipedia.org/wiki/NaN#Signaling_NaN for an exception to that. The same Wikipedia article has interesting things to say about 0**0 and infinity**0 also.
Interesting. It says:
The 2008 version of the IEEE 754 standard says that pow(1,qNaN) and pow(qNaN,0) should both return 1 since they return 1 whatever else is used instead of quiet NaN.
Should we add another definition?
(1:+0) ** _ = 1:+0
To fix this:
Prelude Data.Complex> 1 ** (0/0) :: Double 1.0 Prelude Data.Complex> 1 ** (0/0) :: Complex Double NaN :+ NaN
comment:23 Changed 5 years ago by
Other floating functions seems to be ok (or almost ok) with infinity:
Prelude Data.Complex> sqrt $ 1/0 :+ 0 Infinity :+ 0.0 Prelude Data.Complex> sqrt $ 1/0 :+ 1/0 Infinity :+ NaN Prelude Data.Complex> sqrt $ (-1/0) :+ 0 0.0 :+ Infinity Prelude Data.Complex> sqrt $ (-1/0) :+ (-1/0) NaN :+ (-Infinity) Prelude Data.Complex> log $ 1/0 :+ 0 Infinity :+ 0.0 Prelude Data.Complex> log $ 1/0 :+ 1/0 Infinity :+ NaN Prelude Data.Complex> log $ (-1/0) :+ 0 Infinity :+ 3.141592653589793 Prelude Data.Complex> log $ 0 :+ 0 (-Infinity) :+ 0.0 Prelude Data.Complex> log $ (-1/0) :+ (-1/0) Infinity :+ NaN Prelude Data.Complex> exp $ 1/0 :+ 1 Infinity :+ Infinity Prelude Data.Complex> exp $ 1/0 :+ 0 Infinity :+ NaN Prelude Data.Complex> exp $ (-1/0) :+ 0 0.0 :+ 0.0
But (**) can not raise to infinity:
Prelude Data.Complex> 2 ** (1/0 :+ 1) NaN :+ NaN Prelude Data.Complex> 2 ** (1/0 :+ 0) NaN :+ NaN Prelude Data.Complex> 2 ** ((-1/0) :+ 0) NaN :+ NaN
Is it need to be fixed too?
comment:24 Changed 5 years ago by
Status: | patch → infoneeded |
---|
I haven't followed this whole thread, but it seems like there are still a few questions to answer before merging anything; if anyone would like to summarize the discussion here and take ownership of this to get a full patch in, that would be great.
comment:25 Changed 5 years ago by
The early answer,
If the real part of y is positive, then the result should be 0. If the real part is negative, then the result should be Infinity. If y is imaginary, the result should be NaN.
Infinity ** y = 0 ** (negate y)
provides the basis for a satisfactory solution. One improvement is to take Yalas's recommendation regarding NaN**0. In code, the implementation is
x ** y = case (x,y) of (_ , (0:+0)) -> 1 :+ 0 ((0:+0), (re:+_)) | re > 0 -> 0 :+ 0 | re < 0 -> inf :+ 0 | otherwise -> nan :+ nan ((re:+im), y) | (isInfinite re || isInfinite im) -> case y of (exp_re:+_) | exp_re > 0 -> inf :+ 0 | exp_re < 0 -> 0 :+ 0 | otherwise -> nan :+ nan (x, y) -> exp (log x * y) where inf = 1/0 nan = 0/0
Regarding the questions raised in the earlier discussion: Q: What to do about the other operations than **? A: That should be a separate bug report.
Q: special case for raising to 1 to increase accuracy A: Addressed in the earlier discussion.
Q: Do we need some tests? A: See examples below.
Q: why we want to handle infinity**0 and don't want to handle infinity**y or x**infinity? A: We do handle infinity**y. For x**infinity, the returned value is NaN:+NaN. The above patch could be refined using a directional interpretation of complex infinities, i.e. Infinity:+0, Infinity:+Infinity, 0:+Infinity, -Infinity:+0, etc. would each yield different results, and the result would also depend on log x. To some extent that might be managed by revising the exp function. It gets fiddly so I recommend making the above improvement, and leaving x**infinity for another day.
Examples
existing | patched | Wolfram Alpha | ||
---|---|---|---|---|
0 ** 2 | NaN :+ NaN | 0.0 :+ 0.0 | 0 | |
0 ** 0 | NaN :+ NaN | 1.0 :+ 0.0 | indeterminate | |
0 ** -2 | NaN :+ NaN | Infinity :+ 0.0 | complex infinity | |
inf:+0 ** 2 | NaN :+ NaN | Infinity :+ 0.0 | infinity | |
inf:+0 ** 0 | NaN :+ NaN | 1.0 :+ 0.0 | indeterminate | |
inf:+0 ** -2 | NaN :+ NaN | 0.0 :+ 0.0 | 0 | |
-1 ** 0.5 | 6.12e-17 :+ (-1.0) | 6.12e-17 :+ (-1.0) | 0:+1 | Note 2 |
0 ** (1:+1) | NaN :+ NaN | 0.0 :+ 0.0 | 0 | |
0 ** (0:+1) | NaN :+ NaN | NaN :+ NaN | indeterminate | |
nan:+0 ** 0 | NaN :+ NaN | 1.0 :+ 0.0 | indeterminate | |
1 ** inf:+0 | NaN :+ NaN | NaN :+ NaN | indeterminate | |
1 ** nan:+0 | NaN :+ NaN | NaN :+ NaN | indeterminate | |
2 ** (inf:+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
2 ** (inf:+0) | NaN :+ NaN | NaN :+ NaN | infinity | |
2 ** ((-inf):+0) | NaN :+ NaN | NaN :+ NaN | 0 | |
2 ** ((-1000000000):+1) | 0.0 :+ 0.0 | 0.0 :+ 0.0 | ..... | |
(0.8:+0.8) ** inf:+0 | NaN :+ NaN | NaN :+ NaN | complex infinity | |
(0.8:+0.6) ** inf:+0 | NaN :+ NaN | NaN :+ NaN | indeterminate | |
(0.8:+0.4) ** inf:+0 | NaN :+ NaN | NaN :+ NaN | 0 | |
(0.8:+0.8) ** (-inf):+0 | NaN :+ NaN | NaN :+ NaN | 0 | |
(0.8:+0.6) ** (-inf):+0 | NaN :+ NaN | NaN :+ NaN | indeterminate | |
(0.8:+0.4) ** (-inf):+0 | NaN :+ NaN | NaN :+ NaN | complex infinity | |
(0.8:+0.8) ** (inf:+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
(0.8:+0.6) ** (inf:+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
(0.8:+0.4) ** (inf:+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
(0.8:+0.8) ** ((-inf):+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
(0.8:+0.6) ** ((-inf):+1) | NaN :+ NaN | NaN :+ NaN | Note 1 | |
(0.8:+0.4) ** ((-inf):+1) | NaN :+ NaN | NaN :+ NaN | Note 1 |
Note 1: Wolfram Alpha simplifies (infinity+i) to infinity, and as a result it misses the point of the indicated examples.
Note 2: The ghc result for -1**1/2 is somewhat unexpected due to the implementation involving complex log and an intermediate value of pi*i which of course is not represented with perfect accuracy. If this is a problem, it is not the one we are addressing here.
Changed 5 years ago by
Attachment: | complex_power.patch added |
---|
comment:26 Changed 5 years ago by
I agree with Scott Turner, we can leave x**infinity
for another day.
comment:27 Changed 5 years ago by
Cc: | ekmett added |
---|---|
Status: | infoneeded → patch |
yalas added a new patch, which needs to be reviewed.
comment:28 Changed 5 years ago by
Scott, that is an amazingly thorough summary.
Yalas, the patch looks sane to me.
It also just helped me catch some issues in my own automatic differentiation code, so that is a sign it is a good thing to have right in base
.
comment:29 Changed 5 years ago by
Milestone: | → 7.10.1 |
---|
I'm not reading the details here, but it sounds as if this is ready for 7.10.1, right?
Simon
comment:31 Changed 5 years ago by
One minor concern: will the optimizer do the right thing checking exp_re < 0
and exp_re > 0
, or will we be better off with case compare 0 exp_re ...
?
comment:32 Changed 5 years ago by
This is almost closed, but I forgot to add a note to the release notes. Incoming soon.
comment:35 Changed 5 years ago by
@david, unless theres benchmarks showing a substantial difference, probably doesnt matter
@thomie, why's it (re)marked new?
comment:37 Changed 5 years ago by
Milestone: | 7.10.1 → 7.12.1 |
---|
Moving to 7.12.1 milestone; if you feel this is an error and should be addressed sooner, please move it back to the 7.10.1 milestone.
comment:38 Changed 5 years ago by
Replying to thomie:
Because it broke validate, see comment:33.
Can I do something to fix it? And why did it brake validate?
comment:39 follow-up: 40 Changed 5 years ago by
Hacking up your patch to change the shadowing x
and two shadowing y
s should do it. e.g.
- ((re:+im), y) - | (isInfinite re || isInfinite im) -> case y of + ((re:+im), y') + | (isInfinite re || isInfinite im) -> case y' of - (x, y) -> exp (log x * y) + (x', y') -> exp (log x' * y')
comment:40 Changed 5 years ago by
Replying to ekmett:
Hacking up your patch to change the shadowing
x
and two shadowingy
s should do it. e.g.- ((re:+im), y) - | (isInfinite re || isInfinite im) -> case y of + ((re:+im), y') + | (isInfinite re || isInfinite im) -> case y' of - (x, y) -> exp (log x * y) + (x', y') -> exp (log x' * y')
Or just without duplicating them:
- ((re:+im), y) + ((re:+im), _) - (x, y) -> exp (log x * y) + _ -> exp (log x * y)
comment:41 Changed 5 years ago by
And we can simplify it a bit:
x ** y = case (x,y) of (_ , (0:+0)) -> 1 :+ 0 ((0:+0), (exp_re:+_)) -> case compare exp_re 0 of GT -> 0 :+ 0 LT -> inf :+ 0 EQ -> nan :+ nan ((re:+im), (exp_re:+_)) | (isInfinite re || isInfinite im) -> case compare exp_re 0 of GT -> inf :+ 0 LT -> 0 :+ 0 EQ -> nan :+ nan | otherwise -> exp (log x * y) where inf = 1/0 nan = 0/0
Changed 5 years ago by
Attachment: | complex_power.2.patch added |
---|
comment:43 Changed 5 years ago by
Status: | new → patch |
---|
That solves
0**nonzero
.0**0
is controversial but often 1 ( https://en.wikipedia.org/wiki/0^0 ) and in Double is 1:We probably want
0**0
to be consistent between Double and Complex Double.