Opened 6 years ago

Closed 5 years ago

# Data.Complex shouldn't use default implementation of (**)

Reported by: Owned by: jjaredsimpson low 8.0.1 Prelude 7.6.3 ekmett Unknown/Multiple Unknown/Multiple Incorrect result at runtime

### 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)

### comment:1 Changed 6 years ago by isaacdupree

That solves 0**nonzero. 0**0 is controversial but often 1 ( https://en.wikipedia.org/wiki/0^0 ) and in Double is 1:

> 0**0 :: Double
1.0

We probably want 0**0 to be consistent between Double and Complex Double.

### comment:2 Changed 6 years ago by Scott Turner

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 yalas

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

### comment:4 Changed 5 years ago by yalas

Status: new → patch

### comment:5 Changed 5 years ago by Scott Turner

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 yalas

(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 Scott Turner

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 in reply to:  7 Changed 5 years ago by yalas

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.

Last edited 5 years ago by yalas (previous) (diff)

### comment:9 Changed 5 years ago by yalas

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

x ** (1:+0) = x

But I'm not sure. Is it necessary? This function already has 3 definitions.

### comment:10 Changed 5 years ago by yalas

By the way, do we need some tests?

### comment:11 Changed 5 years ago by yalas

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.

Last edited 5 years ago by yalas (previous) (diff)

### comment:12 Changed 5 years ago by yalas

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?

Last edited 5 years ago by yalas (previous) (diff)

### comment:13 Changed 5 years ago by Scott Turner

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 Scott Turner

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 Scott Turner

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 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.

### comment:17 in reply to:  16 Changed 5 years ago by yalas

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?

Last edited 5 years ago by yalas (previous) (diff)

### comment:18 follow-ups:  19  20 Changed 5 years ago by Scott Turner

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 in reply to:  18 Changed 5 years ago by yalas

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 in reply to:  18 Changed 5 years ago by yalas

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.

Last edited 5 years ago by yalas (previous) (diff)

### comment:21 follow-ups:  22  23 Changed 5 years ago by 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.

### comment:22 in reply to:  21 Changed 5 years ago by yalas

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.

(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 in reply to:  21 Changed 5 years ago by yalas

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 thoughtpolice

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 Scott Turner

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.

### comment:26 Changed 5 years ago by yalas

I agree with Scott Turner, we can leave x**infinity for another day.

### comment:27 Changed 5 years ago by thomie

Cc: ekmett added infoneeded → patch

yalas added a new patch, which needs to be reviewed.

### comment:28 Changed 5 years ago by ekmett

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 simonpj

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:30 Changed 5 years ago by Austin Seipp <austin@…>

base: Fix (**) instance for Data.Complex (#8539)

Reviewed-by: Edward Kmett <ekmett@gmail.com>
Authored-by: Yalas, Scott Turner

Signed-off-by: Austin Seipp <austin@well-typed.com>

### comment:31 Changed 5 years ago by dfeuer

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 thoughtpolice

This is almost closed, but I forgot to add a note to the release notes. Incoming soon.

### comment:33 Changed 5 years ago by Austin Seipp <austin@…>

Revert "base: Fix (**) instance for Data.Complex (#8539)"

This broke validate due to name shadowing warnings.

This reverts commit 1f6b1ab4b6d7203481bfaf374b014972f7756fb2.

### comment:34 Changed 5 years ago by thomie

Status: patch → new

Needs work.

### comment:35 Changed 5 years ago by carter

@david, unless theres benchmarks showing a substantial difference, probably doesnt matter

@thomie, why's it (re)marked new?

### comment:36 follow-up:  38 Changed 5 years ago by thomie

Because it broke validate, see comment:33.

### comment:37 Changed 5 years ago by thoughtpolice

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 in reply to:  36 Changed 5 years ago by yalas

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 ekmett

Hacking up your patch to change the shadowing x and two shadowing ys 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')
Last edited 5 years ago by ekmett (previous) (diff)

### comment:40 in reply to:  39 Changed 5 years ago by lukyanov

Hacking up your patch to change the shadowing x and two shadowing ys 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 lukyanov

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

Even better.

### comment:43 Changed 5 years ago by thomie

Status: new → patch

### comment:44 Changed 5 years ago by Austin Seipp <austin@…>

base: Fix (**) implementation for Data.Complex

See the extensive discussion in #8539.

Signed-off-by: Austin Seipp <austin@well-typed.com>

### comment:45 Changed 5 years ago by thoughtpolice

Resolution: → fixed patch → closed

Merged, thanks!

### comment:46 Changed 4 years ago by thoughtpolice

Milestone: 7.12.1 → 8.0.1

Milestone renamed

Note: See TracTickets for help on using tickets.