Specialize multiplication for Fixed (#220)

This specializes most of the multiplication for `Fixed` and avoids floating point operations.
A major change is that the rounding mode is changed from `RoundNearestTiesUp` to `RoundNearest`.
The existing `RoundNearestTiesUp` and `RoundDown` modes are now supported by the new unexported function `mul_with_rounding`.
This also improves `rem`.
Unlike multiplication for `Normed`, the wrapping arithmetic is the default for `Fixed`.

Author: kimikage

src/fixed.jl

@@ -95,19 +95,20 @@ function _convert(::Type{F}, x::Rational) where {T, f, F <: Fixed{T,f}}
     end
 end
 
-# unchecked arithmetic
-
-# with truncation:
-#*(x::Fixed{T,f}, y::Fixed{T,f}) = Fixed{T,f}(Base.widemul(x.i,y.i)>>f,0)
-# with rounding up:
-*(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}((Base.widemul(x.i,y.i) + (one(widen(T)) << (f-1)))>>f,0)
-
-/(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(div(convert(widen(T), x.i) << f, y.i), 0)
-
-
-rem(x::Integer, ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(x,T)<<f,0)
-rem(x::Real,    ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(Integer(trunc(x)),T)<<f + rem(Integer(round(rem(x,1)*(one(widen1(T))<<f))),T),0)
-
+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}}
+    y = _unsafe_trunc(promote_type(Int64, T), round(x * @exp2(f)))
+    reinterpret(F, _unsafe_trunc(T, y))
+end
+function rem(x::BigFloat, ::Type{F}) where {T, f, F <: Fixed{T,f}}
+    reinterpret(F, _unsafe_trunc(T, round(x * @exp2(f))))
+end
 
 (::Type{Tf})(x::Fixed{T,f}) where {Tf <: AbstractFloat, T, f} = Tf(Tf(x.i) * Tf(@exp2(-f)))
 Base.Float16(x::Fixed{T,f}) where {T, f} = Float16(Float32(x))
@@ -118,6 +119,60 @@ function Base.Rational(x::Fixed{T,f}) where {T, f}
     f < bitwidth(T)-1 ? x.i//rawone(x) : x.i//(one(widen1(T))<<f)
 end
 
+function div_2f(x::T, ::Val{f}) where {T, f}
+    xf = x & (T(-1) >>> (bitwidth(T) - f - 1))
+    half = oneunit(T) << (f - 1)
+    c = half - (xf === half)
+    (x + c) >> f
+end
+div_2f(x::T, ::Val{0}) where {T} = x
+
+# wrapping arithmetic
+function wrapping_mul(x::F, y::F) where {T <: Union{Int8,Int16,Int32,Int64}, f, F <: Fixed{T,f}}
+    z = widemul(x.i, y.i)
+    F(div_2f(z, Val(Int(f))) % T, 0)
+end
+
+# saturating arithmetic
+function saturating_mul(x::F, y::F) where {T <: Union{Int8,Int16,Int32,Int64}, f, F <: Fixed{T,f}}
+    if f >= bitwidth(T) - 1
+        z = min(widemul(x.i, y.i), widen1(typemax(T)) << f)
+    else
+        z = clamp(widemul(x.i, y.i), widen1(typemin(T)) << f, widen1(typemax(T)) << f)
+    end
+    F(div_2f(z, Val(Int(f))) % T, 0)
+end
+
+# checked arithmetic
+function checked_mul(x::F, y::F) where {T <: Union{Int8,Int16,Int32,Int64}, f, F <: Fixed{T,f}}
+    z = widemul(x.i, y.i)
+    if f < 1
+        m = widen1(typemax(T)) + 0x1
+        n = widen1(typemin(T))
+    else
+        half = widen1(oneunit(T)) << (f - 1)
+        m = widen1(typemax(T)) << f + half
+        n = widen1(typemin(T)) << f - half
+    end
+    (n <= z) & (z < m) || throw_overflowerror(:*, x, y)
+    F(div_2f(z, Val(Int(f))) % T, 0)
+end
+
+function mul_with_rounding(x::F, y::F, ::RoundingMode{:Nearest}) where {F <: Fixed}
+    wrapping_mul(x, y)
+end
+function mul_with_rounding(x::F, y::F, ::RoundingMode{:NearestTiesUp}) where
+                          {T <: Union{Int8,Int16,Int32,Int64}, f, F <: Fixed{T, f}}
+    z = widemul(x.i, y.i)
+    F(((z + (oftype(z, 1) << f >>> 1)) >> f) % T, 0)
+end
+function mul_with_rounding(x::F, y::F, ::RoundingMode{:Down}) where
+                          {T <: Union{Int8,Int16,Int32,Int64}, f, F <: Fixed{T, f}}
+    F((widemul(x.i, y.i) >> f) % T, 0)
+end
+
+/(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(div(convert(widen(T), x.i) << f, y.i), 0)
+
 function trunc(x::Fixed{T,f}) where {T, f}
     f == 0 && return x
     f == bitwidth(T) && return zero(x) # TODO: remove this line

test/fixed.jl

@@ -263,6 +263,10 @@ end
 end
 
 @testset "type modulus" begin
+    @test  Q0f7(0.2) % Q0f7  === Q0f7(0.2)
+    @test Q1f14(1.2) % Q0f15 === Q0f15(-0.8)
+    @test Q1f14(1.2) % Q0f7  === Q0f7(-0.8)
+
     T = Fixed{Int8,7}
     for i = -1.0:0.1:typemax(T)
         @test i % T === T(i)
@@ -276,6 +280,9 @@ end
     end
     @test ( 65.2 % T).i == round(Int,  65.2*512) % Int16
     @test (-67.2 % T).i == round(Int, -67.2*512) % Int16
+
+    @test -1 % Q0f7 === Q0f7(-1)
+    @test -2 % Q0f7 === Q0f7(0)
 end
 
 @testset "neg" begin
@@ -379,9 +386,6 @@ end
         @test saturating_mul(typemax(F), zero(F)) === zero(F)
         @test    checked_mul(typemax(F), zero(F)) === zero(F)
 
-        # FIXME: Both the rhs and lhs of the following tests may be inaccurate due to `rem`
-        F === Fixed{Int128,127} && continue
-
         @test   wrapping_mul(F(-1), typemax(F)) === -typemax(F)
         @test saturating_mul(F(-1), typemax(F)) === -typemax(F)
         @test    checked_mul(F(-1), typemax(F)) === -typemax(F)
@@ -405,6 +409,13 @@ end
         @test all(((x, y),) -> !(typemin(F) <= fmul(x, y) <= typemax(F)) ||
                                wrapping_mul(x, y) === checked_mul(x, y), xys)
     end
+
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1,  0.5Q6f1, RoundNearest) ===  1.0Q6f1
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1, -0.5Q6f1, RoundNearest) === -1.0Q6f1
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1,  0.5Q6f1, RoundNearestTiesUp) ===  1.0Q6f1
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1, -0.5Q6f1, RoundNearestTiesUp) === -0.5Q6f1
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1,  0.5Q6f1, RoundDown) ===  0.5Q6f1
+    FixedPointNumbers.mul_with_rounding(1.5Q6f1, -0.5Q6f1, RoundDown) === -1.0Q6f1
 end
 
 @testset "rounding" begin

test/normed.jl

@@ -268,8 +268,8 @@ end
     @test (65.2 % N6f10).i == round(Int, 65.2*1023) % UInt16
     @test (-0.3 % N6f10).i == round(Int, -0.3*1023) % UInt16
 
-    @test 1 % N0f8 == 1
-    @test 2 % N0f8 == N0f8(0.996)
+    @test 1 % N0f8 === N0f8(1)
+    @test 2 % N0f8 === N0f8(0.996)
 
     # issue #150
     @test all(f -> 1.0f0 % Normed{UInt32,f} == oneunit(Normed{UInt32,f}), 1:32)