Sarnold 2005-12-22
Here is a package for decimal arithmetic with math::bignum in pure Tcl.
Usage:
% namespace import ::math::decimal::* % set a [fromstr 0.150] % tostr [round $a 1] 0.2 % set b [fromstr 23945.59] % tostr [- $b $a] 23945.44 % set product [* $a $b] % tostr $product 3591.83850 % tostr [round $product [precision $b]] 3591.84
The source
package require math::bignum
# this line helps when I want to source this file again and again
catch {namespace delete ::math::decimal}
# private namespace
# this software works only with Tcl v8.4 and higher
# it is using the package math::bignum
namespace eval ::math::decimal {
# cached constants
# 10 as a bignum
variable Ten
set Ten [::math::bignum::fromstr 10]
variable One
set One [::math::bignum::fromstr 1]
}
################################################################################
# procedures that handle decimal numbers
# these procedures are sorted by name (after eventually removing the underscores)
#
# Decimal numbers are stored as a list:
# {"decimal" Mantissa Exponent Power} where "decimal" is a tag
# Mantissa and Power are bignums
# Exponent is an integer
#
# The resulting value is Mantissa*10^Exponent or Mantissa*Power, both results are good.
#
################################################################################
#
# returns the absolute value
#
proc ::math::decimal::abs {number} {
lset number 1 [::math::bignum::abs [lindex $number 1]]
return $number
}
#
# returns A + B
#
proc ::math::decimal::+ {a b} {
if {[lindex $a 2]<[lindex $b 2]} {
return [+ $b $a]
}
# retrieving parameters from A and B
foreach {dummy integerA expA powA} $a {break}
foreach {dummy integerB expB powB} $b {break}
if {$expA>$expB} {
set pow [::math::bignum::div $powA $powB]
lset b 1 [::math::bignum::add [::math::bignum::div $integerA $pow] $integerB]
return $b
} elseif {$expA==$expB} {
# nothing to shift
lset a 1 [::math::bignum::add $integerA $integerB]
return $a
}
error "internal error"
}
#
# rounds to ceiling the number, with a given precision
#
proc ::math::decimal::ceil {number {precision 0}} {
foreach {dummy integer exp pow} $number {break}
if {$exp<$precision} {
error "not enough precision to perform ceil"
}
if {$exp==$precision} {
return $number
}
lset number 2 $precision
variable Ten
set divider [::math::bignum::pow $Ten [::math::bignum::fromstr [expr {$exp-$precision}]]]
lset number 3 [::math::bignum::div $pow $divider]
# saving the sign ...
set sign [::math::bignum::sign $integer]
if {$sign} {
# negative number: return the floor
lset number 1 [::math::bignum::div $integer $divider]
return $number
}
# decimal part
foreach {integer remainder} [::math::bignum::divqr $integer $divider] {break}
# fractional part
if {![::math::bignum::iszero $remainder]} {
lset number 1 [::math::bignum::add 1 $integer]
} else {
lset number 1 $integer
}
return $number
}
#
# returns 0 if A = B , otherwise returns the sign of A-B
# IN : a, b (Decimals)
#
proc ::math::decimal::compare {a b} {
if {[lindex $a 2]<[lindex $b 2]} {
return [expr {-[compare $b $a]}]
}
# retrieving data
foreach {dummy aint aexp apow} $a {break}
foreach {dummy bint bexp bpow} $b {break}
if {$aexp==$bexp} {
return [::math::bignum::compare $aint $bint]
}
set diff [expr {$aexp-$bexp}]
set aint [::math::bignum::div $aint [::math::bignum::div $apow $bpow]]
if {[::math::bignum::compare $aint $bint]>=0} {
return 1
}
return -1
}
#
# divide A by B and returns the result
# throw error : when dividing by zero
#
proc ::math::decimal::/ {a b} {
foreach {dummy integerA expA powA} $a {break}
foreach {dummy integerB expB powB} $b {break}
# dividing any number by zero throws an error
if {[::math::bignum::iszero $integerB]} {
error "divide by zero"
}
# aE-ea / bE-eb = (a/b)E(eb-ea)
# a/pa / b/pb = (a/b)/(pa/pb)
# now ea=expA, eb=expB, pa=powA and pb=powB
if {$expB>$expA} {
error "not enough precision in the numerator"
}
lset a 2 [expr {$expA-$expB}]
lset a 3 [::math::bignum::div $powA $powB]
lset a 1 [::math::bignum::div $integerA $integerB]
return $a
}
#
# procedure floor : rounds to floor with a given precision (number of decimals)
# returns a Decimal number
#
proc ::math::decimal::floor {number {precision 0}} {
return [opp [ceil [opp $number] $precision]]
}
#
# returns a list formed by an integer and a precision (decimals number)
# x = A * 10 power (-B) = A/C
# return [list "decimal" A B C] (where "decimal" is the tag)
#
proc ::math::decimal::fromstr {number} {
set number [string trimleft $number +]
set sign 0
if {[string index $number 0]=="-"} {
set sign 1
set number [string range $number 1 end]
}
# a floating-point number may have a dot
set tab [split $number .]
set exp 0
if {[llength $tab]>2} {error "syntax error in number"}
if {[llength $tab]==2} {
set number [join $tab ""]
if {![string is digit $number]} {
error "not a number"
}
# increment by the number of decimals (after the dot)
set exp [string length [lindex $tab 1]]
}
set number [::math::bignum::fromstr $number]
::math::bignum::setsign number $sign
if {!$exp} {
variable One
return [list decimal $number 0 $One]
}
variable Ten
return [list decimal $number $exp [::math::bignum::pow $Ten [::math::bignum::fromstr $exp]]]
}
#
# returns 1 if A is null, 0 otherwise
#
proc ::math::decimal::iszero {a} {
return [::math::bignum::iszero [lindex $a 1]]
}
#
# returns A modulo B (like with fmod() math function)
#
proc ::math::decimal::% {a b} {
set ratio [/ $a $b]
# examples : fmod(3,2)=1 ratio=1.5
# fmod(1,2)=1 ratio=0.5
# ratio>0 and b>0 : get floor(ratio)
# fmod(-3,-2)=-1 ratio=1.5
# fmod(-1,-2)=-1 ratio=0.5
# ratio>0 and b<0 : get floor(ratio)
# fmod(-3,2)=-1 ratio=-1.5
# fmod(-1,2)=-1 ratio=-0.5
# ratio<0 and b>0 : get ceil(ratio)
# fmod(3,-2)=1 ratio=-1.5
# fmod(1,-2)=1 ratio=-0.5
# ratio<0 and b<0 : get ceil(ratio)
if {[sign $ratio]} {
set ratio [ceil $ratio]
} else {
set ratio [floor $ratio]
}
return [- $a [* $ratio $b]]
}
#
# returns the product of A by B
#
proc ::math::decimal::* {a b} {
foreach {dummy integerA expA powA} $a {break}
foreach {dummy integerB expB powB} $b {break}
return [list decimal [::math::bignum::mul $integerA $integerB] \
[expr {$expA+$expB}] [::math::bignum::mul $powA $powB]]
}
#
# given a decimal number A, returns -A
#
proc ::math::decimal::opp {a} {
set integer [lindex $a 1]
::math::bignum::setsign integer [expr {![::math::bignum::sign $integer]}]
lset a 1 $integer
return $a
}
#
# returns the precision of some decimal number
#
proc ::math::decimal::precision {a} {
return [lindex $a 2]
}
#
# returns $number rounded to the given precision
#
proc ::math::decimal::round {number {precision 0}} {
if {$precision<0} {
error "negative second argument, cannot round"
}
foreach {dummy integer exp power} $number {break}
if {$exp<$precision} {
error "not enough precision to round"
}
if {$exp==$precision} {
return $number
}
lset number 2 $precision
variable Ten
# divider := 10^(old_precision-new_precision)/2
# := $power/10^$precision/2
# := $power/$new_power >> 1
# old_integer/divider := new_integer*2 + carry
# if carry=1 then increment new_integer
lset number 3 [set new_power [::math::bignum::pow $Ten [::math::bignum::fromstr $precision]]]
set divider [::math::bignum::rshift [::math::bignum::div $power $new_power] 1]
# saving the sign, ...
set sign [::math::bignum::sign $integer]
set integer [::math::bignum::abs $integer]
set integer [::math::bignum::div $integer $divider]
# first bit
set carry [::math::bignum::testbit $integer 0]
set integer [::math::bignum::rshift $integer 1]
if {$carry} {
set integer [::math::bignum::add 1 $integer]
}
# ... restore the sign now
::math::bignum::setsign integer $sign
lset number 1 $integer
return $number
}
#
# gets the sign of a decimal number
# we keep the bignum convention : 0 for positive, 1 for negative
#
proc ::math::decimal::sign {a} {
return [::math::bignum::sign [lindex $a 1]]
}
#
# substracts B to A
#
proc ::math::decimal::- {a b} {
return [+ $a [opp $b]]
}
#
# converts a decimal number stored as a list to a string
#
proc ::math::decimal::tostr {number} {
foreach {dummy integer exp delta} $number {break}
if {$exp<0} {
error "negative decimals number"
}
set integer [::math::bignum::tostr $integer]
set sign ""
if {[string index $integer 0]=="-"} {
set sign -
set integer [string range $integer 1 end]
}
set len [string length $integer]
if {$len<$exp+1} {
set integer [string repeat 0 [expr {$exp+1-$len}]]$integer
}
return $sign[string range $integer 0 end-$exp].[string range $integer end-[incr exp -1] end]
}
# exporting public interface
namespace eval ::math::decimal {
foreach function {
+ - * / % abs opp precision
compare fromstr tostr iszero
round ceil floor
} {
namespace export $function
}
}
package provide math::decimal 0.1