Add checked, wrapping and saturating arithmetic for rem/mod (#230)

The default arithmetic for `rem` is still checked arithmetic.
This also adds the support for 3-arg `rem`.
In `RoundUp` mode, `rem` causes `OverflowError` with unsigned (i.e. `Normed`) values.
Otherwise, `rem`/`mod` does not cause an overflow.
This also adds dummy checked, wrapping and saturating rem for type modulus 
to prevent `MethodError` due to lexical replacement of `x % X`.

Author: kimikage

src/FixedPointNumbers.jl

@@ -13,7 +13,7 @@ import Statistics   # for _mean_promote
 import Random: Random, AbstractRNG, SamplerType, rand!
 
 import Base.Checked: checked_neg, checked_abs, checked_add, checked_sub, checked_mul,
-                     checked_div, checked_fld, checked_cld
+                     checked_div, checked_fld, checked_cld, checked_rem, checked_mod
 
 using Base: @pure
 
@@ -38,10 +38,10 @@ export
 # Functions
     scaledual,
     wrapping_neg, wrapping_abs, wrapping_add, wrapping_sub, wrapping_mul,
-    wrapping_fdiv, wrapping_div, wrapping_fld, wrapping_cld,
+    wrapping_div, wrapping_fld, wrapping_cld, wrapping_rem, wrapping_mod,
     saturating_neg, saturating_abs, saturating_add, saturating_sub, saturating_mul,
-    saturating_fdiv, saturating_div, saturating_fld, saturating_cld,
-    checked_fdiv
+    saturating_div, saturating_fld, saturating_cld, saturating_rem, saturating_mod,
+    wrapping_fdiv, saturating_fdiv, checked_fdiv
 
 include("utilities.jl")
 
@@ -64,6 +64,12 @@ wrapping_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = x % X
 saturating_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = clamp(x, X)
 checked_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = _convert(X, x)
 
+# type modulus
+rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
+wrapping_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
+saturating_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
+checked_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
+
 # constructor-style conversions
 (::Type{X})(x::X) where {X <: FixedPoint}      = x
 (::Type{X})(x::Number) where {X <: FixedPoint} = _convert(X, x)
@@ -225,6 +231,9 @@ function wrapping_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <:
 end
 wrapping_fld(x::X, y::X) where {X <: FixedPoint} = wrapping_div(x, y, RoundDown)
 wrapping_cld(x::X, y::X) where {X <: FixedPoint} = wrapping_div(x, y, RoundUp)
+wrapping_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}} =
+    X(x.i - wrapping_div(x, y, r) * y.i, 0)
+wrapping_mod(x::X, y::X) where {X <: FixedPoint} = wrapping_rem(x, y, RoundDown)
 
 # saturating arithmetic
 saturating_neg(x::X) where {X <: FixedPoint} = X(~min(x.i - true, x.i), 0)
@@ -257,6 +266,11 @@ function saturating_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <
 end
 saturating_fld(x::X, y::X) where {X <: FixedPoint} = saturating_div(x, y, RoundDown)
 saturating_cld(x::X, y::X) where {X <: FixedPoint} = saturating_div(x, y, RoundUp)
+function saturating_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
+    T <: Unsigned && r isa RoundingMode{:Up} && return zero(X)
+    X(x.i - saturating_div(x, y, r) * y.i, 0)
+end
+saturating_mod(x::X, y::X) where {X <: FixedPoint} = saturating_rem(x, y, RoundDown)
 
 # checked arithmetic
 checked_neg(x::X) where {X <: FixedPoint} = checked_sub(zero(X), x)
@@ -301,6 +315,16 @@ function checked_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: F
 end
 checked_fld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundDown)
 checked_cld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundUp)
+function checked_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
+    y === zero(X) && throw(DivideError())
+    fx, fy = floattype(X)(x.i), floattype(X)(y.i)
+    z = fx - round(fx / fy, r) * fy
+    if T <: Unsigned && r isa RoundingMode{:Up}
+        z >= zero(z) || throw_overflowerror_rem(r, x, y)
+    end
+    X(_unsafe_trunc(T, z), 0)
+end
+checked_mod(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)
 
 # default arithmetic
 const DEFAULT_ARITHMETIC = :wrapping
@@ -322,6 +346,10 @@ end
 div(x::X, y::X, r::RoundingMode = RoundToZero) where {X <: FixedPoint} = checked_div(x, y, r)
 fld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundDown)
 cld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundUp)
+rem(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundToZero)
+rem(x::X, y::X, ::RoundingMode{:Down}) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)
+rem(x::X, y::X, ::RoundingMode{:Up}) where {X <: FixedPoint} = checked_rem(x, y, RoundUp)
+mod(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)
 
 function minmax(x::X, y::X) where {X <: FixedPoint}
     a, b = minmax(reinterpret(x), reinterpret(y))
@@ -377,7 +405,7 @@ for f in (:~, )
         $f(x::X) where {X <: FixedPoint} = X($f(x.i), 0)
     end
 end
-for f in (:rem, :mod, :mod1, :min, :max)
+for f in (:mod1, :min, :max)
     @eval begin
         $f(x::X, y::X) where {X <: FixedPoint} = X($f(x.i, y.i), 0)
     end
@@ -544,6 +572,11 @@ end
     print(io, op, x, ", ", y, ") overflowed for type ", rawtype(x))
     throw(OverflowError(String(take!(io))))
 end
+@noinline function throw_overflowerror_rem(r::RoundingMode, @nospecialize(x), @nospecialize(y))
+    io = IOBuffer()
+    print(io, "rem(", x, ", ", y, ", ", r, ") overflowed for type ", typeof(x))
+    throw(OverflowError(String(take!(io))))
+end
 
 function Random.rand(r::AbstractRNG, ::SamplerType{X}) where X <: FixedPoint
     X(rand(r, rawtype(X)), 0)

src/fixed.jl

@@ -95,14 +95,14 @@ function _convert(::Type{F}, x::Rational) where {T, f, F <: Fixed{T,f}}
     end
 end
 
-rem(x::F, ::Type{F}) where {F <: Fixed} = x
-function rem(x::Fixed, ::Type{F}) where {T, f, F <: Fixed{T,f}}
+_rem(x::F, ::Type{F}) where {F <: Fixed} = x
+function _rem(x::Fixed, ::Type{F}) where {T, f, F <: Fixed{T,f}}
     f2 = nbitsfrac(typeof(x))
     y = round(@exp2(f - f2) * reinterpret(x))
     reinterpret(F, _unsafe_trunc(T, y))
 end
-rem(x::Integer, ::Type{F}) where {T, f, F <: Fixed{T,f}} = F(_unsafe_trunc(T, x) << f, 0)
-function rem(x::Real, ::Type{F}) where {T, f, F <: Fixed{T,f}}
+_rem(x::Integer, ::Type{F}) where {T, f, F <: Fixed{T,f}} = F(_unsafe_trunc(T, x) << f, 0)
+function _rem(x::Real, ::Type{F}) where {T, f, F <: Fixed{T,f}}
     if bitwidth(T) < 32
         Ti = T
     else
@@ -113,7 +113,7 @@ function rem(x::Real, ::Type{F}) where {T, f, F <: Fixed{T,f}}
     y = _unsafe_trunc(Ti, round(x * Tf(@exp2(f))))
     reinterpret(F, _unsafe_trunc(T, y))
 end
-function rem(x::BigFloat, ::Type{F}) where {T, f, F <: Fixed{T,f}}
+function _rem(x::BigFloat, ::Type{F}) where {T, f, F <: Fixed{T,f}}
     isfinite(x) || return zero(F)
     reinterpret(F, _unsafe_trunc(T, round(x * @exp2(f))))
 end

src/normed.jl

@@ -108,24 +108,25 @@ function _convert(::Type{N}, x::Rational) where {T, f, N <: Normed{T,f}}
     end
 end
 
-rem(x::N, ::Type{N}) where {N <: Normed} = x
-rem(x::Normed, ::Type{N}) where {T, N <: Normed{T}} = reinterpret(N, _unsafe_trunc(T, round((rawone(N)/rawone(x))*reinterpret(x))))
-function rem(x::Real, ::Type{N}) where {T, N <: Normed{T}}
+_rem(x::N, ::Type{N}) where {N <: Normed} = x
+_rem(x::Normed, ::Type{N}) where {T, N <: Normed{T}} =
+    reinterpret(N, _unsafe_trunc(T, round((rawone(N)/rawone(x))*reinterpret(x))))
+function _rem(x::Real, ::Type{N}) where {T, N <: Normed{T}}
     bitwidth(T) < 32 || isfinite(x) || return zero(N)
     reinterpret(N, _unsafe_trunc(T, round(rawone(N) * x)))
 end
-rem(x::Float16, ::Type{N}) where {N <: Normed} = rem(Float32(x), N)  # avoid overflow
+_rem(x::Float16, ::Type{X}) where {X <: Normed} = _rem(Float32(x), X)  # avoid overflow
 # Float32 and Float64 cannot exactly represent `rawone(N)` with `f` greater than
 # the number of their significand bits, resulting in rounding errors (issue #150).
 # So, we use another strategy for the large `f`s explained in:
 # https://github.com/JuliaMath/FixedPointNumbers.jl/pull/166#issuecomment-574135643
-function rem(x::Float32, ::Type{N}) where {f, N <: Normed{UInt32,f}}
+function _rem(x::Float32, ::Type{N}) where {f, N <: Normed{UInt32,f}}
     isfinite(x) || return zero(N)
     f <= 24 && return reinterpret(N, _unsafe_trunc(UInt32, round(rawone(N) * x)))
     r = _unsafe_trunc(UInt32, round(x * @f32(0x1p24)))
     reinterpret(N, r << UInt8(f - 24) - unsigned(signed(r) >> 0x18))
 end
-function rem(x::Float64, ::Type{N}) where {f, N <: Normed{UInt64,f}}
+function _rem(x::Float64, ::Type{N}) where {f, N <: Normed{UInt64,f}}
     isfinite(x) || return zero(N)
     f <= 53 && return reinterpret(N, _unsafe_trunc(UInt64, round(rawone(N) * x)))
     r = _unsafe_trunc(UInt64, round(x * 0x1p53))

test/common.jl

@@ -157,6 +157,7 @@ function test_rem_type(TX::Type)
     @testset "% $X" for X in target(TX, :i8, :i16; ex = :thin)
         xs = typemin(X):0.1:typemax(X)
         @test all(x -> x % X === X(x), xs)
+        @test wrapping_rem(2, X) === saturating_rem(2, X) === checked_rem(2, X) === 2 % X
     end
 end
 
@@ -271,6 +272,33 @@ function test_div_3arg(TX::Type)
     end
 end
 
+function test_rem(TX::Type)
+    for X in target(TX, :i8; ex = :thin)
+        T = rawtype(X)
+        xys = xypairs(X)
+        frem(x, y) = y === zero(y) ? float(x) : x - float(wrapping_div(x, y)) * y
+        fmod(x, y) = y === zero(y) ? float(x) : x - float(wrapping_fld(x, y)) * y
+        frems(x, y) = y === zero(y) ? float(x) : x - float(saturating_div(x, y)) * y
+        fmods(x, y) = y === zero(y) ? float(x) : x - float(saturating_fld(x, y)) * y
+        @test all(((x, y),) -> wrapping_rem(x, y) === frem(x, y) % X, xys)
+        @test all(((x, y),) -> wrapping_mod(x, y) === fmod(x, y) % X, xys)
+        @test all(((x, y),) -> saturating_rem(x, y) === frems(x, y) % X, xys)
+        @test all(((x, y),) -> saturating_mod(x, y) === fmods(x, y) % X, xys)
+        @test all(((x, y),) -> y === zero(y) ||
+                               wrapping_rem(x, y) === checked_rem(x, y), xys)
+        @test all(((x, y),) -> y === zero(y) ||
+                               wrapping_mod(x, y) === checked_mod(x, y), xys)
+    end
+end
+
+function test_rem_3arg(TX::Type)
+    for X in target(TX; ex = :thin)
+        @test rem(eps(X), typemax(X), RoundToZero) === rem(eps(X), typemax(X))
+        @test rem(eps(X), typemax(X), RoundDown)   === mod(eps(X), typemax(X))
+        @test rem(eps(X), eps(X), RoundUp) === zero(X)
+    end
+end
+
 function test_fld1_mod1(TX::Type)
     for X in target(TX, :i8, :i16; ex = :thin)
         T = rawtype(X)

test/fixed.jl

@@ -448,6 +448,51 @@ end
     test_div_3arg(Fixed)
 end
 
+@testset "rem/mod" begin
+    for F in target(Fixed; ex = :thin)
+        fm, fn, fz, fe = typemax(F), typemin(F), zero(F), eps(F)
+        T = rawtype(F)
+        @test   wrapping_rem(fm, fm) ===   wrapping_mod(fm, fm) === fz
+        @test saturating_rem(fm, fm) === saturating_mod(fm, fm) === fz
+        @test    checked_rem(fm, fm) ===    checked_mod(fm, fm) === fz
+
+        @test   wrapping_rem(fz, fe) ===   wrapping_mod(fz, fe) === fz
+        @test saturating_rem(fz, fe) === saturating_mod(fz, fe) === fz
+        @test    checked_rem(fz, fe) ===    checked_mod(fz, fe) === fz
+
+        @test   wrapping_rem(fm, fe) ===   wrapping_mod(fm, fe) === fz
+        @test saturating_rem(fm, fe) === saturating_mod(fm, fe) === fz
+        @test    checked_rem(fm, fe) ===    checked_mod(fm, fe) === fz
+
+        @test   wrapping_rem(fz, fz) ===   wrapping_mod(fz, fz) === fz
+        @test saturating_rem(fz, fz) === saturating_mod(fz, fz) === fz
+        @test_throws DivideError checked_rem(fz, fz)
+        @test_throws DivideError checked_mod(fz, fz)
+
+        @test   wrapping_rem(fe, fz) ===   wrapping_mod(fe, fz) === fe
+        @test saturating_rem(fe, fz) === saturating_mod(fe, fz) === fe
+        @test_throws DivideError checked_rem(fe, fz)
+        @test_throws DivideError checked_mod(fe, fz)
+
+        @test   wrapping_rem(fn, -fe) === wrapping_mod(fn, -fe) === fz
+        @test saturating_rem(fn, -fe) === saturating_mod(fn, -fe) === -fe
+        @test    checked_rem(fn, -fe) ===  checked_mod(fn, -fe) === fz
+
+        @test wrapping_rem(fe, fm) === saturating_rem(fe, fm) === checked_rem(fe, fm) === fe
+        @test wrapping_mod(fe, fm) === saturating_mod(fe, fm) === checked_mod(fe, fm) === fe
+
+        @test wrapping_rem(fe, fn) === saturating_rem(fe, fn) === checked_rem(fe, fn) === fe
+        @test wrapping_mod(fe, fn) === saturating_mod(fe, fn) === checked_mod(fe, fn) === fn + fe
+
+        @test wrapping_rem(fn, fm) === saturating_rem(fn, fm) === checked_rem(fn, fm) === -fe
+        @test wrapping_mod(fn, fm) === saturating_mod(fn, fm) === checked_mod(fn, fm) === fm - fe
+    end
+    test_rem(Fixed)
+    test_rem_3arg(Fixed)
+
+    @test rem(0.5Q0f7, 0.75Q0f7, RoundUp) === -0.25Q0f7
+end
+
 @testset "fld1/mod1" begin
     test_fld1_mod1(Fixed)
 end

test/normed.jl

@@ -442,8 +442,39 @@ end
 end
 
 @testset "rem/mod" begin
-    @test mod(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == rem(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == 0
-    @test mod(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == rem(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == reinterpret(N0f8, 0x01)
+    for N in target(Normed; ex = :thin)
+        nm, nz, ne = typemax(N), zero(N), eps(N)
+        T = rawtype(N)
+        @test   wrapping_rem(nm, nm) ===   wrapping_mod(nm, nm) === nz
+        @test saturating_rem(nm, nm) === saturating_mod(nm, nm) === nz
+        @test    checked_rem(nm, nm) ===    checked_mod(nm, nm) === nz
+
+        @test   wrapping_rem(nz, ne) ===   wrapping_mod(nz, ne) === nz
+        @test saturating_rem(nz, ne) === saturating_mod(nz, ne) === nz
+        @test    checked_rem(nz, ne) ===    checked_mod(nz, ne) === nz
+
+        @test   wrapping_rem(nm, ne) ===   wrapping_mod(nm, ne) === nz
+        @test saturating_rem(nm, ne) === saturating_mod(nm, ne) === nz
+        @test    checked_rem(nm, ne) ===    checked_mod(nm, ne) === nz
+
+        @test   wrapping_rem(nz, nz) ===   wrapping_mod(nz, nz) === nz
+        @test saturating_rem(nz, nz) === saturating_mod(nz, nz) === nz
+        @test_throws DivideError checked_rem(nz, nz)
+        @test_throws DivideError checked_mod(nz, nz)
+
+        @test   wrapping_rem(ne, nz) ===   wrapping_mod(ne, nz) === ne
+        @test saturating_rem(ne, nz) === saturating_mod(ne, nz) === ne
+        @test_throws DivideError checked_rem(ne, nz)
+        @test_throws DivideError checked_mod(ne, nz)
+
+        @test wrapping_rem(ne, nm) === saturating_rem(ne, nm) === checked_rem(ne, nm) === ne
+        @test wrapping_mod(ne, nm) === saturating_mod(ne, nm) === checked_mod(ne, nm) === ne
+    end
+    test_rem(Normed)
+    test_rem_3arg(Normed)
+
+    @test_throws OverflowError rem(0.5N0f8, 1N0f8, RoundUp)
+    @test saturating_rem(0.5N0f8, 1N0f8, RoundUp) === zero(N0f8)
 end
 
 @testset "fld1/mod1" begin

test/runtests.jl

@@ -1,6 +1,6 @@
 using FixedPointNumbers, Test
 
-if VERSION >= v"1.6.0-DEV.816" # JuliaLang/julia #36962
+if VERSION >= v"1.6.0-DEV.816" # JuliaLang/julia #36962 # FIXME
     @test isempty(detect_ambiguities(FixedPointNumbers))
 else
     @test isempty(detect_ambiguities(FixedPointNumbers, Base, Core))