I am pleased to announce that today, almost four years after the last version, a new version of Decimals.jl has been released!
Decimals.jl is a package for arbitrary-precision decimal floating point arithmetics written in Julia.
The new version, number 0.5.0, comes with a lot of changes, yet almost no new features. This release primarily fixes severe bugs in arithmetic operations, ensures compliance with the floating-point arithmetic standard, and generally stabilizes the package. The correctness of the implemented operations is now guarded by almost 5,000 verified test cases.
One of the major changes is that users can now set the precision, among other parameters, through so-called “context”. This feature is conveniently implemented via ScopedValues.
As of now, we are out of issues! Put the package under load and bring us new ones!
Most of the API is the same to what we are used to with Base floats in Julia (e.g., hash(x), round(x, sigdigits=3), or x + y works for Decimals). I am not worried about the stability (as in “interface might still change”) of this code.
What I am concerned about is the Context type, which is new. People might dislike the with_context/@with_context way of setting it. The introduction of context has consequences in the arithmetic operations too, because it influences Overflow/Underflow behavior. Right now, the operations throw an exception in case overflow occurs, but underflow is silently ignored (we should probably make this behavior parametrized through Context too). I, personally, would like to hear what the users think of all this before releasing 1.0.
Thank you! Although is is not mine. I have joined the project only recently. Many people contributed before me, and several people contributed to this release too.
Nice! Your representation is very similar to what Decimals.jl uses (we represent mantissa with a BigInt though). If you have gained any experience that you’d like to share, you are more than welcome to contribute to Decimals.jl too. For example, I’d be interested to hear what motivated you to make precision a parameter of the type, and what the consequences were.
Raising DFP to a power is difficult, so I used log10 logarithm to do so
# Integer Plus Fraction is used to perform calculations
# with base 10 logarithmns
struct IpfType
int::BigInt
d::DFP
end
@inline function takelog10ToIpfType(a::DFP)
if a.s == 1
throw(DomainError(-1,"Cannot take the log10 of a negative number in takelog10ToIpfType"))
end
prec = DFP_lengthMantissa(a)
newprec = prec + GUARD_DIGITS
numofbits = Int64( ceil(3.322 * newprec) )
v = setprecision(numofbits) do
log10( parse(BigFloat,MantissaStringInScientificFormat(a)) )
end
IpfType(a.expo,DFP(v,prec))
end
@inline function straightToIpfType(a::DFP)
frac = DFP_fraction(a)
expo = parse(BigInt,DFP_toIntegerString(a))
IpfType(expo,frac)
end
@inline function straightToIpfType(a::Integer)
frac = DFP_createZero(DFP_getDefaultPrecision())
IpfType(a,frac)
end
@inline function straightToIpfType(a::Integer,prec::Int64)
frac = DFP_createZero(prec)
IpfType(a,frac)
end
function +(a::IpfType,b::IpfType)::IpfType
c = a.d + b.d
if c.expo < 0 || DFP_isZero(c)
return IpfType(a.int + b.int,c)
end
if c.expo == 0
newexpo = a.int + b.int + c.m[1]
c = c << 1
return IpfType(newexpo,c)
end
end
function -(a::IpfType,b::IpfType)::IpfType
c = a.d - b.d
if c.s == 1
prec = DFP_lengthMantissa(c)
one = DFP_createOne(prec)
newc = one + c
return IpfType(a.int - b.int - 1,newc)
else
return IpfType(a.int - b.int,c)
end
end
function *(a::IpfType,b::IpfType)::IpfType
# Stage 1
stage1_int = a.int * b.int
# Stage 2
stringa = string(a.int)
aintprec = length(stringa)
bfracprec = DFP_lengthMantissa(b.d)
newprec = max(aintprec,bfracprec) + GUARD_DIGITS
aa = DFP(stringa,newprec)
bb = DFP_setPrecision(b.d,newprec)
rr = aa * bb
intrr = DFP_integer(rr)
stage2_int = parse(BigInt,DFP_toIntegerString(intrr))
stage2_frac = DFP_setPrecision(DFP_fraction(rr),bfracprec)
# Stage 3
afracprec = DFP_lengthMantissa(a.d)
stringb = string(b.int)
bintprec = length(stringb)
newprec = max(afracprec,bintprec) + GUARD_DIGITS
aa = DFP_setPrecision(a.d,newprec)
bb = DFP(stringb,newprec)
rr = aa * bb
intrr = DFP_integer(rr)
stage3_int = parse(BigInt,DFP_toIntegerString(intrr))
stage3_frac = DFP_setPrecision(DFP_fraction(rr),afracprec)
# Stage 4
afracprec = DFP_lengthMantissa(a.d)
bfracprec = DFP_lengthMantissa(b.d)
newprec = max(afracprec,bfracprec) + GUARD_DIGITS
aa = DFP_setPrecision(a.d,newprec)
bb = DFP_setPrecision(b.d,newprec)
rr = aa * bb
intrr = DFP_integer(rr)
stage4_int = parse(BigInt,DFP_toIntegerString(intrr))
stage4_frac = DFP_setPrecision(DFP_fraction(rr),afracprec)
# Final Stage
r = stage2_frac + stage3_frac + stage4_frac
intr = DFP_integer(r)
finalfrac = DFP_fraction(r)
finalint = stage1_int + stage2_int +
stage3_int + stage4_int +
parse(BigInt,DFP_toIntegerString(intr))
return IpfType(finalint,finalfrac)
end
function /(a::IpfType,b::IpfType)::IpfType
prec = DFP_lengthMantissa(a.d)
a_int_prec = length(string(a.int))
b_int_prec = length(string(b.int))
newprec = max(a_int_prec + 2 * (prec + GUARD_DIGITS), b_int_prec + 2 * (prec + GUARD_DIGITS))
newa = DFP(a.int,newprec) + DFP(a.d,newprec)
newb = DFP(b.int,newprec) + DFP(b.d,newprec)
newc = newa / newb
frac = DFP_setPrecision(DFP_fraction(newc),prec)
expo = parse(BigInt,DFP_toIntegerString(newc))
IpfType(expo,frac)
end
@inline function takeexp10OfIpfType(a::IpfType)
prec = DFP_lengthMantissa(a.d)
newprec = prec + GUARD_DIGITS
numofbits = Int64( ceil(3.322 * newprec) )
r = setprecision(numofbits) do
exp10( parse(BigFloat,MantissaStringInLog10Format(a.d)) )
end
temp = DFP(r,newprec)
DFP_setPrecision(DFP{newprec}(0,a.int,temp.m),prec)
end
function DFP_power(a::DFP,b::Integer)
# Check if a is an Error
if DFP_isError(a)
return a
end
prec = DFP_lengthMantissa(a)
if iszero(b)
return DFP_createOne(prec)
end
if DFP_isZero(a)
return DFP_createZero(prec)
end
if DFP_isOne(a)
return a
end
newprec = prec + GUARD_DIGITS
if DFP_isPositive(a)
aa = DFP_setPrecision(a,newprec)
ipfaa = takelog10ToIpfType(aa)
ipfbb = straightToIpfType(b,newprec)
result = ipfbb * ipfaa
# Handle the special case when the logarithmn
# is negative by taking the fractional part
# and add one to it and take the integer part
# and subtract one from it
if result.d.s == 1
f = result.d
len = DFP_lengthMantissa(f)
newf = DFP_createOne(len) + f
newint = result.int - 1
result = IpfType(newint,newf)
end
return DFP_setPrecision( takeexp10OfIpfType( result ),prec )
else
# Since it is not Zero then it must be negative
positive_a = DFP{prec}(0,a.expo,a.m)
result = DFP_power(positive_a,b)
if isodd(b)
# If the power is an odd integer then the
# result is negative
result = DFP{prec}(1,result.expo,result.m)
end
return result
end
end
function DFP_power(a::DFP{N},b::DFP{N}) where {N}
# Check if a or b is an Error
if DFP_isError(a)
return a
end
if DFP_isError(b)
return b
end
if DFP_isZero(b)
return DFP_createOne(N)
end
if DFP_isZero(a)
return DFP_createZero(N)
end
if DFP_isOne(a)
return a
end
if DFP_isPositive(a)
newprec = N + GUARD_DIGITS
aa = DFP_setPrecision(a,newprec)
bb = DFP_setPrecision(b,newprec)
ipfaa = takelog10ToIpfType(aa)
ipfbb = straightToIpfType(bb)
result = ipfbb * ipfaa
# Handle the special case when the logarithmn
# is negative by taking the fractional part
# and add one to it and take the integer part
# and subtract one from it
if result.d.s == 1
f = result.d
len = DFP_lengthMantissa(f)
newf = DFP_createOne(len) + f
newint = result.int - 1
result = IpfType(newint,newf)
end
return DFP_setPrecision( takeexp10OfIpfType( result ),N )
else
# Since it is not Zero then a must be negative
if DFP_isInteger(b)
positive_a = DFP{N}(0,a.expo,a.m)
result = DFP_power(positive_a,b)
unitdigit = DFP_getUnitDigit(b)
if unitdigit != nothing && isodd(unitdigit)
# If the power is an odd integer then the
# result is negative
result = DFP{N}(1,result.expo,result.m)
end
return result
else
# b is not an integer so the result is complex
# return a domain error as we only deal in the Real domain
return DFP_createError(N,3)
end
end # if DFP_isPositive(a)
end
Simple, I couldn’t decide how many digits of precision to make the DFP so I thought why not let the user decide and make N the number of decimal digits a part of DFP.
Then I had to figure out how to add together two DFP of different precisions and the simplest answer is to disallow it and force the user to change the precision to make both arguments (to addition) have the same precision.
This solves a lot of problems.
I don’t know what the maximum number of precision N could be. I only tested for N up to 2000.
Most users will just stick to one type. The default precision is 16 decimal places which was choosen to emulate Float64. I like my version because it allows the user to have multiple precision at the same time. Warning: a 2000 digit precision DFP takes a long time to do calculations.
I have upload the source via Git Hub so that you can download it to play around with it.
sample code:
include("/Users/ssiew/juliascript/SimpleDFP/SimpleDFPs.jl")
using .SimpleDFPs
a = DFP(12345.6789)
println(DFP_toCommonString(a))
cs = SimpleDFPs.DFP_toCommonString
println(cs(a^2))
