Support complex rational infinity

[eschnett]

Rational numbers (Rational{Int}) support infinity in the sense that inv(0//1) is handled well. This is not the case for complex rational numbers, i.e. inv(Complex(0//1)) is not handled; instead, it reports an ERROR: DivideError: integer division error.

I understand the reason for this. In the standard conformal completion of the complex number plan, there is just a single point at infinity. This is different from Julia's choice for completing the real number line, where there are two infinities (-1//0 and 1//0).

I suggest to allow complex infinity for rational numbers. It doesn't really matter for my application how complex infinity is handled – it could be a single point, or it could be handled the same as for complex floating-point numbers – but I want to avoid thrown an exception.

Why?

The problem I'm encountering is with Julia's linear algebra. To test certain algorithms I don't want to have to deal with floating-point errors, and I thus use Rational{BigInt} and Complex{Rational{BigInt}} as element types. This works well in general, but the corner cases are all haywire.

For example, let's consider inverting each element of a sparse vector: map(inv, SparseVector(...)). The sparse vector map function check whether the result of inv applied to the element type is zero. If the element type is Complex{Rational{BigInt}}, then inv throws an exception, even if none of the elements are zero (i.e. if the sparse vector happens to be dense):

using LinearAlgebra
using SparseArrays
julia> inv.(sparsevec([1], [1//1]))
1-element SparseVector{Rational{Int64}, Int64} with 1 stored entry:
  [1]  =  1//1
julia> inv.(sparsevec([1], [Complex(1//1)]))
ERROR: DivideError: integer division error

Note that this exception is spurious: The result is well-defined, Julia should return sparsevec([1], [Complex(1//1)]) again.

Handling this (and many similar cases in other types/functions) is quite tedious. One would need to add new methods for Base.inv, as in

function Base.inv(D::Diagonal{T,SparseVector{T,I}}) where {T,I}
    if length(D.diag.nzval) == D.diag.n
        # Omit check `iszero(inv(zero(T)))`
        return Diagonal(SparseVector(D.diag.n, D.diag.nzind, map(inv, D.diag.nzval)))
    else
        # Fall back to original function
        return Diagonal(map(inv, D.diag))
    end
end

I was considering submitting respective pull requests to Julia (for quite a few linear algebra algorithms) – inv, ldiv!, rdiv!, map, etc., for Diagonal, SparseVector, SparseArrayCSC, and more. Essentially, everywhere the value zero(T) is used in a way that could involve an inverse or a division, that respective check needs to be delayed.

Allowing some kind of complex rational infinity instead would simplify this a lot.

[JeffreySarnoff]

+1