Opened 5 months ago

Last modified 7 days ago

#8539 new bug

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

Reported by: jjaredsimpson Owned by:
Priority: low Milestone:
Component: Prelude Version: 7.6.3
Keywords: Cc:
Operating System: Unknown/Multiple Architecture: Unknown/Multiple
Type of failure: Incorrect result at runtime Difficulty: Easy (less than 1 hour)
Test Case: Blocked By:
Blocking: Related Tickets:

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 (2)

complex_power.path (520 bytes) - added by yalas 7 days ago.
complex_power.patch (520 bytes) - added by yalas 7 days ago.

Download all attachments as: .zip

Change History (5)

comment:1 Changed 5 months 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 5 months 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 7 days 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

Changed 7 days ago by yalas

Changed 7 days ago by yalas

Note: See TracTickets for help on using tickets.