Backport partially "Optimize multiplication for Normed (JuliaMath#213)"

This adds `wrapping_mul` and `checked_mul` binary operations for `Normed`.

This replaces most of Normed's implementation of multiplication with integer operations.
This improves the speed in many cases and the accuracy in some cases.

Author: kimikage

src/FixedPointNumbers.jl

@@ -180,6 +180,9 @@ floattype(::Type{Base.TwicePrecision{T}}) where T<:Union{Float16,Float32} = wide
 
 float(x::FixedPoint) = convert(floattype(x), x)
 
+wrapping_mul(x::X, y::X) where {X <: FixedPoint} = (float(x) * float(y)) % X
+*(x::X, y::X) where {X <: FixedPoint} = wrapping_mul(x, y)
+
 function minmax(x::X, y::X) where {X <: FixedPoint}
     a, b = minmax(reinterpret(x), reinterpret(y))
     X(a,0), X(b,0)

src/normed.jl

@@ -127,7 +127,6 @@ function rem(x::Float64, ::Type{N}) where {f, N <: Normed{UInt64,f}}
     reinterpret(N, r << UInt8(f - 53) - unsigned(signed(r) >> 0x35))
 end
 
-
 function (::Type{T})(x::Normed) where {T <: AbstractFloat}
     # The following optimization for constant division may cause rounding errors.
     # y = reinterpret(x)*(one(rawtype(x))/convert(T, rawone(x)))
@@ -248,8 +247,37 @@ Base.BigFloat(x::Normed) = reinterpret(x) / BigFloat(rawone(x))
 
 Base.Rational(x::Normed) = reinterpret(x)//rawone(x)
 
-# unchecked arithmetic
-*(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)*convert(floattype(T), y))
+# Division by `2^f-1` with RoundNearest. The result would be in the lower half bits.
+div_2fm1(x::T, ::Val{f}) where {T, f} = (x + (T(1)<<(f - 1) - 0x1)) ÷ (T(1) << f - 0x1)
+div_2fm1(x::T, ::Val{1}) where T = x
+div_2fm1(x::UInt16,  ::Val{8})  = (((x + 0x80) >> 0x8) + x + 0x80) >> 0x8
+div_2fm1(x::UInt32,  ::Val{16}) = (((x + 0x8000) >> 0x10) + x + 0x8000) >> 0x10
+div_2fm1(x::UInt64,  ::Val{32}) = (((x + 0x80000000) >> 0x20) + x + 0x80000000) >> 0x20
+div_2fm1(x::UInt128, ::Val{64}) = (((x + 0x8000000000000000) >> 0x40) + x + 0x8000000000000000) >> 0x40
+
+# wrapping arithmetic
+function wrapping_mul(x::N, y::N) where {T <: Union{UInt8,UInt16,UInt32,UInt64}, f, N <: Normed{T,f}}
+    z = widemul(x.i, y.i)
+    N(div_2fm1(z, Val(Int(f))) % T, 0)
+end
+
+# checked arithmetic
+function checked_mul(x::N, y::N) where {N <: Normed}
+    z = float(x) * float(y)
+    z < typemax(N) + eps(N)/2 || throw_overflowerror(:*, x, y)
+    z % N
+end
+function checked_mul(x::N, y::N) where {T <: Union{UInt8,UInt16,UInt32,UInt64}, f, N <: Normed{T,f}}
+    f == bitwidth(T) && return wrapping_mul(x, y)
+    z = widemul(x.i, y.i)
+    m = widemul(typemax(N).i, rawone(N)) + (rawone(N) >> 0x1)
+    z < m || throw_overflowerror(:*, x, y)
+    N(div_2fm1(z, Val(Int(f))) % T, 0)
+end
+
+# Override the default arithmetic with `checked` for backward compatibility
+*(x::N, y::N) where {N <: Normed} = checked_mul(x, y)
+
 /(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)/convert(floattype(T), y))
 
 # Functions

test/fixed.jl

@@ -256,6 +256,28 @@ end
     @test (-67.2 % T).i == round(Int, -67.2*512) % Int16
 end
 
+@testset "mul" begin
+    wrapping_mul(x::F, y::F) where {F <: Fixed} = x * y
+    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
+
+        @test wrapping_mul(typemin(F), typemin(F)) === big(typemin(F))^2 % F
+    end
+    for F in target(Fixed, :i8; ex = :thin)
+        xs = typemin(F):eps(F):typemax(F)
+        xys = ((x, y) for x in xs, y in xs)
+        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
+end
+
 @testset "rounding" begin
     for sym in (:i8, :i16, :i32, :i64)
         T = symbol_to_inttype(Fixed, sym)

test/normed.jl

@@ -283,6 +283,26 @@ end
     end
 end
 
+@testset "mul" begin
+    checked_mul = FixedPointNumbers.checked_mul
+    for N in target(Normed; ex = :thin)
+        @test checked_mul(typemax(N), zero(N)) === zero(N)
+
+        @test checked_mul(one(N), typemax(N)) === typemax(N)
+
+        if typemax(N) != 1
+            @test_throws OverflowError checked_mul(typemax(N), typemax(N))
+        end
+    end
+    for N in target(Normed, :i8; ex = :thin)
+        xs = typemin(N):eps(N):typemax(N)
+        xys = ((x, y) for x in xs, y in xs)
+        fmul(x, y) = float(x) * float(y) # note that precision(Float32) < 32
+        @test all(((x, y),) -> !(typemin(N) <= fmul(x, y) <= typemax(N)) ||
+                               (fmul(x, y) % N) === checked_mul(x, y), xys)
+    end
+end
+
 @testset "div/fld1" begin
     @test div(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == fld(reinterpret(N0f8, 0x10), reinterpret(N0f8, 0x02)) == 8
     @test div(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == fld(reinterpret(N0f8, 0x0f), reinterpret(N0f8, 0x02)) == 7