Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

-im*Inf should probably be equal to -Inf*im #9531

Open
dlfivefifty opened this issue Jan 1, 2015 · 7 comments
Open

-im*Inf should probably be equal to -Inf*im #9531

dlfivefifty opened this issue Jan 1, 2015 · 7 comments
Labels
complex Complex numbers equality Issues relating to equality relations: ==, ===, isequal maths Mathematical functions

Comments

@dlfivefifty
Copy link
Contributor

This behaviour is surprising and probably wrong:

julia> -Infim
-0.0 - Inf
im

julia> -imInf
NaN - Inf
im

@Ismael-VC
Copy link
Contributor

The issue is that -im is actually 0 - 1im, and the product of Inf by 0 is NaN:

julia> 0Inf
NaN

julia> -im
0 - 1im

julia> ans * Inf
NaN - Inf*im

@jiahao
Copy link
Member

jiahao commented Jan 1, 2015

This is a special case of well-known problems with nonfinite complex floating point operations:

#*(z::Complex, w::Complex) = Complex(real(z) * real(w) - imag(z) * imag(w),
#                                    real(z) * imag(w) + imag(z) * real(w))
julia> Complex(5.0, Inf) * Complex(0.0, 1.0)
-Inf + NaN*im

Arguably, the issue here is that we promote integral 0 to floating-point 0.0, then multiply by Inf, and IEEE rules define 0.0 * Inf == NaN.

One alternative is to try to make 0*Inf == 0; this is not currently the case:

julia> 0*Inf
NaN

@jiahao
Copy link
Member

jiahao commented Jan 1, 2015

  1. Another thought to consider:

    julia> -(im*Inf)
    -0.0 - Inf*im
    
    julia> (-im)*Inf
    NaN - Inf*im
  2. The Handbook of Floating-Point Arithmetic, p 144 suggests using Karatsuba's algorithm, which for Complex*Real can be simplified to:

    function times(a::Complex, a1::Real) #Karatsuba
         a0, b0 = reim(a)
         p1 = (a0+b0)*(a1) #a1+b1
         p2 = a0*a1
        #p3 = b0*b1
         a = p2 #-p3
         b = (p1-p2) #-p3
         Complex(a, b)
       end

However, the answer is even worse for this example:

julia> times(-im, Inf)
p1 = (a0 + b0) * a1 = -Inf
p2 = a0 * a1 = NaN
a = p2 = NaN
b = p1 - p2 = NaN
NaN + NaN*im

@jiahao
Copy link
Member

jiahao commented Jan 1, 2015

Another thought was to do complex multiplication in polar form, but there is no floating point number x near pi/2 such that cos(x) == 0.0 exactly, and the answer is also not useful:

julia> type ComplexPolar
         r
         θ
       end

julia> *(a::ComplexPolar, b::ComplexPolar) = ComplexPolar(a.r*b.r, a.θ+b.θ)
* (generic function with 120 methods)

julia> Base.convert(::Type{Complex},a::ComplexPolar) = Complex(a.r*cos(a.θ), a.r*sin(a.θ));
convert (generic function with 493 methods)

julia> Base.convert(::Type{ComplexPolar},a::Complex) = ComplexPolar(abs(a), angle(a))
convert (generic function with 494 methods)

julia> a = convert(ComplexPolar, -im)
ComplexPolar(1.0,-1.5707963267948966)

julia> b = convert(ComplexPolar, Complex(Inf,0))
ComplexPolar(Inf,0.0)

julia> a*b
ComplexPolar(Inf,-1.5707963267948966)

julia> convert(Complex, ans)
Inf - Inf*im

julia> cos(-1.5707963267948966)
6.123233995736766e-17

julia> cos(prevfloat(-1.5707963267948966))
-1.6081226496766366e-16

@dlfivefifty
Copy link
Contributor Author

May store theta as a number between [-1,1) that is multiplied by pi?

Sent from my iPad

On 2 Jan 2015, at 5:37 am, Jiahao Chen notifications@github.com wrote:

Another thought was to do complex multiplication in polar form, but there is no floating point number x near pi/2 such that cos(x) == 0.0 exactly, and the answer is also not useful:

julia> type ComplexPolar
r
θ
end

julia> _(a::ComplexPolar, b::ComplexPolar) = ComplexPolar(a.r_b.r, a.θ+b.θ)

  • (generic function with 120 methods)

julia> Base.convert(::Type{Complex},a::ComplexPolar) = Complex(a.r_cos(a.θ), a.r_sin(a.θ));
convert (generic function with 493 methods)

julia> Base.convert(::Type{ComplexPolar},a::Complex) = ComplexPolar(abs(a), angle(a))
convert (generic function with 494 methods)

julia> a = convert(ComplexPolar, -im)
ComplexPolar(1.0,-1.5707963267948966)

julia> b = convert(ComplexPolar, Complex(Inf,0))
ComplexPolar(Inf,0.0)

julia> a*b
ComplexPolar(Inf,-1.5707963267948966)

julia> convert(Complex, ans)
Inf - Inf*im

julia> cos(-1.5707963267948966)
6.123233995736766e-17

julia> cos(prevfloat(-1.5707963267948966))
-1.6081226496766366e-16

Reply to this email directly or view it on GitHub.

@jiahao
Copy link
Member

jiahao commented Jan 1, 2015

If you do that, the answer is still NaN - Inf*im, because the real part is Inf * cospi(0.5) == Inf * 0.0 == NaN.

@dlfivefifty
Copy link
Contributor Author

Maybe this is an argument for an imaginary number type?

Sent from my iPad

On 2 Jan 2015, at 7:42 am, Jiahao Chen notifications@github.com wrote:

If you do that, the answer is still NaN - Inf*im, because the real part is Inf * cospi(0.5) == Inf * 0.0 == NaN.


Reply to this email directly or view it on GitHub.

@jiahao jiahao added the complex Complex numbers label Feb 1, 2015
@StefanKarpinski StefanKarpinski added equality Issues relating to equality relations: ==, ===, isequal maths Mathematical functions labels May 10, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
complex Complex numbers equality Issues relating to equality relations: ==, ===, isequal maths Mathematical functions
Projects
None yet
Development

No branches or pull requests

4 participants