Implemented rounded division

NEWS.md

@@ -26,6 +26,8 @@ New library functions
 Standard library changes
 ------------------------
 
+* `div` now accepts a rounding mode as the third argument, consistent with the corresponding argument to `rem`. Support for rounding division, by passing one of the RoundNearest modes to this function, was added. For future compatibility, library authors should now extend this function, rather than extending the two-argument `fld`/`cld`/`div` directly. ([#33040])
+
 
 #### Libdl
 

base/div.jl

@@ -1,10 +1,77 @@
 # Div is truncating by default
+
+"""
+    div(x, y, r::RoundingMode=RoundToZero)
+
+Compute the remainder of `x` after integer division by `y`, with the quotient rounded
+according to the rounding mode `r`. In other words, the quantity
+
+    y*round(x/y,r)
+
+without any intermediate rounding.
+
+See also: [`fld`](@ref), [`cld`](@ref) which are special cases of this function
+
+# Examples:
+```jldoctest
+julia> div(4, 3, RoundDown) # Matches fld(4, 3)
+1
+julia> div(4, 3, RoundUp) # Matches cld(4, 3)
+2
+julia> div(5, 2, RoundNearest)
+2
+julia> div(5, 2, RoundNearestTiesAway)
+3
+julia> div(-5, 2, RoundNearest)
+-2
+julia> div(-5, 2, RoundNearestTiesAway)
+-3
+julia> div(-5, 2, RoundNearestTiesUp)
+-2
+```
+"""
+div(x, y, r::RoundingMode)
+
 div(a, b) = div(a, b, RoundToZero)
 
+"""
+    rem(x, y, r::RoundingMode)
+
+Compute the remainder of `x` after integer division by `y`, with the quotient rounded
+according to the rounding mode `r`. In other words, the quantity
+
+    x - y*round(x/y,r)
+
+without any intermediate rounding.
+
+- if `r == RoundNearest`, then the result is exact, and in the interval
+  ``[-|y|/2, |y|/2]``. See also [`RoundNearest`](@ref).
+
+- if `r == RoundToZero` (default), then the result is exact, and in the interval
+  ``[0, |y|)`` if `x` is positive, or ``(-|y|, 0]`` otherwise. See also [`RoundToZero`](@ref).
+
+- if `r == RoundDown`, then the result is in the interval ``[0, y)`` if `y` is positive, or
+  ``(y, 0]`` otherwise. The result may not be exact if `x` and `y` have different signs, and
+  `abs(x) < abs(y)`. See also[`RoundDown`](@ref).
+
+- if `r == RoundUp`, then the result is in the interval `(-y,0]` if `y` is positive, or
+  `[0,-y)` otherwise. The result may not be exact if `x` and `y` have the same sign, and
+  `abs(x) < abs(y)`. See also [`RoundUp`](@ref).
+
+"""
+rem(x, y, r::RoundingMode)
+
+# TODO: Make these primitive and have the two-argument version call these
+rem(x, y, ::RoundingMode{:ToZero}) = rem(x,y)
+rem(x, y, ::RoundingMode{:Down}) = mod(x,y)
+rem(x, y, ::RoundingMode{:Up}) = mod(x,-y)
+
 """
     fld(x, y)
 
-Largest integer less than or equal to `x/y`.
+Largest integer less than or equal to `x/y`. Equivalent to `div(x, y, RoundDown)`.
+
+See also: [`div`](@ref)
 
 # Examples
 ```jldoctest
@@ -17,7 +84,9 @@ fld(a, b) = div(a, b, RoundDown)
 """
     cld(x, y)
 
-Smallest integer larger than or equal to `x/y`.
+Smallest integer larger than or equal to `x/y`. Equivalent to `div(x, y, RoundUp)`.
+
+See also: [`div`](@ref)
 
 # Examples
 ```jldoctest
@@ -27,22 +96,101 @@ julia> cld(5.5,2.2)
 """
 cld(a, b) = div(a, b, RoundUp)
 
+# divrem
+"""
+    divrem(x, y)
+
+The quotient and remainder from Euclidean division. Equivalent to `(div(x,y), rem(x,y))` or
+`(x÷y, x%y)`.
+
+# Examples
+```jldoctest
+julia> divrem(3,7)
+(0, 3)
+
+julia> divrem(7,3)
+(2, 1)
+```
+"""
+divrem(x, y) = divrem(x, y, RoundToZero)
+divrem(a, b, r::RoundingMode) = (div(a, b, r), rem(a, b, r))
+
+"""
+    fldmod(x, y)
+
+The floored quotient and modulus after division. A convenience wrapper for
+`divrem(x, y, RoundDown)`. Equivalent to `(fld(x,y), mod(x,y))`.
+"""
+fldmod(x,y) = divrem(x, y, RoundDown)
+
 # We definite generic rounding methods for other rounding modes in terms of
 # RoundToZero.
-div(x::Signed, y::Unsigned, ::typeof(RoundDown)) = div(x, y, RoundToZero) - (signbit(x) & (rem(x, y) != 0))
-div(x::Unsigned, y::Signed, ::typeof(RoundDown)) = div(x, y, RoundToZero) - (signbit(y) & (rem(x, y) != 0))
+function div(x::Signed, y::Unsigned, ::typeof(RoundDown))
+    (q, r) = divrem(x, y)
+    q - (signbit(x) & (r != 0))
+end
+function div(x::Unsigned, y::Signed, ::typeof(RoundDown))
+    (q, r) = divrem(x, y)
+    q - (signbit(y) & (r != 0))
+end
 
-div(x::Signed, y::Unsigned, ::typeof(RoundUp)) = div(x, y, RoundToZero) + (!signbit(x) & (rem(x, y) != 0))
-div(x::Unsigned, y::Signed, ::typeof(RoundUp)) = div(x, y, RoundToZero) + (!signbit(y) & (rem(x, y) != 0))
+function div(x::Signed, y::Unsigned, ::typeof(RoundUp))
+    (q, r) = divrem(x, y)
+    q + (!signbit(x) & (r != 0))
+end
+function div(x::Unsigned, y::Signed, ::typeof(RoundUp))
+    (q, r) = divrem(x, y)
+    q + (!signbit(y) & (r != 0))
+end
+
+function div(x::Integer, y::Integer, rnd::Union{typeof(RoundNearest),
+                                              typeof(RoundNearestTiesAway),
+                                              typeof(RoundNearestTiesUp)})
+    (q, r) = divrem(x, y)
+    # No remainder, we're done
+    iszero(r) && return q
+    if isodd(y)
+        # The divisior is odd - no ties are possible, just
+        # round to nearest. N.B. 2r == y is impossible because
+        # y is odd.
+        return abs(2r) > abs(y) ? q + copysign(one(q), q) : q
+    end
+    # y is even, divide y by two and check whether we have a tie.
+    # If not, as above we just round to the nearest value.
+    halfy = copysign(y >> 1, r)
+    c = cmp(r, halfy)
+    c == -1 && return q
+    c == 1 && return q + copysign(one(q), q)
+    # We have a tie (r == y/2). Select tie behavior according to
+    # rounding mode.
+    if rnd == RoundNearest
+        return iseven(q) ? q : q + copysign(one(q), q)
+    elseif rnd == RoundNearestTiesAway
+        return q + copysign(one(q), q)
+    else
+        @assert rnd == RoundNearestTiesUp
+        return sign(q) == -1 ? q : q + one(q)
+    end
+end
 
 # For bootstrapping purposes, we define div for integers directly. Provide the
 # generic signature also
 div(a::T, b::T, ::typeof(RoundToZero)) where {T<:Union{BitSigned, BitUnsigned64}} = div(a, b)
 div(a::Bool, b::Bool, r::RoundingMode) = div(a, b)
+fld(a::T, b::T) where {T<:Union{Integer,AbstractFloat}} = div(a, b, RoundDown)
+cld(a::T, b::T) where {T<:Union{Integer,AbstractFloat}} = div(a, b, RoundUp)
 
-# For compatibility
-fld(a::T, b::T) where {T<:Integer} = div(a, b, RoundDown)
-cld(a::T, b::T) where {T<:Integer} = div(a, b, RoundDown)
+# These are kept for compatibility with external packages overriding fld/cld.
+# In 2.0, packages should extend div(a,b,r) instead, in which case, these can
+# be removed.
+fld(x::Real, y::Real) = div(promote(x,y)..., RoundDown)
+cld(x::Real, y::Real) = div(promote(x,y)..., RoundUp)
+fld(x::Signed, y::Unsigned) = div(x, y, RoundDown)
+fld(x::Unsigned, y::Signed) = div(x, y, RoundDown)
+cld(x::Signed, y::Unsigned) = div(x, y, RoundUp)
+cld(x::Unsigned, y::Signed) = div(x, y, RoundUp)
+fld(x::T, y::T) where {T<:Real} = throw(MethodError(div, (x, y, RoundDown)))
+cld(x::T, y::T) where {T<:Real} = throw(MethodError(div, (x, y, RoundUp)))
 
 # Promotion
 div(x::Real, y::Real, r::RoundingMode) = div(promote(x, y)..., r)
@@ -66,9 +214,5 @@ function div(x::T, y::T, ::typeof(RoundUp)) where T<:Integer
 end
 
 # Real
-div(x::T, y::T, ::typeof(RoundDown)) where {T<:Real} = convert(T,round((x-mod(x,y))/y))
-
-div(x::T, y::T, ::typeof(RoundUp)) where {T<:Real} = convert(T,round((x-modCeil(x,y))/y))
-#rem(x::T, y::T) where {T<:Real} = convert(T,x-y*trunc(x/y))
-#mod(x::T, y::T) where {T<:Real} = convert(T,x-y*floor(x/y))
-modCeil(x::T, y::T) where {T<:Real} = convert(T,x-y*ceil(x/y))
+div(x::T, y::T, r::RoundingMode) where {T<:Real} = convert(T,round((x-rem(x,y,r))/y))
+rem(x::T, y::T, ::typeof(RoundUp)) where {T<:Real} = convert(T,x-y*ceil(x/y))

base/gmp.jl

@@ -5,7 +5,7 @@ module GMP
 export BigInt
 
 import .Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), xor,
-             binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
+             binomial, cmp, convert, div, divrem, factorial, cld, fld, gcd, gcdx, lcm, mod,
              ndigits, promote_rule, rem, show, isqrt, string, powermod,
              sum, trailing_zeros, trailing_ones, count_ones, tryparse_internal,
              bin, oct, dec, hex, isequal, invmod, _prevpow2, _nextpow2, ndigits0zpb,
@@ -146,7 +146,8 @@ sizeinbase(a::BigInt, b) = Int(ccall((:__gmpz_sizeinbase, :libgmp), Csize_t, (mp
 
 for (op, nbits) in (:add => :(BITS_PER_LIMB*(1 + max(abs(a.size), abs(b.size)))),
                     :sub => :(BITS_PER_LIMB*(1 + max(abs(a.size), abs(b.size)))),
-                    :mul => 0, :fdiv_q => 0, :tdiv_q => 0, :fdiv_r => 0, :tdiv_r => 0,
+                    :mul => 0, :fdiv_q => 0, :tdiv_q => 0, :cdiv_q => 0,
+                    :fdiv_r => 0, :tdiv_r => 0, :cdiv_r => 0,
                     :gcd => 0, :lcm => 0, :and => 0, :ior => 0, :xor => 0)
     op! = Symbol(op, :!)
     @eval begin
@@ -479,6 +480,11 @@ for (r, f) in ((RoundToZero, :tdiv_q),
     @eval div(x::BigInt, y::BigInt, ::typeof($r)) = MPZ.$f(x, y)
 end
 
+# For compat only. Remove in 2.0.
+div(x::BigInt, y::BigInt) = div(x, y, RoundToZero)
+fld(x::BigInt, y::BigInt) = div(x, y, RoundDown)
+cld(x::BigInt, y::BigInt) = div(x, y, RoundUp)
+
 /(x::BigInt, y::BigInt) = float(x)/float(y)
 
 function invmod(x::BigInt, y::BigInt)

base/math.jl

@@ -779,35 +779,6 @@ function frexp(x::T) where T<:IEEEFloat
     return reinterpret(T, xu), k
 end
 
-"""
-    rem(x, y, r::RoundingMode)
-
-Compute the remainder of `x` after integer division by `y`, with the quotient rounded
-according to the rounding mode `r`. In other words, the quantity
-
-    x - y*round(x/y,r)
-
-without any intermediate rounding.
-
-- if `r == RoundNearest`, then the result is exact, and in the interval
-  ``[-|y|/2, |y|/2]``. See also [`RoundNearest`](@ref).
-
-- if `r == RoundToZero` (default), then the result is exact, and in the interval
-  ``[0, |y|)`` if `x` is positive, or ``(-|y|, 0]`` otherwise. See also [`RoundToZero`](@ref).
-
-- if `r == RoundDown`, then the result is in the interval ``[0, y)`` if `y` is positive, or
-  ``(y, 0]`` otherwise. The result may not be exact if `x` and `y` have different signs, and
-  `abs(x) < abs(y)`. See also[`RoundDown`](@ref).
-
-- if `r == RoundUp`, then the result is in the interval `(-y,0]` if `y` is positive, or
-  `[0,-y)` otherwise. The result may not be exact if `x` and `y` have the same sign, and
-  `abs(x) < abs(y)`. See also [`RoundUp`](@ref).
-
-"""
-rem(x, y, ::RoundingMode{:ToZero}) = rem(x,y)
-rem(x, y, ::RoundingMode{:Down}) = mod(x,y)
-rem(x, y, ::RoundingMode{:Up}) = mod(x,-y)
-
 rem(x::Float64, y::Float64, ::RoundingMode{:Nearest}) =
     ccall((:remainder, libm),Float64,(Float64,Float64),x,y)
 rem(x::Float32, y::Float32, ::RoundingMode{:Nearest}) =

base/number.jl

@@ -87,30 +87,6 @@ first(x::Number) = x
 last(x::Number) = x
 copy(x::Number) = x # some code treats numbers as collection-like
 
-"""
-    divrem(x, y)
-
-The quotient and remainder from Euclidean division. Equivalent to `(div(x,y), rem(x,y))` or
-`(x÷y, x%y)`.
-
-# Examples
-```jldoctest
-julia> divrem(3,7)
-(0, 3)
-
-julia> divrem(7,3)
-(2, 1)
-```
-"""
-divrem(x,y) = (div(x,y),rem(x,y))
-
-"""
-    fldmod(x, y)
-
-The floored quotient and modulus after division. Equivalent to `(fld(x,y), mod(x,y))`.
-"""
-fldmod(x,y) = (fld(x,y),mod(x,y))
-
 """
     signbit(x)
 

base/promotion.jl

@@ -351,12 +351,6 @@ div(x::Real, y::Real) = div(promote(x,y)...)
 rem(x::Real, y::Real) = rem(promote(x,y)...)
 mod(x::Real, y::Real) = mod(promote(x,y)...)
 
-# These are kept for compatibility with external packages overriding fld/cld.
-# In 2.0, packages should extend div(a,b,r) instead, in which case, these can
-# be removed.
-fld(x::Real, y::Real) = fld(promote(x,y)...)
-cld(x::Real, y::Real) = cld(promote(x,y)...)
-
 mod1(x::Real, y::Real) = mod1(promote(x,y)...)
 fld1(x::Real, y::Real) = fld1(promote(x,y)...)
 

base/rational.jl

@@ -364,6 +364,18 @@ function div(x::Rational, y::Rational, r::RoundingMode)
     div(checked_mul(xn,yd), checked_mul(xd,yn), r)
 end
 
+# For compatibility - to be removed in 2.0 when the generic fallbacks
+# are removed from div.jl
+div(x::T, y::T, r::RoundingMode) where {T<:Rational} =
+    invoke(div, Tuple{Rational, Rational, RoundingMode}, x, y, r)
+for (S, T) in ((Rational, Integer), (Integer, Rational), (Rational, Rational))
+    @eval begin
+        div(x::$S, y::$T) = div(x, y, RoundToZero)
+        fld(x::$S, y::$T) = div(x, y, RoundDown)
+        cld(x::$S, y::$T) = div(x, y, RoundUp)
+    end
+end
+
 trunc(::Type{T}, x::Rational) where {T} = convert(T,div(x.num,x.den))
 floor(::Type{T}, x::Rational) where {T} = convert(T,fld(x.num,x.den))
 ceil(::Type{T}, x::Rational) where {T} = convert(T,cld(x.num,x.den))

test/int.jl

@@ -305,3 +305,26 @@ let i = MyInt26779(1)
     @test_throws MethodError i << 1
     @test_throws MethodError i >>> 1
 end
+
+@testset "rounding division" begin
+    for (a, b, nearest, away, up) in (
+            (3, 2, 2, 2, 2),
+            (5, 3, 2, 2, 2),
+            (-3, 2, -2, -2, -1),
+            (5, 2, 2, 3, 3),
+            (-5, 2, -2, -3, -2),
+            (-5, 3, -2, -2, -2),
+            (5, -3, -2, -2, -2))
+        for sign in (+1, -1)
+            (a, b) = (a*sign, b*sign)
+            @test div(a, b, RoundNearest) == nearest
+            @test div(a, b, RoundNearestTiesAway) == away
+            @test div(a, b, RoundNearestTiesUp) == up
+        end
+    end
+
+    @test div(typemax(Int64), typemax(Int64)-1, RoundNearest) == 1
+    @test div(-typemax(Int64), typemax(Int64)-1, RoundNearest) == -2
+    @test div(typemax(Int64), 2, RoundNearest) == 4611686018427387904
+    @test div(-typemax(Int64), 2, RoundNearestTiesUp) == -4611686018427387903
+end