Merge pull request #31 from JeffBezanson/teh/reductions

Fix reductions so they don't overflow

Author: timholy

src/FixedPointNumbers.jl

@@ -4,13 +4,15 @@ module FixedPointNumbers
 
 using Compat
 
+using Base: IdFun, AddFun, MulFun, reducedim_initarray
+
 import Base: ==, <, <=, -, +, *, /, ~,
              convert, promote_rule, show, showcompact, isinteger, abs, decompose,
              isnan, isinf, isfinite,
              zero, one, typemin, typemax, realmin, realmax, eps, sizeof, reinterpret,
              trunc, round, floor, ceil, bswap,
              div, fld, rem, mod, mod1, rem1, fld1, min, max,
-             start, next, done
+             start, next, done, r_promote, reducedim_init
 # T => BaseType
 # f => Number of Bytes reserved for fractional part
 abstract FixedPoint{T <: Integer, f} <: Real
@@ -60,6 +62,19 @@ include("ufixed.jl")
 include("deprecations.jl")
 
 
+# Promotions for reductions
+const Treduce = Float64
+for F in (AddFun, MulFun)
+    @eval r_promote{T}(::$F, x::FixedPoint{T}) = Treduce(x)
+end
+
+reducedim_init{T<:FixedPoint}(f::IdFun, op::AddFun,
+                              A::AbstractArray{T}, region) =
+    reducedim_initarray(A, region, zero(Treduce))
+reducedim_init{T<:FixedPoint}(f::IdFun, op::MulFun,
+                              A::AbstractArray{T}, region) =
+    reducedim_initarray(A, region, one(Treduce))
+
 # TODO: rewrite this by @generated
 for T in tuple(Fixed16, UF...)
     R = rawtype(T)

test/fixed.jl

@@ -53,3 +53,16 @@ for (TI, f) in [(Int8, 8), (Int16, 8), (Int16, 10), (Int32, 16)]
     println("  Testing $T")
     test_fixed(T, f)
 end
+
+# reductions
+F8 = Fixed{Int8,8}
+a = F8[0.498, 0.1]
+acmp = Float64(a[1]) + Float64(a[2])
+@test sum(a) == acmp
+@test sum(a, 1) == [acmp]
+
+F6 = Fixed{Int8,6}
+a = F6[1.2, 1.4]
+acmp = Float64(a[1])*Float64(a[2])
+@test prod(a) == acmp
+@test prod(a, 1) == [acmp]

test/ufixed.jl

@@ -153,3 +153,13 @@ b, ad = scaledual(0.5, ad)
 generic_scale!(rfloat, a, 0.5)
 generic_scale!(rfixed, ad, b)
 @test rfloat == rfixed
+
+# reductions
+a = UFixed8[0xffuf8, 0xffuf8]
+@test sum(a) == 2.0
+@test sum(a, 1) == [2.0]
+
+a = UFixed14[3.2, 2.4]
+acmp = Float64(a[1])*Float64(a[2])
+@test prod(a) == acmp
+@test prod(a, 1) == [acmp]