I prototyped a domain implementing "decimal floating point" with very basic arithmetics, which is similar to REXX's one and Java's BigDecimal class, and I've taken some advices from the paper http://www.speleotrove.com/decimal/IEEE-cowlishaw-arith16.pdf.
One feature of decimal float over binary is the printed representation (is supposed to be) exactly matches the internal one. no more 1.2-1.1 ~= 0.1.
(1) -> a := float(12,-1)$Decimal
(1) 12E-1
Type: Decimal
(2) -> b := float(11,-1)$Decimal
(2) 11E-1
Type: Decimal
(3) -> c := float(1,-1)$Decimal
(3) 1E-1
Type: Decimal
(4) -> (a - b = c)@Boolean
(4) true
Type: Boolean
(5) -> (1.2 - 1.1 = 0.1)@Boolean
(5) false
Type: Boolean
)abbrev domain DECIMAL Decimal
B ==> Boolean
I ==> Integer
S ==> String
PI ==> PositiveInteger
RN ==> Fraction Integer
SF ==> DoubleFloat
N ==> NonNegativeInteger
++ Author: LdBeth
++ Date Created: November 2021
++ Basic Operations:
++ Keywords: float, floating point, number, decimal
++ Description: \spadtype{Decimal} implements arbitrary precision floating
++ point arithmetic.
Decimal():
Join(FloatingPointSystem,
ConvertibleTo InputForm,
arbitraryPrecision, arbitraryExponent)
== add
Rep := Record( mantissa : I, exponent : I )
DIGS : PI := 20
inc ==> increasePrecision
dec ==> decreasePrecision
LENGTH(n) ==> ilog10p1(n)
-- local utility operations
ilog10p1 : I -> I -- integer log10 plus 1
chop : (%, PI) -> % -- chop x at p bits of precision
power : (%, I) -> % -- x ^ n with chopping
plus : (%, %) -> % -- addition with no rounding
sub : (%, %) -> % -- subtraction with no rounding
negate : % -> % -- negation with no rounding
times : (%, %) -> % -- multiply x and y with no rounding
dvide : (%, %) -> % -- divide x by y with special normalizing
ceillog10base2 : PI -> PI -- rational approximation
ilog10p1(x) ==
i : N := (4004 * ((length(x)$I)-1) quo 13301) :: N
n : I := abs(x) quo (10^i)$I
while n > 0 repeat
n := n quo 10
i := i + 1
i
chop(x, p) ==
e : I := LENGTH x.mantissa - p
if e > 0 then x := [(x.mantissa quo (10^e::N)), x.exponent+e]
x
float(m, e) == [m, e]
float(m, e, b) ==
m = 0 => 0
inc 1; r := m * [b, 0] ^ e; dec 1
round r
order(a) == LENGTH(a.mantissa) + a.exponent - 1
ceillog10base2 n == ((13301 * n + 4003) quo 4004) :: PI
digits() == DIGS
digits(n) == (t := DIGS; DIGS := n; t)
precision() == digits()
precision(n) == digits(n)
bits() == 1 + ceillog10base2 digits()
bits(n) == (t := bits(); digits (max(1, 4004 * (n-1) quo 13301)::PI); t)
increasePrecision n == (b := digits(); digits((b + n)::PI); b)
decreasePrecision n == (b := digits(); digits((b - n)::PI); b)
round x ==
m := x.mantissa
m = 0 => 0
e := LENGTH m - digits()
if e > 0 then
y := m quo 10^((e-1)::N) +5
y := y quo 10
if LENGTH y > digits() then
y := y quo 10
e := e+1
x := [y, x.exponent+e]
x
0 == [0, 0]
1 == [1, 0]
base() == 10
mantissa x == x.mantissa
exponent x == x.exponent
one? a == a = 1
zero? a == zero?(a.mantissa)
negative? a == negative?(a.mantissa)
positive? a == positive?(a.mantissa)
x = y ==
x.exponent = y.exponent =>
x.mantissa = y.mantissa
order x = order y and sign x = sign y and zero? (x - y)
x < y ==
y.mantissa = 0 => x.mantissa < 0
x.mantissa = 0 => y.mantissa > 0
negative? x and positive? y => true
negative? y and positive? x => false
order x < order y => positive? x
order x > order y => negative? x
negative? (x-y)
abs x == if negative? x then -x else round x
sign x == sign x.mantissa
- x == negate x
negate x == [-x.mantissa, x.exponent]
x + y == round plus(x, y)
x - y == round plus(x, negate y)
plus(x, y) ==
mx := x.mantissa; my := y.mantissa
mx = 0 => y
my = 0 => x
ex := x.exponent; ey := y.exponent
ex = ey => [mx+my, ex]
de := ex + LENGTH mx - ey - LENGTH my
de > digits()+1 => x
de < -(digits()+1) => y
if ex < ey then (mx, my, ex, ey) := (my, mx, ey, ex)
mw := my + mx * (10^(ex-ey)::N)
[mw, ey]
x : % * y : % == round times(x, y)
x : I * y : % ==
if LENGTH x > digits() then round [x, 0] * y
else round [x * y.mantissa, y.exponent]
x : % / y : % == dvide(x, y)
x : % / y : I ==
if LENGTH y > digits() then x / round [y, 0] else x / [y, 0]
times(x : %, y : %) == [x.mantissa * y.mantissa, x.exponent + y.exponent]
dvide(x, y) ==
(q, r) := divide(x.mantissa, y.mantissa)
if r = 0 then
[q, x.exponent - y.exponent]
else
ew :I := digits() + 1
mw := (r * 10^(ew::N)) quo y.mantissa
-- possible to use binary search for optimization
while (mw rem 10) = 0 repeat
mw := mw quo 10
ew := ew - 1
[q, x.exponent - y.exponent] +
[mw, x.exponent - y.exponent - ew]
power(x, n) ==
y : % := 1; z : % := x
repeat
if odd? n then y := chop( times(y, z), digits() )
if (n := n quo 2) = 0 then return y
z := chop( times(z, z), digits() )
x : % ^ n : I ==
x = 0 =>
n = 0 => 1
n < 0 => error "division by 0"
0
n = 0 => 1
n = 1 => x
x = 1 => 1
p := LENGTH(n) + 2
inc p
y := power(x, abs n)
if n < 0 then y := dvide(1, y)
dec p
round y
--very basic formating
convert(f) : S ==
concat(concat(convert(mantissa f)@S, "E"), convert(exponent f)@S)
coerce(f) : OutputForm ==
message(convert(f)@S)
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent