Context

Author: barucden

src/Decimals.jl

@@ -21,6 +21,7 @@ struct Decimal <: AbstractFloat
 end
 
 include("bigint.jl")
+include("context.jl")
 
 # Convert between Decimal objects, numbers, and strings
 include("decimal.jl")

src/arithmetic.jl

@@ -4,7 +4,8 @@ Base.promote_rule(::Type{Decimal}, ::Type{<:Real}) = Decimal
 Base.promote_rule(::Type{BigFloat}, ::Type{Decimal}) = Decimal
 Base.promote_rule(::Type{BigInt}, ::Type{Decimal}) = Decimal
 
-const BigTen = BigInt(10)
+Base.:(+)(x::Decimal) = fix(x)
+Base.:(-)(x::Decimal) = fix(Decimal(!x.s, x.c, x.q))
 
 # Addition
 # To add, convert both decimals to the same exponent.
@@ -24,9 +25,6 @@ function Base.:(+)(x::Decimal, y::Decimal)
     return normalize(Decimal(s, abs(c), y.q))
 end
 
-# Negation
-Base.:(-)(x::Decimal) = Decimal(!x.s, x.c, x.q)
-
 # Subtraction
 Base.:(-)(x::Decimal, y::Decimal) = +(x, -y)
 

src/bigint.jl

@@ -4,6 +4,8 @@ else
     const libgmp = Base.GMP.libgmp
 end
 
+const BigTen = BigInt(10)
+
 function isdivisible(x::BigInt, n::Int)
     r = ccall((:__gmpz_divisible_ui_p, libgmp), Cint,
               (Base.GMP.MPZ.mpz_t, Culong), x, n)

src/context.jl

@@ -0,0 +1,106 @@
+Base.@kwdef mutable struct Context
+    precision::Int
+    rounding::RoundingMode
+    Emax::Int
+    Emin::Int
+end
+
+const CONTEXT = Context(precision=28,
+                        rounding=RoundNearest,
+                        Emax=999999,
+                        Emin=-999999)
+
+function Base.setprecision(::Type{Decimal}, precision::Int)
+    CONTEXT.precision = precision
+    return precision
+end
+
+Base.precision(::Type{Decimal}) = CONTEXT.precision
+
+function Base.setrounding(::Type{Decimal}, rounding::RoundingMode)
+    CONTEXT.rounding = rounding
+    return rounding
+end
+
+Base.rounding(::Type{Decimal}) = CONTEXT.rounding
+
+"""
+    fix(x)
+
+Round and fix the exponent of `x` to keep it within the precision and exponent
+limits as given by the current `CONTEXT`.
+"""
+function fix(x::Decimal)
+    prec = precision(Decimal)
+    rmod = rounding(Decimal)
+
+    Emin, Emax = CONTEXT.Emin, CONTEXT.Emax
+    Etiny = Emin - prec + 1
+    Etop = Emax - prec + 1
+
+    if iszero(x)
+        return Decimal(x.s, x.c, clamp(x.q, Etiny, Etop))
+    end
+
+    clen = ndigits(x.c)
+    exp_min = clen + x.q - prec
+
+    # Equivalent to `clen + x.q - 1 > Emax`
+    if exp_min > Etop
+        throw(OverflowError("Exponent limit ($Emax) exceeded: $x"))
+    end
+
+    subnormal = exp_min < Etiny
+    if subnormal
+        exp_min = Etiny
+    end
+
+    # Number of digits and exponent within bounds
+    if x.q ≥ exp_min
+        return x
+    end
+
+    # Signed coefficient for rounding modes like RoundToZero
+    c = (-1)^x.s * x.c
+    q = exp_min
+
+    # Number of digits of the resulting coefficient
+    digits = clen + x.q - exp_min
+    if digits < 0
+        c = big(1)
+        q = exp_min - 1
+        digits = 0
+    end
+
+    # Number of least significant digits to remove from `c`
+    trun_len = clen - digits
+
+    # Split `c` into `digits` most significant digits and `trun_len` least
+    # significant digits
+    # This is like round(c, rmod, sigdigits=digits), except here we can
+    # tell from `rem` if the rounding was lossless
+    c, rem = divrem(c, BigTen ^ trun_len, rmod)
+
+    # Rounding is exact if the truncated digits were zero
+    exact = iszero(rem)
+
+    # If the number of digits exceeded `digits` after rounding,
+    # it means that `c` was like 99...9 and was rounded up,
+    # becoming                  100...0, so `c` is divisible by 10
+    if ndigits(c) > prec
+        c = exactdiv(c, 10)
+        q += 1
+    end
+
+    # Exponent might have exceeded due to rounding
+    if q > Etop
+        throw(OverflowError("Exponent limit ($Emax) exceeded: $x"))
+    end
+
+    if subnormal && !exact
+        #            throw(ErrorException("Underflow"))
+    end
+
+    return Decimal(signbit(c), abs(c), q)
+end
+