div with rounding modes [+ rounded division] (JuliaLang#33040)

* Re-arrange fld/cld code

In preparation for supporting other rounding modes in div, create
a three-argument div function that takes a rounding mode as the last
argument and make this the fundamental fallback for fld/cld.

* Implemented rounded division

* add various divrem combinations to avoid overflow

* Whitespace/test fixes

* Small tweaks to docstrings

* Bugfixes

* Add the exhaustive test

* Tigthen up types for div fallback

I think it's better to give a MethodError here than an approximate
answer for non-AbstractFloat reals (e.g. custom integer types).

NEWS.md

@@ -36,6 +36,8 @@ Standard library changes
 
 * The methods of `mktemp` and `mktempdir` which take a function body to pass temporary paths to no longer throw errors if the path is already deleted when the function body returns ([#33091]).
 
+* `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])
+
 * Verbose `display` of `Char` (`text/plain` output) now shows the codepoint value in standard-conforming `"U+XXXX"` format ([#33291]).
 
 #### Libdl

base/Base.jl

@@ -116,6 +116,7 @@ include("namedtuple.jl")
 include("hashing.jl")
 include("rounding.jl")
 using .Rounding
+include("div.jl")
 include("float.jl")
 include("twiceprecision.jl")
 include("complex.jl")

base/bool.jl

@@ -112,7 +112,5 @@ end
 *(y::AbstractFloat, x::Bool) = x * y
 
 div(x::Bool, y::Bool) = y ? x : throw(DivideError())
-fld(x::Bool, y::Bool) = div(x,y)
-cld(x::Bool, y::Bool) = div(x,y)
 rem(x::Bool, y::Bool) = y ? false : throw(DivideError())
 mod(x::Bool, y::Bool) = rem(x,y)

base/div.jl

@@ -0,0 +1,258 @@
+# 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=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
+
+ 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`. Equivalent to `div(x, y, RoundDown)`.
+
+See also: [`div`](@ref)
+
+# Examples
+```jldoctest
+julia> fld(7.3,5.5)
+1.0
+```
+"""
+fld(a, b) = div(a, b, RoundDown)
+
+"""
+ cld(x, y)
+
+Smallest integer larger than or equal to `x/y`. Equivalent to `div(x, y, RoundUp)`.
+
+See also: [`div`](@ref)
+
+# Examples
+```jldoctest
+julia> cld(5.5,2.2)
+3.0
+```
+"""
+cld(a, b) = div(a, b, RoundUp)
+
+# divrem
+"""
+ divrem(x, y, r::RoundingMode=RoundToZero)
+
+The quotient and remainder from Euclidean division.
+Equivalent to `(div(x,y,r), rem(x,y,r))`. Equivalently, with the the default
+value of `r`, this call is equivalent to `(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))
+function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearest))
+ (q, r) = divrem(x, y)
+ if x >= 0
+ if y >= 0
+ r >= (y÷2) + (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
+ else
+ r >= -(y÷2) + (isodd(y) | iseven(q)) ? (q-true, r+y) : (q, r)
+ end
+ else
+ if y >= 0
+ r <= -signed(y÷2) - (isodd(y) | iseven(q)) ? (q-true, r+y) : (q, r)
+ else
+ r <= (y÷2) - (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
+ end
+ end
+end
+function divrem(x::Integer, y::Integer, rnd:: typeof(RoundNearestTiesAway))
+ (q, r) = divrem(x, y)
+ if x >= 0
+ if y >= 0
+ r >= (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
+ else
+ r >= -(y÷2) + isodd(y) ? (q-true, r+y) : (q, r)
+ end
+ else
+ if y >= 0
+ r <= -signed(y÷2) - isodd(y) ? (q-true, r+y) : (q, r)
+ else
+ r <= (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
+ end
+ end
+end
+function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearestTiesUp))
+ (q, r) = divrem(x, y)
+ if x >= 0
+ if y >= 0
+ r >= (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
+ else
+ r >= -(y÷2) + true ? (q-true, r+y) : (q, r)
+ end
+ else
+ if y >= 0
+ r <= -signed(y÷2) - true ? (q-true, r+y) : (q, r)
+ else
+ r <= (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
+ end
+ end
+end
+
+"""
+ 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.
+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
+
+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)})
+ divrem(x,y,rnd)[1]
+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)
+# Prevent ambiguities
+for rm in (RoundUp, RoundDown, RoundToZero)
+ @eval div(a::Bool, b::Bool, r::$(typeof(rm))) = div(a, b)
+end
+function div(x::Bool, y::Bool, rnd::Union{typeof(RoundNearest),
+ typeof(RoundNearestTiesAway),
+ typeof(RoundNearestTiesUp)})
+ div(x, y)
+end
+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)
+div(a::Int128, b::Int128, ::typeof(RoundToZero)) = div(a, b)
+div(a::UInt128, b::UInt128, ::typeof(RoundToZero)) = div(a, b)
+rem(a::Int128, b::Int128, ::typeof(RoundToZero)) = rem(a, b)
+rem(a::UInt128, b::UInt128, ::typeof(RoundToZero)) = rem(a, b)
+
+# 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)
+
+# Integers
+# fld(x,y) == div(x,y) - ((x>=0) != (y>=0) && rem(x,y) != 0 ? 1 : 0)
+div(x::T, y::T, ::typeof(RoundDown)) where {T<:Unsigned} = div(x,y)
+function div(x::T, y::T, ::typeof(RoundDown)) where T<:Integer
+ d = div(x, y, RoundToZero)
+ return d - (signbit(x ⊻ y) & (d * y != x))
+end
+
+# cld(x,y) = div(x,y) + ((x>0) == (y>0) && rem(x,y) != 0 ? 1 : 0)
+function div(x::T, y::T, ::typeof(RoundUp)) where T<:Unsigned
+ d = div(x, y, RoundToZero)
+ return d + (d * y != x)
+end
+function div(x::T, y::T, ::typeof(RoundUp)) where T<:Integer
+ d = div(x, y, RoundToZero)
+ return d + (((x > 0) == (y > 0)) & (d * y != x))
+end
+
+# Real
+# NOTE: C89 fmod() and x87 FPREM implicitly provide truncating float division,
+# so it is used here as the basis of float div().
+div(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = convert(T,round((x-rem(x,y,r))/y))
+rem(x::T, y::T, ::typeof(RoundUp)) where {T<:AbstractFloat} = 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
@@ -465,14 +466,25 @@ big(n::Integer) = convert(BigInt, n)
 
 # Binary ops
 for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
- (:fld, :fdiv_q), (:div, :tdiv_q), (:mod, :fdiv_r), (:rem, :tdiv_r),
+ (:mod, :fdiv_r), (:rem, :tdiv_r),
  (:gcd, :gcd), (:lcm, :lcm),
  (:&, :and), (:|, :ior), (:xor, :xor))
  @eval begin
  ($fJ)(x::BigInt, y::BigInt) = MPZ.$fC(x, y)
  end
 end
 
+for (r, f) in ((RoundToZero, :tdiv_q),
+ (RoundDown, :fdiv_q),
+ (RoundUp, :cdiv_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/int.jl

@@ -175,8 +175,15 @@ div(x::Unsigned, y::BitSigned) = unsigned(flipsign(signed(div(x, unsigned(abs(y)
 rem(x::BitSigned, y::Unsigned) = flipsign(signed(rem(unsigned(abs(x)), y)), x)
 rem(x::Unsigned, y::BitSigned) = rem(x, unsigned(abs(y)))
 
-fld(x::Signed, y::Unsigned) = div(x, y) - (signbit(x) & (rem(x, y) != 0))
-fld(x::Unsigned, y::Signed) = div(x, y) - (signbit(y) & (rem(x, y) != 0))
+function divrem(x::BitSigned, y::Unsigned)
+ q, r = divrem(unsigned(abs(x)), y)
+ flipsign(signed(q), x), flipsign(signed(r), x)
+end
+
+function divrem(x::Unsigned, y::BitSigned)
+ q, r = divrem(x, unsigned(abs(y)))
+ unsigned(flipsign(signed(q), y)), r
+end
 
 
 """
@@ -220,34 +227,13 @@ mod(x::BitSigned, y::Unsigned) = rem(y + unsigned(rem(x, y)), y)
 mod(x::Unsigned, y::Signed) = rem(y + signed(rem(x, y)), y)
 mod(x::T, y::T) where {T<:Unsigned} = rem(x, y)
 
-cld(x::Signed, y::Unsigned) = div(x, y) + (!signbit(x) & (rem(x, y) != 0))
-cld(x::Unsigned, y::Signed) = div(x, y) + (!signbit(y) & (rem(x, y) != 0))
-
 # Don't promote integers for div/rem/mod since there is no danger of overflow,
 # while there is a substantial performance penalty to 64-bit promotion.
 div(x::T, y::T) where {T<:BitSigned64} = checked_sdiv_int(x, y)
 rem(x::T, y::T) where {T<:BitSigned64} = checked_srem_int(x, y)
 div(x::T, y::T) where {T<:BitUnsigned64} = checked_udiv_int(x, y)
 rem(x::T, y::T) where {T<:BitUnsigned64} = checked_urem_int(x, y)
 
-
-# fld(x,y) == div(x,y) - ((x>=0) != (y>=0) && rem(x,y) != 0 ? 1 : 0)
-fld(x::T, y::T) where {T<:Unsigned} = div(x,y)
-function fld(x::T, y::T) where T<:Integer
- d = div(x, y)
- return d - (signbit(x ⊻ y) & (d * y != x))
-end
-
-# cld(x,y) = div(x,y) + ((x>0) == (y>0) && rem(x,y) != 0 ? 1 : 0)
-function cld(x::T, y::T) where T<:Unsigned
- d = div(x, y)
- return d + (d * y != x)
-end
-function cld(x::T, y::T) where T<:Integer
- d = div(x, y)
- return d + (((x > 0) == (y > 0)) & (d * y != x))
-end
-
 ## integer bitwise operations ##
 
 """