Fix behaviour of + and - with rational infinity

[TotalVerb]

The current behaviour of rational infinities under + and - differs from both floats and mathematical intuition:

julia> 1//0 + 1//0
ERROR: DivideError: integer division error
 [inlined code] from ./rational.jl:19
 in +(::Rational{Int64}, ::Rational{Int64}) at ./rational.jl:179
 in eval(::Module, ::Any) at ./boot.jl:236

I have created a fix and tests to make sure rational infinities behave like floating point ones. This is not a very elegant fix, but I could not figure out a better way.

[tkelman]

better to be as specific as possible

[TotalVerb]

What kind of exception is most appropriate? Currently it throws either ArgumentError or DivideError depending on what goes wrong first, and I think that could be improved (if the overhead of a try-catch is acceptable, of course). I am thinking that DomainError makes most sense, since if the result would be NaN for floats that would mean the inputs are not in the domain of the function.

[TotalVerb]

Upon further inspection, it looks like DivideError is currently thrown everywhere, because almost every rational number exception occurs when divgcd is called with two zeroes. The only exception is the constructor for Rational, which does a sort of divgcd operation without actually calling divgcd. It throws ArgumentError instead.

I looked at the uses for divgcd and in almost each case, it was used immediately before constructing a rational number. So I think it is safe to simply have divgcd throw an ArgumentError which will propagate through all the other operations.

The few exceptions are:

julia> 0//1 ÷ 0

This now results in ArgumentError: operation results in invalid rational, even though the return type isn't rational. Whereas

julia> 1//1 ÷ 0

still results in DivideError. Same goes for 0 ÷ 0//1 and 1 ÷ 0//1, which now return different errors.

They were also returning different errors before (as in, different stack frame causing the error), but this time the error type has changed also.

[TotalVerb]

Here's another conundrum.

julia> 1//0 % 0//1

DivideError or ArgumentError? We're dividing by 0, but it's not an integer in this case.

julia> 1//0 % 1//0

DivideError or ArgumentError? The docs only say DivideError occurs when dividing by 0 or typemin by -1. Here the division is by infinity.

Maybe all operations on rationals should be DivideError. This would require the Rational constructor itself to throw DivideError instead of ArgumentError.

[StefanKarpinski]

Just as an aside, the fact that it's often so difficult to figure out the "right" exception type is one of the things that makes them sort of suspicious to me.

[cstjean]

I'm surprised that +-infinity is allowed for rationals without a corresponding +-0. Why should 1//0 yield +infinity and not -infinity? inv(inv(-1//0)) == -1//0 is false, because the LHS evaluates to +infinity.

[pkofod]

I'm surprised that +-infinity is allowed for rationals without a corresponding +-0.

A compromise I guess. It's to mimick floats (double) , such that conversion between floats and rationals are lossless. Then we can convert between -Inf (float) and -1//0, and Inf and 1//0. It is true we cannot convert -0.0 to -0 //1 and back to -0.0, so there is loss of information there. To fix it, I guess we would have to use positive integers in the denom and num, and then add a sign-field.

[TotalVerb]

For what it's worth, I think that the signed 0 is a good idea. But using an additional byte of memory just for a signed zero is excessive in my opinion.

[nalimilan]

The loss on conversion isn't specific to rationals:

julia> Int(-0.0)
0

Though it looks like it could be a avoided for rationals: instead adding a field for the sign, which would be wasteful, why not use the convention that 0//x is converted to sign(x) * 0.0? Not sure that's worth it, though.

[pkofod]

why not use the convention that 0//x is converted to sign(x) * 0.0?

What do you mean? The numerator holds the sign such that -1//0 is "minus infinity" and 1//0 is plus infinity. In y//x we always have sign(x) == 1. If we chose to have the denominator hold the sign, then we can distinguish -0//1 from 0//1, but no longer infinity. Do you propose to sometimes hold the sign in num and sometimes in den ?

[TotalVerb]

Allowing negative numbers in the denominator only if the numerator numerator is zero is a bit of a hack. We could dispense of the normalization for rationals, and only compute gcds when required (on overflow or when displaying). That could make rational arithmetic faster also.

[nalimilan]

Yes, that's a hack and I'm not recommending it. I just mentioned it as an alternative to adding a separate field, which is clearly overkill.

[StefanKarpinski]

I think what @nalimilan was suggesting is to have

Rational(+0.0) === 0//+1
Rational(-0.0) === 0//-1

Currently we normalize rationals so that the denominator is always non-negative – this would change that.

[TotalVerb]

I have a solution to the error problem that I like:

• If the result of the operation would be an integer, then the problem is assumed to be an integer division, and is DivideError.
• If the result of the operation would be a rational, then an integer division may not have occurred at all, so the correct error is ArgumentError.

I updated the tests with these semantics.

---

Edit: Actually, on second thought, these semantics are incorrect (and the tests for them are incomplete). They fail to cover the rem(1//0, 0//1) case... though this returns a rational, it is clearly a division by 0 that is the problem, and not the creation of an invalid 0//0 rational.

I feel like the "correct" semantics are very similar to those in the code right now... if 0//0 is constructed, then the result is ArgumentError, and if an integer division by 0 occurs, then the result is DivideError. But these semantics are a bit hard to think about unless one reasons about how exactly the computation is done.

[TotalVerb]

This is still an issue on Julia 0.6, and regrettably the answer to "what's the right exception type to throw" is holding it up.

[StefanKarpinski]

Throwing an inconsistent exception seems like a fairly separate issue. The most important question is exception or not.

[TotalVerb]

For what case? There is a pretty good argument that 1//0 + 1//0 should not throw an exception.

[StefanKarpinski]

Agree 1//0 + 1//0 should not be an error, nor should -1//0 + -1//0.

[StefanKarpinski]

You should probably rebase this.

[TotalVerb]

I tried to do that using the GitHub interface; I suppose that failed CI. I'll rebase it properly tomorrow.

[TotalVerb]

I've rebased this. I modified a test to use @testsets for convenience reasons; it was difficult otherwise to track down a failing test. I can revert that part if too unrelated.

[TotalVerb]

Is this good to go? It's a bugfix so should be OK for 0.6.

[TotalVerb]

Bump.

[StefanKarpinski]

Looks good to me. I left some comments, but I think you should ignore them; we can make those kinds of changes at some point in the future since they're largely stylistic.

[StefanKarpinski]

Would it make more sense to promote to a common type first:

divgcd(x::Integer, y::Integer) = divgcd(promote(x, y)...)

function divgcd(x::T, y::T) where T<:Integer
    if !iszero(x) || !iszero(y)
        g = gcd(x, y)
        x ÷= g
        y ÷= g
    end
    return x, y
end

[TotalVerb]

It's possibly intentional that things are done this way, since div has special cases for unsigned/signed combinations.

julia> @which 0x1 ÷ -1
div(x::Unsigned, y::Union{Int128, Int16, Int32, Int64, Int8}) in Base at int.jl:160

This behavior might require reviewing for desirability/consistency.

[StefanKarpinski]

Likewise, maybe promote first, but then again, maybe that can be another PR.

[nalimilan]

I'm really not the best person to review this, I never use rationals. My suggestion above was just an idea...

[musm]

Fixed on master (I can't find the specific PR where these were fixed but running the added test here passes, as well as the error in the PR's first message)