Export sigexp (#138)

Author: jmkuhn

src/DecFP.jl

@@ -5,7 +5,7 @@ using DecFP_jll
 
 import Printf, Random, SpecialFunctions
 
-export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str, exponent10, ldexp10
+export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str, exponent10, ldexp10, sigexp
 
 const _buffer = Vector{Vector{UInt8}}()
 
@@ -238,6 +238,19 @@ function isnanstr(s::AbstractString)
     return false
 end
 
+"""
+    sigexp(x::DecFP.DecimalFloatingPoint)
+
+Return `(sign, significand, exponent)` such that `x` is equal to `sign * significand * 10^exponent`.
+Throws `DomainError` for infinite or `NaN` arguments.
+
+# Examples
+```jldoctest
+julia> sigexp(Dec64(1.25))
+(1, 125, -2)
+```
+"""
+sigexp(x::DecimalFloatingPoint)
 
 """
     exponent10(x::DecFP.DecimalFloatingPoint)
@@ -268,6 +281,7 @@ ldexp10(x::DecFP.DecimalFloatingPoint, n::Integer)
 for w in (32,64,128)
     BID = Symbol(string("Dec",w))
     Ti = eval(Symbol(string("UInt",w)))
+    Tsi = eval(Symbol(string("Int",w)))
     T = eval(BID)
 
     # hack: we need an internal parsing function that doesn't check exceptions, since
@@ -439,6 +453,26 @@ for w in (32,64,128)
             return Int32(n), Int32(pt), neg
         end
 
+        function sigexp(x::$BID)
+            isnan(x) && throw(DomainError(x, "sigexp(x) is not defined for NaN."))
+            isinf(x) && throw(DomainError(x, "sigexp(x) is only defined for finite x."))
+            p = 9 * $w ÷ 32 - 2
+            emax = 3 * 2^($w ÷ 16 + 3)
+            bias = emax + p - 2
+            ebits = $w ÷ 16 + 6
+            sbits = $w - ebits - 1
+            n = reinterpret($Ti, x)
+            if n & $Ti(0x3) << ($w - 3) == $Ti(0x3) << ($w - 3)
+                s = ((n & typemax($Ti) >> (ebits + 3)) | one($Ti) << sbits) % $Tsi
+                e = ((n & typemax($Ti) << (sbits + 1) >> 3) >> (sbits - 2)) % Int
+            else
+                s = (n & typemax($Ti) >> (ebits + 1)) % $Tsi
+                e = ((n & typemax($Ti) << (sbits + 1) >> 1) >> sbits) % Int
+            end
+            e -= bias
+            return signbit(x) ? -1 : 1, s, e
+        end
+
         Base.fma(x::$BID, y::$BID, z::$BID) = nox(ccall(($(bidsym(w,"fma")), libbid), $BID, ($BID,$BID,$BID,Cuint,Ref{Cuint}), x, y, z, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
         Base.muladd(x::$BID, y::$BID, z::$BID) = fma(x,y,z) # faster than x+y*z
 
@@ -583,7 +617,7 @@ for (f) in (:trunc, :floor, :ceil)
     @eval function Base.$f(::Type{I}, x::DecimalFloatingPoint) where {I<:Integer}
         x′ = $f(x)
         typemin(I) <= x′ <= typemax(I) || throw(InexactError(Symbol($f), I, x))
-        s, e = sigexp(x′)
+        _, s, e = sigexp(x′)
         return I(flipsign(s * I(10)^e, x))
     end
 end
@@ -595,7 +629,7 @@ Base.convert(::Type{Integer}, x::DecimalFloatingPoint) = convert(Int, x)
 function Base.convert(::Type{I}, x::DecimalFloatingPoint) where {I<:Integer}
     x != trunc(x) && throw(InexactError(:convert, I, x))
     typemin(I) <= x <= typemax(I) || throw(InexactError(:convert, I, x))
-    s, e = sigexp(x)
+    _, s, e = sigexp(x)
     return I(flipsign(s * I(10)^e, x))
 end
 
@@ -641,12 +675,12 @@ function Base.decompose(x::DecimalFloatingPoint)::Tuple{BigInt, Int, BigInt}
     isnan(x) && return 0, 0, 0
     isinf(x) && return ifelse(signbit(x), -1, 1), 0, 0
     iszero(x) && return 0, 0, ifelse(signbit(x), -1, 1)
-    s, e = sigexp(x)
+    sign, s, e = sigexp(x)
     if e >= 0
         if e <= 27
-            return s * BigInt(Int64(5)^e), e, ifelse(signbit(x), -1, 1)
+            return s * BigInt(Int64(5)^e), e, sign
         else
-            return s * BigInt(5)^e, e, ifelse(signbit(x), -1, 1)
+            return s * BigInt(5)^e, e, sign
         end
     else
         e2 = -e
@@ -664,45 +698,6 @@ function Base.decompose(x::DecimalFloatingPoint)::Tuple{BigInt, Int, BigInt}
     end
 end
 
-function sigexp(x::Dec32)
-    n = reinterpret(UInt32, x)
-    if n & 0x60000000 == 0x60000000
-        s = ((n & 0x001fffff) | 0x00800000) % Int32
-        e = ((n & 0x1fe00000) >> 21) % Int
-    else
-        s = (n & 0x007fffff) % Int32
-        e = ((n & 0x7f800000) >> 23) % Int
-    end
-    e -= 101
-    return s, e
-end
-
-function sigexp(x::Dec64)
-    n = reinterpret(UInt64, x)
-    if n & 0x6000000000000000 == 0x6000000000000000
-        s = ((n & 0x0007ffffffffffff) | 0x0020000000000000) % Int64
-        e = ((n & 0x1ff8000000000000) >> 51) % Int
-    else
-        s = (n & 0x001fffffffffffff) % Int64
-        e = ((n & 0x7fe0000000000000) >> 53) % Int
-    end
-    e -= 398
-    return s, e
-end
-
-function sigexp(x::Dec128)
-    n = reinterpret(UInt128, x)
-    if n & 0x60000000000000000000000000000000 == 0x60000000000000000000000000000000
-        s = ((n & 0x00007fffffffffffffffffffffffffff) | 0x00020000000000000000000000000000) % Int128
-        e = ((n & 0x1fff8000000000000000000000000000) >> 111) % Int
-    else
-        s = (n & 0x0001ffffffffffffffffffffffffffff) % Int128
-        e = ((n & 0x7ffe0000000000000000000000000000) >> 113) % Int
-    end
-    e -= 6176
-    return s, e
-end
-
 function Base.:(==)(dec::DecimalFloatingPoint, rat::Rational)
     if isfinite(dec)
         rat.den == 1 && return dec == rat.num

test/runtests.jl

@@ -183,6 +183,16 @@ for T in (Dec32, Dec64, Dec128)
     x,y,z = 1.5, -3.25, 0.0625 # exactly represented in binary
     xd = T(x); yd = T(y); zd = T(z)
 
+    @test_throws DomainError sigexp(T(NaN))
+    @test_throws DomainError sigexp(T(Inf))
+    @test_throws DomainError sigexp(T(-Inf))
+    @test sigexp(T(0)) == (1, 0, 0)
+    @test sigexp(T("-0")) == (-1, 0, 0)
+    @test sigexp(T(1)) == (1, 1, 0)
+    @test sigexp(T(-1)) == (-1, 1, 0)
+    @test sigexp(T(-1.25)) == (-1, 125, -2)
+    @test sigexp(maxintfloat(T) - 1) == (1, Int128(maxintfloat(T) - 1), 0)
+
     @test fma(xd,yd,zd)::T == muladd(xd,yd,zd)::T == xd*yd+zd == fma(x,y,z)
 
     @test one(T)::T == 1