Merge branch 'cv/limits' (#8)

Author: omus

src/FixedPointDecimals.jl

@@ -29,10 +29,10 @@ export FixedDecimal, RoundThrows
 
 using Compat
 
-import Base: reinterpret, zero, one, abs, sign, ==, <, <=, +, -, /, *, div,
-             rem, divrem, fld, mod, fldmod, fld1, mod1, fldmod1, isinteger,
-             typemin, typemax, realmin, realmax, print, show, string, convert, parse,
-             promote_rule, min, max, trunc, round, floor, ceil, eps, float, widemul
+import Base: reinterpret, zero, one, abs, sign, ==, <, <=, +, -, /, *, div, rem, divrem,
+             fld, mod, fldmod, fld1, mod1, fldmod1, isinteger, typemin, typemax,
+             realmin, realmax, print, show, string, convert, parse, promote_rule, min, max,
+             trunc, round, floor, ceil, eps, float, widemul
 
 const IEEEFloat = Union{Float16, Float32, Float64}
 
@@ -70,17 +70,26 @@ for fn in [:trunc, :floor, :ceil]
 end
 
 """
-    FixedDecimal{I <: Integer, f::Int}
+    FixedDecimal{T <: Integer, f::Int}
 
-A fixed-point decimal type backed by integral type `I`, with `f` digits after
+A fixed-point decimal type backed by integral type `T`, with `f` digits after
 the decimal point stored.
 """
 immutable FixedDecimal{T <: Integer, f} <: Real
     i::T
 
-    # internal constructor
-    Base.reinterpret{T, f}(::Type{FixedDecimal{T, f}}, i::Integer) =
-        new{T, f}(i % T)
+    # inner constructor
+    function Base.reinterpret{T, f}(::Type{FixedDecimal{T, f}}, i::Integer)
+        n = max_exp10(T)
+        if f >= 0 && (n < 0 || f <= n)
+            new{T, f}(i % T)
+        else
+            throw(ArgumentError(
+                "Requested number of decimal places $f exceeds the max allowed for the " *
+                "storage type $T: [0, $(max_exp10(T))]"
+            ))
+        end
+    end
 end
 
 const FD = FixedDecimal
@@ -131,8 +140,8 @@ _round_to_even(q, r, d) = _round_to_even(promote(q, r, d)...)
 # TODO: can we use floating point to speed this up? after we build a
 # correctness test suite.
 function *{T, f}(x::FD{T, f}, y::FD{T, f})
-    powt = T(10)^f
-    quotient, remainder = fldmod(Base.widemul(x.i, y.i), powt)
+    powt = coefficient(FD{T, f})
+    quotient, remainder = fldmod(widemul(x.i, y.i), powt)
     reinterpret(FD{T, f}, _round_to_even(quotient, remainder, powt))
 end
 
@@ -142,14 +151,15 @@ end
 *{T, f}(x::FD{T, f}, y::Integer) = reinterpret(FD{T, f}, T(x.i * y))
 
 function /{T, f}(x::FD{T, f}, y::FD{T, f})
-    powt = T(10)^f
+    powt = coefficient(FD{T, f})
     quotient, remainder = fldmod(widemul(x.i, powt), y.i)
     reinterpret(FD{T, f}, T(_round_to_even(quotient, remainder, y.i)))
 end
 
-# these functions are needed to avoid InexactError when converting from the integer type
+# These functions allow us to perform division with integers outside of the range of the
+# FixedDecimal.
 function /{T, f}(x::Integer, y::FD{T, f})
-    powt = T(10)^f
+    powt = coefficient(FD{T, f})
     powtsq = widemul(powt, powt)
     quotient, remainder = fldmod(widemul(x, powtsq), y.i)
     reinterpret(FD{T, f}, T(_round_to_even(quotient, remainder, y.i)))
@@ -161,16 +171,17 @@ function /{T, f}(x::FD{T, f}, y::Integer)
 end
 
 # integerification
-trunc{T, f}(x::FD{T, f}) = FD{T, f}(div(x.i, T(10)^f))
-floor{T, f}(x::FD{T, f}) = FD{T, f}(fld(x.i, T(10)^f))
+trunc{T, f}(x::FD{T, f}) = FD{T, f}(div(x.i, coefficient(FD{T, f})))
+floor{T, f}(x::FD{T, f}) = FD{T, f}(fld(x.i, coefficient(FD{T, f})))
+
 # TODO: round with number of digits; should be easy
 function round{T, f}(x::FD{T, f}, ::RoundingMode{:Nearest}=RoundNearest)
-    powt = T(10)^f
+    powt = coefficient(FD{T, f})
     quotient, remainder = fldmod(x.i, powt)
     FD{T, f}(_round_to_even(quotient, remainder, powt))
 end
 function ceil{T, f}(x::FD{T, f})
-    powt = T(10)^f
+    powt = coefficient(FD{T, f})
     quotient, remainder = fldmod(x.i, powt)
     if remainder > 0
         FD{T, f}(quotient + one(quotient))
@@ -185,49 +196,48 @@ for fn in [:trunc, :floor, :ceil]
     # round/trunc/ceil/flooring to FD; generic
     # TODO. we may need to check overflow and boundary conditions here.
     @eval function $fn{T, f}(::Type{FD{T, f}}, x::Real)
-        powt = T(10)^f
+        powt = coefficient(FD{T, f})
         val = trunc(T, $(Symbol(fn, "mul"))(x, powt))
         reinterpret(FD{T, f}, val)
     end
 end
-round{TI <: Integer}(::Type{TI}, x::FD,
-                     ::RoundingMode{:Nearest}=RoundNearest)::TI = round(x)
-function round{T, f}(::Type{FD{T, f}}, x::Real,
-                     ::RoundingMode{:Nearest}=RoundNearest)
-    reinterpret(FD{T, f}, round(T, T(10)^f * x))
+function round{TI <: Integer}(::Type{TI}, x::FD, ::RoundingMode{:Nearest}=RoundNearest)::TI
+    round(x)
+end
+function round{T, f}(::Type{FD{T, f}}, x::Real, ::RoundingMode{:Nearest}=RoundNearest)
+    reinterpret(FD{T, f}, round(T, x * coefficient(FD{T, f})))
 end
 
 # needed to avoid ambiguity
-function round{T, f}(::Type{FD{T, f}}, x::Rational,
-                     ::RoundingMode{:Nearest}=RoundNearest)
-    reinterpret(FD{T, f}, round(T, T(10)^f * x))
+function round{T, f}(::Type{FD{T, f}}, x::Rational, ::RoundingMode{:Nearest}=RoundNearest)
+    reinterpret(FD{T, f}, round(T, x * coefficient(FD{T, f})))
 end
 
 # conversions and promotions
-convert{T, f}(::Type{FD{T, f}}, x::Integer) =
-    reinterpret(FD{T, f}, round(T, Base.widemul(T(x), T(10)^f)))
+function convert{T, f}(::Type{FD{T, f}}, x::Integer)
+    reinterpret(FD{T, f}, T(widemul(x, coefficient(FD{T, f}))))
+end
 
 convert{T <: FD}(::Type{T}, x::AbstractFloat) = round(T, x)
 
 function convert{T, f}(::Type{FD{T, f}}, x::Rational)::FD{T, f}
-    powt = T(10)^f
-    num::T, den::T = numerator(x), denominator(x)
-    g = gcd(powt, den)
-    powt = div(powt, g)
-    den = div(den, g)
-    reinterpret(FD{T, f}, powt * num) / FD{T, f}(den)
+    powt = coefficient(FD{T, f})
+    reinterpret(FD{T, f}, T(x * powt))
 end
 
 function convert{T, f, U, g}(::Type{FD{T, f}}, x::FD{U, g})
     if f ≥ g
-        reinterpret(FD{T, f}, convert(T, Base.widemul(T(10)^(f-g), x.i)))
+        # Compute `10^(f - g)` without overflow
+        powt = div(coefficient(FD{T, f}), coefficient(FD{U, g}))
+        reinterpret(FD{T, f}, T(widemul(x.i, powt)))
     else
-        sf = T(10)^(g - f)
-        q, r = divrem(x.i, sf)
-        if r ≠ 0
-            throw(InexactError())
+        # Compute `10^(g - f)` without overflow
+        powt = div(coefficient(FD{U, g}), coefficient(FD{T, f}))
+        q, r = divrem(x.i, powt)
+        if r == 0
+            reinterpret(FD{T, f}, T(q))
         else
-            reinterpret(FD{T, f}, convert(T, q))
+            throw(InexactError())
         end
     end
 end
@@ -240,16 +250,20 @@ for divfn in [:div, :fld, :fld1]
 end
 
 convert(::Type{AbstractFloat}, x::FD) = convert(floattype(typeof(x)), x)
-convert{TF <: AbstractFloat, T, f}(::Type{TF}, x::FD{T, f})::TF =
-    x.i / TF(10)^f
+function convert{TF <: AbstractFloat, T, f}(::Type{TF}, x::FD{T, f})::TF
+    x.i / coefficient(FD{T, f})
+end
+
+function convert{TF <: BigFloat, T, f}(::Type{TF}, x::FD{T, f})::TF
+    BigInt(x.i) / BigInt(coefficient(FD{T, f}))
+end
 
 function convert{TI <: Integer, T, f}(::Type{TI}, x::FD{T, f})::TI
     isinteger(x) || throw(InexactError())
-    div(x.i, T(10)^f)
+    div(x.i, coefficient(FD{T, f}))
 end
 
-convert{TR<:Rational,T,f}(::Type{TR}, x::FD{T, f})::TR =
-    x.i // T(10)^f
+convert{TR <: Rational, T, f}(::Type{TR}, x::FD{T, f})::TR = x.i // coefficient(FD{T, f})
 
 promote_rule{T, f, TI <: Integer}(::Type{FD{T, f}}, ::Type{TI}) = FD{T, f}
 promote_rule{T, f, TF <: AbstractFloat}(::Type{FD{T, f}}, ::Type{TF}) = TF
@@ -266,7 +280,7 @@ Base.@pure promote_rule{T, f, U, g}(::Type{FD{T, f}}, ::Type{FD{U, g}}) =
 <={T <: FD}(x::T, y::T) = x.i <= y.i
 
 # predicates and traits
-isinteger{T, f}(x::FD{T, f}) = rem(x.i, T(10)^f) == 0
+isinteger{T, f}(x::FD{T, f}) = rem(x.i, coefficient(FD{T, f})) == 0
 typemin{T, f}(::Type{FD{T, f}}) = reinterpret(FD{T, f}, typemin(T))
 typemax{T, f}(::Type{FD{T, f}}) = reinterpret(FD{T, f}, typemax(T))
 eps{T <: FD}(::Type{T}) = reinterpret(T, 1)
@@ -283,7 +297,7 @@ function print{T, f}(io::IO, x::FD{T, f})
 
     # note: a is negative if x.i == typemin(x.i)
     s, a = sign(x.i), abs(x.i)
-    integer, fractional = divrem(a, T(10)^f)
+    integer, fractional = divrem(a, coefficient(x))
     integer = abs(integer)  # ...but since f > 0, this is positive
     fractional = abs(fractional)
 
@@ -382,4 +396,38 @@ function parse_round{T}(::Type{T}, fractional::AbstractString, ::RoundingMode{:N
     return T(0)
 end
 
+
+"""
+    max_exp10(T)
+
+The highest value of `x` which does not result in an overflow when evaluating `T(10)^x`. For
+types of `T` that do not overflow -1 will be returned.
+"""
+function max_exp10{T <: Integer}(::Type{T})
+    applicable(typemax, T) || return -1
+    W = widen(T)
+    type_max = W(typemax(T))
+
+    powt = one(W)
+    ten = W(10)
+    exponent = 0
+
+    while type_max > powt
+        powt *= ten
+        exponent += 1
+    end
+
+    exponent - 1
+end
+
+"""
+    coefficient(::Type{FD{T, f}}) -> T
+
+Compute `10^f` as an Integer without overflow. Note that overflow will not occur for any
+constructable `FD{T, f}`.
+"""
+coefficient{T, f}(::Type{FD{T, f}}) = T(10)^f
+coefficient{T, f}(fd::FD{T, f}) = coefficient(FD{T, f})
+value(fd::FD) = fd.i
+
 end