Should complex numbers be closed at infinity by default?

[scheinerman]

Dividing a nonzero value by 0 should give an infinite result. This works fine for floating point numbers, but not for complex:

julia> 1. / 0.
Inf

julia> 1. / (0.+0.im)
NaN + NaN*im

Similarly, dividing by infinity (with a finite numerator) should give 0. This isn't working entirely properly:

julia> 1/Inf
0.0

julia> 1/(Inf*im)
0.0 - 0.0im

julia> 1/(Inf + 0*im)
0.0 - 0.0im

julia> 1/(Inf + Inf*im)
NaN + NaN*im         # this is wrong

Finally, in my opinion, there should be only a single ComplexInf value and not allow this:

julia> a = 3+Inf*im
3.0 + Inf*im

julia> b = 4+Inf*im
4.0 + Inf*im

julia> a==b
false

Both a and b should have evaluated to a common ComplexInf value. (We may also need a ComplexNaN.) Separate +Inf and -Inf is fine for Real values.

[andreasnoack]

Probably related: #5234

[jiahao]

There are two issues. The first is entirely a duplicate of #5234. The second is not a problem with how complex numbers are implemented, but rather a question of rather we want complex numbers to be closed at infinity by default.

[jiahao]

It is also good to cross-reference relevant threads on the mailing list.

Ref: julia-users.

[scheinerman]

I've posted an idea for how to handle here:
https://github.com/scheinerman/RiemannComplexNumbers.jl.git

[simonbyrne]

I think that both of the following are implementation bugs, and should be fixed:

julia> 1/(Inf + Inf*im)
NaN + NaN*im

julia> 1. / (0.+0.im)
NaN + NaN*im

From the C 2011 standard (with which we are generally trying to be compatible):

• if the first operand is a finite number and the second operand is an infinity, then the result of the / operator is a zero;
• if the first operand is a nonzero finite number or an infinity and the second operand is
  a zero, then the result of the / operator is an infinity.

@scheinerman I'm still not really convinced of the advantages of a single complex infinity: do you have a particular application in mind for which this would be better than the current behaviour?

[scheinerman]

Hi Simon,

It's a fair question: do I have an application in mind? Very short answer: No.

My perspective is that of a mathematician and how we might extend the (finite) complex numbers to include infinity. For the real numbers, extending with separate positive and negative infinities is fine. Having those two, and only those two, infinities for the complex numbers doesn't make sense. This leads us to a few reasonable choices. But for whatever we decide, having Inf+2im and Inf+3im as different values is not reasonable; they should compare as equal just like Inf+2 and Inf+3.

A single ComplexInf turns the complex plane into the Riemann sphere. This ComplexInf is the value returned by dividing any nonzero complex number by 0. Mathematically and aesthetically, I think this is the cleanest.

But if we wish to emulate +Inf and -Inf from the real case, we can imagine that every ray in the complex plane points to an infinite value (for the reals, there are only two ways a ray can point). Parallel rays point to the same infinity. So if we write complex numbers in polar coordinates, (r,theta) then the infinities are of the form (Inf, theta) where two complex infinities are equal iff their angles are the same modulo 2pi. Thus we have an entire circle of infinities added to the complex plane.

I also think that having a single ComplexNaN is a good idea. We could, of course use the existing NaN but that is of type Float64. I think the result of complex operations should always be a Complex type. This is what I have attempted to do here: [https://github.com/scheinerman/RiemannComplexNumbers.jl.git]

But I concede this: Efficiency and speed may well trump mathematical elegance. Having to do checks for every complex operation will slow things down. Folks that just want to do a bunch of complex arithmetic quickly in which infinities and/or NaN's are not going to crop up would not want this overhead. Hence my proposed solution is to create a package that a user may load that redefines Complex operations to behave well in the face of infinities and NaN's. Those who only want speed and for which this is an academic exercise simply wouldn't load the package, but those who are concerned about this can then ask for the more careful operations. And we could even have more than one way to do this: one with a single ComplexInfinity as I favor and another that would have an infinity for every parallel class of rays in the complex plane. (For those who like projective geometry, we could even have a third in which there is an infinity for every parallel class of lines.)

[eschnett]

There are probably arguments to be made either way. The Riemann sphere is very elegant, but treating complex numbers as a Cartesian pair of real numbers is also elegant.

I favour consistency between languages. Since e.g. C, C++, or maybe even Lisp define how infinities and nans should behave, we should follow suit. If we want to diverge, then this should be a large community-wide decision, and we may even want comments from the C and C++ standards committees to see why they decided their way.

The argument regarding "don't care about infinities and nans" is not valid. The IEEE standard defines very closely how these are handled, and by default, Julia follows this standard. We should attempt something equivalent for complex numbers, and not say "efficiency trumps everything for complex numbers". For example, there's always @fastmath or similar options for people who don't care about infinities and nans and signed zeros.

[jiahao]

In the julia-users thread referenced, I suggested yet another possibility where equality comparison could be redefined for complex numbers so that complex infinity could have multiple representations, but they would all compare equal.

==(z::Complex, w::Complex) = (real(z) == real(w)) & (imag(z) == imag(w)) #Argand plane

==(z::Complex, w::Complex) = if isinf(z) && isinf(w)
    return true #Riemann closure at infinity
else
    return (real(z) == real(w)) & (imag(z) == imag(w))
end

It might seem a little crazy, but redefining equality is actually not an unreasonable way to capture the notion of topological closure. This option just happens to be not usually possible in most languages.

Ref: JeffBezanson/phdthesis#4

[simonbyrne]

Hmm, this is an intriguing idea: the main downside would be the overhead of the branch.

The C standard seems fairly flexible when it comes to complex infinities: as above, it uses the phrase "an infinity" and "a zero" a lot. In case anyone is wondering, the C standard (like Julia) regards NaN + Inf*im as an infinity.

Any discussion of complex infinities cannot go without a reference to Kahan's article "Branch Cuts for Complex Elementary Functions or Much Ado About Nothing's Sign Bit".

We should probably also try to conform to ISO/IEC 10967-3:2006 "Language independent arithmetic -- Part 3: Complex integer and floating point arithmetic and complex elementary numerical functions" (free link from here), which specifies "general" behaviour of complex numbers in software (though it doesn't seem to say too much about infinities).

[jiahao]

@simonbyrne The language in the ISO standard mirrors much of the C99 spec. I found the C9X N557 draft more enlightening, as it was annotated with comments and discussions that were omitted from the final spec.

[simonbyrne]

Yes, I've looked briefly through other LIA specs and generally they are not all that enlightening, as they seem to be written so that almost all existing software would be conformant.

[scheinerman]

Hello again folks.

First, I was curious to see how C++ handles complex division by zero. So I wrote this:

#include <iostream>
#include <complex>
using namespace std;

main() {
  complex<double> zero = complex<double>(0,0);
  complex<double> one  = complex<double>(1,0);
  complex<double> eye  = complex<double>(0,1);
  cout << "1/0 = " << one/zero << endl;
  cout << "i/0 = " << eye/zero << endl;
  cout << "0/0 = " << zero/zero << endl;
  cout << "(1+i)/0 = " << (one+eye)/zero << endl;
}

and the output when run on my Mac is this

1/0 = (inf,nan)
i/0 = (nan,inf)
0/0 = (nan,nan)
(1+i)/0 = (inf,inf)

but when run on my Ubuntu linux machine I get something slightly different:

1/0 = (inf,-nan)
i/0 = (-nan,inf)
0/0 = (-nan,-nan)
(1+i)/0 = (inf,inf)

The negative nan values really surprised me.

So the fact that Julia evaluates (1+0im)/0 as Inf + NaN*im is consistent with C++.

Now Jiahao's idea of redefining == so all the manifestations of complex infinity merge together sounds good, but I think we'd need a bit more work because if z is any flavor of infinity, I'd want 1/z to be zero. This won't happen without more work:

julia> z = (1+0im)/0
Inf + NaN*im

julia> isinf(z)
true

julia> 1/z
NaN + NaN*im

By the way, this is something C++ gets right. Dividing (1,0) by (inf,nan) does yield (0,0).

[eschnett]

Re-defining equality to close complex numbers at infinity as one disadvantage:

julia> x=Inf;
julia> y=-x;
julia> x==y
false
julia> complex(x)==complex(y)
true

That is, numbers that compare unequal as Float then suddenly compare equal when converted to Complex.

[scheinerman]

Yup. We can’t have everything. One might argue (I won’t) that there should only be a single inf for float values too (one-point closure of the real number system).

Thinking about this more over the day leads me to think we should leave this to user preference. Provide a default that conforms to C++ or some other good standard, but then people can load a module to redefine the behavior to their preference/application.

By the way, when I redefine (say) + for Complex numbers, I get a warning that I’ve redefined + for Complex numbers … anyway to suppress that?

On Jan 18, 2015, at 5:12 PM, Erik Schnetter [email protected] wrote:

Re-defining equality to close complex numbers at infinity as one disadvantage:

julia> x=Inf;
julia> y=-x;
julia> x==y
false
julia> complex(x)==complex(y)
true
That is, numbers that compare unequal as Float then suddenly compare equal when converted to Complex.

—
Reply to this email directly or view it on GitHub #9790 (comment).

[jiahao]

There is no reason to expect that the real number line can be closed the same way as the complex numbers can.

Agreed that closure should be user opt-in behavior.

[jiahao]

Unfortunately it looks like division by complex zero has to be special cased; none of the commonly used algorithms compute the complex infinity correctly. (The last two are taken directly from the base library.)

#Naive
function naive(z::Complex, w::Complex)
    a, b = reim(z)
    c, d = reim(w)
    den = c*c + d*d
    Complex( (a*c + b*d)/den, (b*c - a*d)/den)
end

#Smith, 1962
function smith(z::Complex, w::Complex)
    a₀, b₀ = reim(z)
    a₁, b₁ = reim(w)

    if abs(a₁) ≥ abs(b₁)
        c = b₁/a₁
        d = a₁ + b₁*c
        Complex((a₀ + b₀*c)/d, (b₀ - a₀*c)/d)
    else
        c = a₁/b₁
        d = a₁*c + b₁
        Complex((a₀*c + b₀)/d, (b₀*c - a₀)/d)
    end
end

#Priest, 2004
function priest{T<:Real}(z::Complex{T}, w::Complex{T})
    a₀, b₀ = reim(z)
    a₁, b₁ = reim(w)

    s = (a₁*a₁ + b₁*b₁)^(-1.5) #Optimal choice; in practice, scale factor is chosen iteratively
    a₁′ = s * a₁
    b₁′ = s * b₁
    t = 1/(a₁′*a₁′ + b₁′*b₁′)
    a₁′′ = s * a₁′
    b₁′′ = s * b₁′

        x = (a₀ * a₁′′ + b₀ * b₁′′ ) * t
        y = (b₀ * a₁′′ - a₀ * b₁′′ ) * t

    Complex(x, y)
end

#Kahan
function kahan{T<:Real}(a::Complex{T}, b::Complex{T})
    are, aim = reim(a)
    bre, bim = reim(b)
    if abs(bre) <= abs(bim)
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        Complex((are*r + aim)/den, (aim*r - are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        Complex((are + aim*r)/den, (aim - are*r)/den)
    end
end

# Baudin, 2012
# arXiv:1210.4539
function baudin{T<:Real}(z::Complex{T}, w::Complex{T})
    a, b = reim(z); c, d = reim(w)
    half = 0.5
    two = 2.0
    ab = max(abs(a), abs(b))
    cd = max(abs(c), abs(d))
    ov = realmax(a)
    un = realmin(a)
    ϵ = eps(T)
    bs = two/(ϵ*ϵ)
    s = 1.0
    ab >= half*ov  && (a=half*a; b=half*b; s=two*s ) # scale down a,b
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    ab <= un*two/ϵ && (a=a*bs; b=b*bs; s=s/bs      ) # scale up a,b
    cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs      ) # scale up c,d
    abs(d)<=abs(c) ? ((p,q)=robust_cdiv1(a,b,c,d)  ) : ((p,q)=robust_cdiv1(b,a,d,c); q=-q)
    return Complex(p*s,q*s) # undo scaling
end
function robust_cdiv1{T<:Real}(a::T, b::T, c::T, d::T)
    r = d/c
    t = 1.0/(c+d*r)
    p = robust_cdiv2(a,b,c,d,r,t)
    q = robust_cdiv2(b,-a,c,d,r,t)
    return p,q
end
function robust_cdiv2{T<:Real}(a::T, b::T, c::T, d::T, r::T, t::T)
    if r != 0
        br = b*r
        return (br != 0 ? (a+br)*t : a*t + (b*t)*r)
    else
        return (a + d*(b/c)) * t
    end
end

@show naive(Complex(1.0, 0.0), Complex(0.0, 0.0)) #NaN + NaN*im
@show smith(Complex(1.0, 0.0), Complex(0.0, 0.0)) #NaN + NaN*im
@show priest(Complex(1.0, 0.0), Complex(0.0, 0.0)) #NaN + NaN*im
@show kahan(Complex(1.0, 0.0), Complex(0.0, 0.0)) #NaN + NaN*im
@show baudin(Complex(1.0, 0.0), Complex(0.0, 0.0)) #NaN + NaN*im

[simonbyrne]

Yes, we would need an explicit branch: there is some sample code for complex division in the C 1999 and 2011 standards (they're slightly different), and they both include an explicit branch to check for a zero denominator, or either argument an infinity.

[toivoh]

@scheinerman: When you redefine + by redefining the method for it that adds eg two complex numbers, you redefine it for all code and not just the current module, so that might not be a good way to go about it. You can have a local+ function instead, but it's a bit tricky at the moment.

[scheinerman]

@toivoh : Thanks. Actually, that's exactly what I had in mind. A user who wants to work with alternative methods for Complex numbers would load a module and from then on, the new + (and other operations) would hold sway. There could be a module that has no infinity/NaN checking (for those who want to go fast), another (as I would like) that would have a single ComplexInf and a single ComplexNaN, and perhaps others that have further ways to deal with infinity for complex numbers. The core, default Julia Complex type can then be some industry standard.

So that's what I'm trying to do (though I'm not a pro!) with my RiemannComplexNumbers module. I think it works, but the initial burst of warnings is unsightly. I've run into the same problem in another project in which I want to redefine show for Set objects. Here's what I've done there:

julia> X = Set([1,2,3])
Set{Int64}({2,3,1})

julia> using ShowSet
Warning: Method definition show(IO,Set{T}) in module Base at set.jl:10 overwritten in module ShowSet at /Users/ers/.julia/v0.3/ShowSet/src/ShowSet.jl:25.
Warning: Method definition show(IO,IntSet) in module Base at intset.jl:176 overwritten in module ShowSet at /Users/ers/.julia/v0.3/ShowSet/src/ShowSet.jl:26.

julia> X
{1,2,3}

I wonder if there's a way to (temporarily) suppress those warnings.

[eschnett]

Instead of redefining complex arithmetic, could you redefine the complex number type? I'm thinking of Complex = ComplexRiemann in the beginning of your module.

Overwriting a routine is global. (I assume this is what the warning is about.) If you then call a library that was designed with the "other" kind of complex numbers, things may break. If you instead overwrite things only in your module, everything should be safe.

[JeffBezanson]

Complex division is already so slow that I wouldn't worry about adding operations to it.

[mattcbro]

What happens if the infinities and NaNs occur as part of a linear algebra computation that is presumably rendered in a Blas or Lapack library? Will you end up with inconsistencies?

From my point of view if all other matters are tied performance should win. After all performance is one of the great notable features of Julia.

[scheinerman]

I'm not sure what happens when a Blas or Lapack function is called ... can you point me to some standard ones that I can try out?

As to the performance issue, I think the resolution is to put the choice in the user's hands. If we're sure that infinities and NaN's won't occur, then speed wins the day. But I also think being correct is important. Consider this

julia> z = (1+0im)/0
Inf + NaN*im

julia> isnan(z)
true

I would argue this is just plain wrong. After using RiemannComplexNumbers we get this:

julia> z = (1+0im)/0
ComplexInf

julia> isnan(z)
false

and that's correct.

[tkelman]

If this hasn't had any discussion or activity in 16 months, I don't think it's release-blocking.

[scheinerman]

I agree that performance top priority, but we have this inconsistency in division by zero:

julia> (1+0im)/0
Inf + NaN*im

julia> (1+0im)/(0+0im)
NaN + NaN*im

Perhaps that can be fixed in a way that doesn't impact performance?

[StefanKarpinski]

This is certainly not a release blocker at this point.

We've worked very hard to make our == transitive, so I consider violating that with something like @eschnett's above example a non-starter, i.e. Inf == Complex(Inf) == -Inf != Inf. In that case Complex(Inf) would have to be considered a different infinity entirely, not equal to either of the real infinities. Personally, I think we should not have Riemannian complex behavior in Base Julia, but we should certainly make it possible to opt into such behavior with a package, which @scheinerman's package does, albeit with annoying warnings.

[PallHaraldsson]

"But for whatever we decide, having Inf+2im and Inf+3im as different values is not reasonable"

I only read this far. I was actually looking into implementing Complex with polar coordinates. There are Pros and cons (maybe I get around the cons..). At least your issue is then non-existent?

[scheinerman]

Yes. There are two natural ways to complete the complex numbers either with a single "Complex Infinity" (my preference) or a "line at infinity" (having a different point at infinity for each ray or line out of the origin). In polar coordinates, then both Inf+2im and Inf+3im would have r=Inf and theta=0, and therefore be the same value.

[StefanKarpinski]

I don't think we're going to change from our current semantics for infinite complex values. The reals are embedded in the complexes as far as equality is concerned – as the values for which the imaginary part is zero. Some of the corner cases should probably be fixed, but we should have individual issues for those.

[StefanKarpinski]

The only issue here is division by zero, which could be argued to be a bug and/or undefined.

[StefanKarpinski]

The mixed type division should be adjusted to be in line with complex-complex division.