Faster Rationals by avoiding unnecessary divgcd

[Liozou]

This PR prevents divgcd from being unnecessarily computed when creating a Rational with numerators and denominators that are known to be coprime. This is entirely non-breaking.

Addresses #11522

The history

I think if we really wanted to tweak performance of rationals, one way would be to have an option to not perform the gcd in the constructor when it isn't required, e.g. for Rational(::Int) or after a multiplication or division. (@simonbyrne)
That seems like it would be an improvement. I guess that would require having a way of constructing Rational values, bypassing the gcd. (@StefanKarpinski)

(#8672)

The implementation

I add two constructors to Rational that take the numerator, the denominator, and either Val(false) or Val(true). The original constructor is left unchanged (that would be breaking).

The two new constructors do not reduce the fraction to its irreducible form. Since the Rational infrastructure relies on that fact however, using any of these two new constructors to create an ill-formed rational (one whose numerator and denominator are not coprime) results in undefined behaviour. I did not document the new constructors for this very reason, but if you believe they should be exported alongside the warning I can add some documentation.

The difference between the two constructors is that the one with Val(true) performs no check whatsoever, whereas the one with Val(false) still checks whether the denominator is negative and errors if it is equal to typemin(T) (to conform to #32569).

EDIT: see discussion below. Intermediate version with unsafe_rational at #35492 (comment). Current version with keyword arguments at #35492 (comment). Other propositions: #35492 (comment) and #35492 (comment)

EDIT: See current state of the PR at #35492 (comment). This PR is non-breaking and does not add any new export.

The new constructor is used to remove some unnecessary checks in the existing code of rational.jl.

The idea for the implementation comes from @JeffreySarnoff's FastRationals.jl. But this PR does not provide any additional speed from relaxing the constraint that the numerator and denominator should be coprime.

The silly benchmark

EDIT: For more serious benchmarking work, see #35492 (comment)

I am using the exact same function as in #35444. This shows the new improvement compared to current master (on which #35444 has been merged).

julia> function foo(::Type{T}) where T<:Integer
           x = Rational{T}(1, 100)
           n = T(10000)
           k = one(T)
           return maximum(i-((i+((x+i)-i))-x) for i in k:n)
       end
foo (generic function with 3 methods)

julia> @btime foo(Int);
  1.147 ms (0 allocations: 0 bytes)   # Currently
  442.850 μs (0 allocations: 0 bytes) # This PR

julia> @btime foo(UInt);
  1.113 ms (0 allocations: 0 bytes)   # Currently
  311.702 μs (0 allocations: 0 bytes) # This PR

julia> @btime foo(BigInt);
  20.428 ms (800012 allocations: 15.56 MiB) # Currently
  12.896 ms (540012 allocations: 10.38 MiB) # This PR

The bonus

This PR also incidentally fixes (what I believe to be) two bugs:

• rationalize(T, x) with T<:Unsigned and x<0 used to yield one(T)//zero(T) aka infinity. It now errors with OverflowError: cannot negate unsigned number.
• typemin(Rational{T}) with T<:Unsigned used to be one(T)//zero(T) (infinity again), it is now zero(T)//one(T).

[KristofferC]

I don't think those actually cause a failure (just scary print out)

[Liozou]

Update on the constructors: the new one and only inner constructor for Rational is unsafe_rational now, which optionally takes a type parameter like the original constructor. It is documented (but not exported, that would require a triage call), here is the documentation:

unsafe_rational(num::Integer, den::Integer; checksign=false)
unsafe_rational{T<:Integer}(num::Integer, den::Integer; checksign::Bool=false)

Create a Rational with numerator num and denominator den.

The unsafe prefix on this function indicates that no check is performed on
num and den to ensure that the resulting Rational is well-formed.
A Rational is well-formed if its numerator and denominator are coprime and if
the denominator is non-negative.

checksign optionally specifies whether to check the sign of den.

Ill-formed Rationals result in undefined behaviour.

The original inner constructor is now an outer constructor with the exact same specification as before. This PR is still entirely non-breaking (apart from the two bugfixes mentionned in the first message).

[Liozou]

Forgot to mention but the idea for unsafe_rational comes from @simeonschaub and thanks a lot to him for the thorough review!

[Liozou]

OK so the only CI error looks unrelated (it passes on 73c5a82 and I only fixed a docstring), and the PR is much better now with the unsafe_rational version. It only provides performance gain (and bugfixes) without breaking or exporting anything new so I guess it should be mergeable, but please let me know if you want me to change or fix something.

[JeffreySarnoff]

To facilitate evaluating this change, (using rational.jl) I created a module Fractions.jl that is this modification with the exported type Fraction rather than Rational and the exported infix op ⨸ (\odiv) rather than \\. The module is given next.

[JeffreySarnoff]

The module is here: https://gist.github.com/JeffreySarnoff/4e13e1cb182908364694a4315b2d10ba

[JeffreySarnoff]

Note: the module requires Julia 1.5.0-DEV.

It appears that this modification provides a useful speedup for Rational multiplication/division (I get 1.15x .. 1.7x, the speedup for Rational{Int64} is more than the speedup for Rational{BigInt}, which is at the lower half of the range). It appears that this modification does not provide a useful speedup for Rational addition/subtraction.

This is somewhat surprising. This probably happens because addition uses the full constructor while multiplication uses the unsafe constructor.

[JeffreySarnoff]

@Liozou try restarting the build

[musm]

The general PR here looks good.

However, we don't generally use lowercase unsafe_rational for types and their constructors.

Was it ever discussed to have leave it implemented via optional keyword argument as in Rational(....; safe=true)?

[JeffreySarnoff]

However, we don't generally use lowercase unsafe_rational for types and their constructors.

Was it ever discussed to have leave it implemented via optional keyword argument as in Rational(....; safe=true)?

Using UnsafeRational in place of unsafe_rational would be cleaner (it is a struct with methods). That way we do not add second keyword into the mix.

[simeonschaub]

I expect using UnsafeRational in place of unsafe_rational would be cleaner (it is a struct with methods).

IMHO, that makes it more confusing, because it's constructor actually creates a Rational. It is basically a hack to specify type parameters to a function.

[JeffreySarnoff]

would just using type parameters be more appropriate?

[simeonschaub]

How? Functions can't accept type parameters.

[JeffreySarnoff]

unsafe_rational is a struct with methods, and it already has a parameter
struct unsafe_rational{T<:Integer} <: Function end
So is Rational.

julia> struct S{A,B} end

julia> S{A,B}() where {A,B} = (A,B)
julia> S{Int,false}()
(Int64, false)

julia> S{A,B}(x) where {A,B} = (A,B,x)
julia> S{true,false}(7)
(true, false, 7)

[simeonschaub]

Yes, but its constructor's only purpose is to create a Rational object. You can still call unsafe_rational{Int}(), but the resulting object isn't of any use.

[JeffreySarnoff]

changing

struct Rational{T<:Integer} <: Real
    num::T
    den::T

    function (::Type{unsafe_rational{T}})(num::Integer, den::Integer; checksign::Bool=false) where T<:Integer
        if checksign && T<:Signed && signbit(den)
            den = -den
            signbit(den) && __throw_rational_argerror_typemin(T)
            num = -num
        end
        return new{T}(num, den)
    end
end

to this seems reasonable

struct Rational{T<:Integer} <: Real
    num::T
    den::T

    function (::Type{unsafe_rational{T}})(num::Integer, den::Integer; checksign::Bool=false) where T<:Signed
        if checksign && signbit(den)
            den = -den
            signbit(den) && __throw_rational_argerror_typemin(T)
            num = -num
        end
        return new{T}(num, den)
    end

    function (::Type{unsafe_rational{T}})(num::Integer, den::Integer; checksign::Bool=false) where T<:Unsigned
        return new{T}(num, den)
    end
end

[sostock]

What about using a third argument instead of the type parameter for unsafe_rational, i.e., unsafe_rational(Int, 7, 2) instead of unsafe_rational{Int}(7, 2)? Then you don’t have to define it as a struct. And one could still have two-argument methods

unsafe_rational(x::Integer, y::Integer) = unsafe_rational(promote(x, y)...)
unsafe_rational(x::T, y::T) where {T<:Integer} = unsafe_rational(T, x, y)

[Liozou]

I added a commit that switches to unsafe_rational as discussed, and rebased on top of the latest changes to rational.jl. The only CI failure looks unrelated.
Allowing unreduced rationals is beyond the point of this PR, so that will be for future work. But nothing in this PR should prevent it anymore.

[Just a note on the implementation detail suggested above: that sign bit does not always exist, typically for Rational{UInt}.]

[StefanKarpinski]

that sign bit does not always exist

I wonder how widespread the use of unsigned rationals actually is in the wild...

[Liozou]

I don't know, but at least one PR was created specifically for Rational{UInt}: #29561.
Anything more I can do for this PR?

[Liozou]

Here are some updated benchmarks for the latest version of the PR, using the benchmarking protocol of @JeffreySarnoff at #35492 (comment). The updated files are RelativeSpeed.jl (no update needed), Fractions.jl and BenchmarkingRationals.jl. The PR is tested against master (commit f1d10e7).

The results are similar to previous except for matrix operations where they are better, and for inv, .+ and .- where they are significantly better:

relspeed(x->x[1] + x[2], rationals2, fractions2)    # 1.0
relspeed(x->x[1] - x[2], rationals2, fractions2)    # 1.0
relspeed(x->x[1] * x[2], rationals2, fractions2)    # 1.6
relspeed(x->x[1] / x[2], rationals2, fractions2)    # 1.6
relspeed(inv, rationals1, fractions1)               # 16772 (!)

relspeed(x -> 2 .+ x, rationals256, fractions256)   # 9.8
relspeed(x -> x .- 1, rationals256, fractions256)   # 9.4
relspeed(x -> 3 .* x, rationals256, fractions256)   # 2.0
relspeed(x -> x ./ 4, rationals256, fractions256)   # 2.7

relspeed(inv, rationals16x16, fractions16x16)       # 1.6
relspeed(inv, rationals4x4, fractions4x4)           # 1.3
relspeed(inv, rationals2x2, fractions2x2)           # 1.3

relspeed(x->inv.(x), rationals4, fractions4)        # 2.7
relspeed(x->inv.(x), rationals16, fractions16)      # 6.3
relspeed(x->inv.(x), rationals64, fractions64)      # 9.8
relspeed(x->inv.(x), rationals256, fractions256)    # 8.2

relspeed(prod, rationals4, fractions4)                                           # 1.6
relspeed(prod, rationals16, fractions16)                                         # 1.7
relspeed(prod, Rational{Int128}.(rationals64), Fraction{Int128}.(fractions64))   # 1.4
relspeed(prod, Rational{BigInt}.(rationals256), Fraction{BigInt}.(fractions256)) # 1.3

relspeed(dot, rationals4, rationals4, fractions4, fractions4)                 # 1.2
relspeed(dot, Rational{BigInt}.(rationals16), Rational{BigInt}.(rationals16),
              Fraction{BigInt}.(fractions16), Fraction{BigInt}.(fractions16)) # 1.1

relspeed(/, rationals16x16, rationals16x16, fractions16x16, fractions16x16)   # 1.3
relspeed(/, rationals4x4, rationals4x4, fractions4x4, fractions4x4)           # 1.3
relspeed(\, rationals16x16, rationals16x16, fractions16x16, fractions16x16)   # 1.3
relspeed(\, rationals4x4, rationals4x4, fractions4x4, fractions4x4)           # 1.2
relspeed(*, rationals16x16, rationals16x16, fractions16x16, fractions16x16)   # 1.1
relspeed(*, rationals4x4, rationals4x4, fractions4x4, fractions4x4)           # 1.1
relspeed(+, rationals16x16, rationals16x16, fractions16x16, fractions16x16)   # 1.0
relspeed(+, rationals4x4, rationals4x4, fractions4x4, fractions4x4)           # 1.0

The difference for the matrix operations comes from an error in the previous benchmark since fractions__x__ were actually defined as matrices of Rational (instead of Fractions). I didn't check the cause for the other improvements, it might also have been a benchmarking issue.

Anyway, it's all above 1 and sometimes significantly so, that's good! Of course I would love for this PR to be merged in time for 1.5, but let me know if there is anything else I can do.

[StefanKarpinski]

Am I correct that as it stands, this PR doesn't change the public API of Rationals at all? It seems like it is, at this point, purely a performance improvement, which is great to have.

I'm not entirely sure about the name unsafe_rational: we usually use the unsafe_ prefix for functions that can lead to Julia crashing somehow, which doesn't seem to be the case here—worst case, you get a wrong answer. Perhaps raw_rational or something like that? Or _Rational?

[sostock]

I'm not entirely sure about the name unsafe_rational

What about reduced_rational or rational_noreduce?

[Liozou]

@StefanKarpinski That's correct, there is no change to the public API of Rationals whatsoever.
@vtjnash had the same concern in #35492 (comment) that unsafe was usually for functions that could make julia crash, but @KristofferC pointed out in #35492 (comment) that there already is unsafe_trunc which is not documented to possibly crash, but only to return potentially garbage value.

[Liozou]

There is also unsafe_length which is not exported and, from what the definitions look like, cannot crash.

[StefanKarpinski]

I've merged. We can discuss the internal name separately. Thanks for the speedup!

[Liozou]

You're welcome, thank you for merging!
It was a long thread, thanks as well for all of you who participated 😀