Merge pull request #158 from kimikage/round

timholy · web-flow · commit 402e637ce784 · 2020-01-11T15:06:50.000-06:00
Add rounding functions for `Fixed` (Fixes #153)
diff --git a/src/FixedPointNumbers.jl b/src/FixedPointNumbers.jl
@@ -139,6 +139,15 @@ for f in (:+, :-, :rem, :mod, :mod1, :min, :max)
         $f(x::X, y::X) where {X <: FixedPoint} = X($f(x.i, y.i), 0)
     end
 end
+for (m, f) in ((:(:Nearest), :round),
+               (:(:ToZero), :trunc),
+               (:(:Up), :ceil),
+               (:(:Down), :floor))
+    @eval begin
+        round(x::FixedPoint, ::RoundingMode{$m}) = $f(x)
+        round(::Type{Ti}, x::FixedPoint, ::RoundingMode{$m}) where {Ti <: Integer} = $f(Ti, x)
+    end
+end
 
 # Printing. These are used to generate type-symbols, so we need them
 # before we include any files.
diff --git a/src/fixed.jl b/src/fixed.jl
@@ -50,6 +50,9 @@ function rawone(::Type{Fixed{T,f}}) where {T, f}
     oneunit(T) << f
 end
 
+intmask(::Fixed{T,f}) where {T, f} = -oneunit(T) << f # Signed
+fracmask(x::Fixed{T,f}) where {T, f} = ~intmask(x) # Signed
+
 # unchecked arithmetic
 
 # with truncation:
@@ -91,6 +94,73 @@ end
 (::Type{TR})(x::Fixed{T,f}) where {TR <: Rational,T,f} =
     TR(x.i>>f + (x.i&(1<<f-1))//(one(widen1(T))<<f))
 
+function trunc(x::Fixed{T,f}) where {T, f}
+    f == 0 && return x
+    f == bitwidth(T) && return zero(x) # TODO: remove this line
+    f == bitwidth(T) - 1 && return x.i == typemin(T) ? x : zero(x)
+    t = x.i & intmask(x)
+    r = x.i & fracmask(x)
+    _rawone = oneunit(T) << f
+    reinterpret(Fixed{T,f}, (x.i < 0) & (r != 0) ? t + _rawone : t)
+end
+function floor(x::Fixed{T,f}) where {T, f}
+    f == bitwidth(T) && x.i < 0 && throw_converterror(Fixed{T,f}, -1) # TODO: remove this line
+    Fixed{T,f}(x.i & intmask(x), 0)
+end
+function ceil(x::Fixed{T,f}) where {T, f}
+    f == 0 && return x
+    upper = typemax(T) & intmask(x)
+    x.i > upper && throw_converterror(Fixed{T,f}, ceil(float(x)))
+    reinterpret(Fixed{T,f}, (x.i + fracmask(x)) & intmask(x))
+end
+function round(x::Fixed{T,f}) where {T, f}
+    f == 0 && return x
+    f == bitwidth(T) && return zero(x) # TODO: remove this line
+    upper = intmask(x) >>> 0x1
+    lower = intmask(x) >> 0x1
+    if f == bitwidth(T) - 1
+        x.i > upper && throw_converterror(Fixed{T,f}, @exp2(bitwidth(T)-f-1))
+        return x.i < lower ? typemin(x) : zero(x)
+    end
+    x.i >= upper && throw_converterror(Fixed{T,f}, @exp2(bitwidth(T)-f-1))
+    y = oneunit(T) << UInt8(f - 1) + x.i
+    m = oneunit(T) << UInt8(f + 1) - oneunit(T)
+    z = y & intmask(x)
+    reinterpret(Fixed{T,f}, z - T(y & m == rawone(x)) << f)
+end
+
+function trunc(::Type{Ti}, x::Fixed{T,f}) where {Ti <: Integer, T, f}
+    f == 0 && return convert(Ti, x.i)
+    f == bitwidth(T) && return zero(Ti) # TODO: remove this line
+    f == bitwidth(T) - 1 && return x.i == typemin(T) ? convert(Ti, -1) : zero(Ti)
+    t = x.i >> f
+    r = x.i & fracmask(x)
+    convert(Ti, (x.i < 0) & (r != 0) ? t + oneunit(T) : t)
+end
+function floor(::Type{Ti}, x::Fixed{T,f}) where {Ti <: Integer, T, f}
+    f == bitwidth(T) && return x.i < 0 ? convert(Ti, -1) : zero(Ti) # TODO: remove this line
+    convert(Ti, x.i >> f)
+end
+function ceil(::Type{Ti}, x::Fixed{T,f}) where {Ti <: Integer, T, f}
+    f == bitwidth(T) && return x.i > 0 ? oneunit(Ti) : zero(Ti) # TODO: remove this line
+    y = x.i + fracmask(x)
+    convert(Ti, x.i >= 0 ? y >>> f : y >> f)
+end
+function round(::Type{Ti}, x::Fixed{T,f}) where {Ti <: Integer, T, f}
+    f == 0 && return convert(Ti, x.i)
+    f == bitwidth(T) && return zero(Ti) # TODO: remove this line
+    upper = intmask(x) >>> 0x1
+    lower = intmask(x) >> 0x1
+    if f == bitwidth(T) - 1
+        x.i < lower && return convert(Ti, -1)
+        return x.i > upper ? oneunit(Ti) : zero(Ti)
+    end
+    y = oneunit(T) << UInt8(f - 1) + x.i
+    m = oneunit(T) << UInt8(f + 1) - oneunit(T)
+    z = x.i >= 0 ? y >>> f : y >> f
+    convert(Ti, z - Ti(y & m == rawone(x)))
+end
+
 promote_rule(ft::Type{Fixed{T,f}}, ::Type{TI}) where {T,f,TI <: Integer} = Fixed{T,f}
 promote_rule(::Type{Fixed{T,f}}, ::Type{TF}) where {T,f,TF <: AbstractFloat} = TF
 promote_rule(::Type{Fixed{T,f}}, ::Type{Rational{TR}}) where {T,f,TR} = Rational{TR}
diff --git a/src/normed.jl b/src/normed.jl
@@ -238,26 +238,35 @@ abs(x::Normed) = x
 /(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)/convert(floattype(T), y))
 
 # Functions
-trunc(x::T) where {T <: Normed} = T(div(reinterpret(x), rawone(T))*rawone(T),0)
-floor(x::T) where {T <: Normed} = trunc(x)
-function round(x::Normed{T,f}) where {T,f}
-    mask = convert(T, 1<<(f-1))
-    y = trunc(x)
-    return convert(T, reinterpret(x)-reinterpret(y)) & mask>0 ?
-            Normed{T,f}(y+oneunit(Normed{T,f})) : y
+trunc(x::N) where {N <: Normed} = floor(x)
+floor(x::N) where {N <: Normed} = reinterpret(N, x.i - x.i % rawone(N))
+function ceil(x::Normed{T,f}) where {T, f}
+    f == 1 && return x
+    if typemax(T) % rawone(x) != 0
+        upper = typemax(T) - typemax(T) % rawone(x)
+        x.i > upper && throw_converterror(Normed{T,f}, ceil(T, typemax(x)))
+    end
+    r = x.i % rawone(x)
+    reinterpret(Normed{T,f}, x.i - r + (r > 0 ? rawone(x) : zero(T)))
 end
-function ceil(x::Normed{T,f}) where {T,f}
-    k = bitwidth(T)-f
-    mask = (typemax(T)<<k)>>k
-    y = trunc(x)
-    return convert(T, reinterpret(x)-reinterpret(y)) & (mask)>0 ?
-            Normed{T,f}(y+oneunit(Normed{T,f})) : y
+function round(x::Normed{T,f}) where {T, f}
+    r = x.i % rawone(x)
+    q = rawone(x) - r
+    reinterpret(Normed{T,f}, r > q ? x.i + q : x.i - r)
 end
 
-trunc(::Type{T}, x::Normed) where {T <: Integer} = convert(T, div(reinterpret(x), rawone(x)))
-round(::Type{T}, x::Normed) where {T <: Integer} = round(T, reinterpret(x)/rawone(x))
-floor(::Type{T}, x::Normed) where {T <: Integer} = trunc(T, x)
- ceil(::Type{T}, x::Normed) where {T <: Integer} =  ceil(T, reinterpret(x)/rawone(x))
+trunc(::Type{Ti}, x::Normed) where {Ti <: Integer} = floor(Ti, x)
+function floor(::Type{Ti}, x::Normed) where {Ti <: Integer}
+    convert(Ti, reinterpret(x) ÷ rawone(x))
+end
+function ceil(::Type{Ti}, x::Normed) where {Ti <: Integer}
+    d, r = divrem(x.i, rawone(x))
+    convert(Ti, r > 0 ? d + oneunit(rawtype(x)) : d)
+end
+function round(::Type{Ti}, x::Normed) where {Ti <: Integer}
+    d, r = divrem(x.i, rawone(x))
+    convert(Ti, r > (rawone(x) >> 0x1) ? d + oneunit(rawtype(x)) : d)
+end
 
 isfinite(x::Normed) = true
 isnan(x::Normed) = false
diff --git a/test/fixed.jl b/test/fixed.jl
@@ -70,6 +70,18 @@ end
     @test reinterpret(Int8, 0.5Q0f7) === signed(0x40)
 end
 
+@testset "masks" begin
+    @test FixedPointNumbers.intmask(0Q7f0) === signed(0xFF)
+    @test FixedPointNumbers.intmask(0Q6f1) === signed(0xFE)
+    @test FixedPointNumbers.intmask(0Q1f6) === signed(0xC0)
+    @test FixedPointNumbers.intmask(0Q0f7) === signed(0x80)
+
+    @test FixedPointNumbers.fracmask(0Q7f0) === signed(0x00)
+    @test FixedPointNumbers.fracmask(0Q6f1) === signed(0x01)
+    @test FixedPointNumbers.fracmask(0Q1f6) === signed(0x3F)
+    @test FixedPointNumbers.fracmask(0Q0f7) === signed(0x7F)
+end
+
 @testset "inexactness" begin
     @test_throws InexactError Q0f7(-2)
     # TODO: change back to InexactError when it allows message strings
@@ -96,6 +108,56 @@ end
     end
 end
 
+@testset "rounding" begin
+    for T in (Int8, Int16, Int32, Int64)
+        rs = vcat([ oneunit(T) << b - oneunit(T) for b = 0:bitwidth(T)-1],
+                  [ oneunit(T) << b              for b = 1:bitwidth(T)-2],
+                  [ oneunit(T) << b + oneunit(T) for b = 2:bitwidth(T)-2],
+                  [-oneunit(T) << b - oneunit(T) for b = 2:bitwidth(T)-2],
+                  [-oneunit(T) << b              for b = 1:bitwidth(T)-1],
+                  [-oneunit(T) << b + oneunit(T) for b = 1:bitwidth(T)-1])
+        @testset "rounding Fixed{$T,$f}" for f = 0:bitwidth(T)-1
+            F = Fixed{T,f}
+            xs = (reinterpret(F, r) for r in rs)
+            @test all(x -> trunc(x) == trunc(float(x)), xs)
+            @test all(x -> floor(float(x)) < typemin(F) || floor(x) == floor(float(x)), xs)
+            @test all(x ->  ceil(float(x)) > typemax(F) ||  ceil(x) ==  ceil(float(x)), xs)
+            @test all(x -> round(float(x)) > typemax(F) || round(x) == round(float(x)), xs)
+            @test all(x -> trunc(Int64, x) === trunc(Int64, float(x)), xs)
+            @test all(x -> floor(Int64, x) === floor(Int64, float(x)), xs)
+            @test all(x ->  ceil(Int64, x) ===  ceil(Int64, float(x)), xs)
+            @test all(x -> round(Int64, x) === round(Int64, float(x)), xs)
+        end
+    end
+    @testset "rounding Fixed{Int16,$f} with overflow" for f = 1:16 # TODO: drop 16
+        F = Fixed{Int16,f}
+        @test_throws ArgumentError ceil(typemax(F))
+        if f == 16
+            @test_throws ArgumentError ceil(eps(F))
+        elseif f == 15
+            @test_throws ArgumentError ceil(eps(F))
+            @test_throws ArgumentError round(typemax(F))
+            @test_throws ArgumentError round(F(0.5) + eps(F))
+        else
+            @test_throws ArgumentError ceil(typemin(F) - oneunit(F) + eps(F))
+            @test_throws ArgumentError round(typemax(F))
+            @test_throws ArgumentError round(typemax(F) - F(0.5) + eps(F))
+        end
+    end
+    @testset "rounding mode" begin
+        @test round(-1.5Q1f6, RoundNearest) === -2Q1f6
+        @test round(-1.5Q1f6, RoundToZero) === -1Q1f6
+        @test round(-1.5Q1f6, RoundUp) === -1Q1f6
+        @test round(-1.5Q1f6, RoundDown) === -2Q1f6
+        @test round(Int, -1.5Q1f6, RoundNearest) === -2
+        @test round(Int, -1.5Q1f6, RoundToZero) === -1
+        @test round(Int, -1.5Q1f6, RoundUp) === -1
+        @test round(Int, -1.5Q1f6, RoundDown) === -2
+    end
+    @test_throws InexactError trunc(UInt, typemin(Q0f7))
+    @test_throws InexactError floor(UInt, -eps(Q0f7))
+end
+
 @testset "modulus" begin
     T = Fixed{Int8,7}
     for i = -1.0:0.1:typemax(T)
diff --git a/test/normed.jl b/test/normed.jl
@@ -240,41 +240,39 @@ end
     end
 end
 
-function testtrunc(inc::T) where {T}
-    incf = convert(Float64, inc)
-    tm = reinterpret(typemax(T))/reinterpret(one(T))
-    local x = zero(T)
-    for i = 0 : min(1e6, reinterpret(typemax(T))-1)
-        xf = incf*i
-        try
-            @test typeof(trunc(x)) == T
-            @test trunc(x) == trunc(xf)
-            @test typeof(round(x)) == T
-            @test round(x) == round(xf)
-            cxf = ceil(xf)
-            if cxf < tm
-                @test typeof(ceil(x)) == T
-                @test ceil(x) == ceil(xf)
-            end
-            @test typeof(floor(x)) == T
-            @test floor(x) == floor(xf)
-            @test trunc(Int,x) == trunc(Int,xf)
-            @test round(Int,x) == round(Int,xf)
-            @test floor(Int,x) == floor(Int,xf)
-            if cxf < tm
-                @test ceil(Int,x) == ceil(Int,xf)
-            end
-        catch err
-            println("Failed on x = ", x, ", xf = ", xf)
-            rethrow(err)
+@testset "rounding" begin
+    for T in (UInt8, UInt16, UInt32, UInt64)
+        rs = vcat([ oneunit(T) << b - oneunit(T) << 1 for b = 1:bitwidth(T)],
+                  [ oneunit(T) << b - oneunit(T)      for b = 1:bitwidth(T)],
+                  [ oneunit(T) << b                   for b = 2:bitwidth(T)-1])
+        @testset "rounding Normed{$T,$f}" for f = 1:bitwidth(T)
+            N = Normed{T,f}
+            xs = (reinterpret(N, r) for r in rs)
+            @test all(x -> trunc(x) == trunc(float(x)), xs)
+            @test all(x -> floor(x) == floor(float(x)), xs)
+            # force `Normed` comparison avoiding rounding errors
+            @test all(x -> ceil(float(x)) > typemax(N) || ceil(x) == N(ceil(float(x))), xs)
+            @test all(x -> round(x) == round(float(x)), xs)
+            @test all(x -> trunc(UInt64, x) === trunc(UInt64, float(x)), xs)
+            @test all(x -> floor(UInt64, x) === floor(UInt64, float(x)), xs)
+            @test all(x ->  ceil(UInt64, x) ===  ceil(UInt64, float(x)), xs)
+            @test all(x -> round(UInt64, x) === round(UInt64, float(x)), xs)
         end
-        x = convert(T, x+inc)
     end
-end
-
-@testset "trunc" begin
-    for T in (FixedPointNumbers.UF..., UF2...)
-        testtrunc(eps(T))
+    @testset "rounding Normed{Int16,$f} with overflow" for f in filter(x->!ispow2(x), 1:16)
+        N = Normed{UInt16,f}
+        @test_throws ArgumentError ceil(typemax(N))
+        @test_throws ArgumentError ceil(floor(typemax(N)) + eps(N))
+    end
+    @testset "rounding mode" begin
+        @test round(1.504N1f7, RoundNearest) === 2N1f7
+        @test round(1.504N1f7, RoundToZero) === 1N1f7
+        @test round(1.504N1f7, RoundUp) === 2N1f7
+        @test round(1.504N1f7, RoundDown) === 1N1f7
+        @test round(Int, 1.504N1f7, RoundNearest) === 2
+        @test round(Int, 1.504N1f7, RoundToZero) === 1
+        @test round(Int, 1.504N1f7, RoundUp) === 2
+        @test round(Int, 1.504N1f7, RoundDown) === 1
     end
 end
 
