Backport partially "Specialize multiplication for Fixed (JuliaMath#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,35 @@ 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
+
+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

@@ -246,6 +246,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)
@@ -259,16 +263,16 @@ 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 "mul" begin
     wrapping_mul = FixedPointNumbers.wrapping_mul
     for F in target(Fixed; ex = :thin)
         @test wrapping_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 wrapping_mul(typemin(F), typemax(F)) === big(typemin(F)) * big(typemax(F)) % F
@@ -281,6 +285,13 @@ end
         fmul(x, y) = float(x) * float(y) # note that precision(Float32) < 32
         @test all(((x, y),) -> wrapping_mul(x, y) === fmul(x, y) % F, 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

@@ -257,8 +257,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)