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.1AM (24 march 2010) As Tcl 8.5 supports arbitrary-precision integers out-of-the-box, the math::bignum package is not needed anymore. We can get the above functionality more directly. Here is an (independent) implementation of decimal arithmetic. Note that it is not well-tested and for serious work, it should borrow stuff from above.
# decmath.tcl --
# Simple package for dealing with decimal arithmetic.
# Numbers are represented as:
# mantisse * 10**exponent
#
# Note:
# Rounding has not been explicitly tested
#
# Decmath --
# Namespace for the decimal arithmetic procedures
#
namespace eval ::Decmath {
variable scale 10
variable maxmantisse [expr {10**$scale}]
variable bndmantisse [expr {10*$maxmantisse}]
namespace export + - / * tostr fromstr setScale
}
# setScale --
# Set the overall scale (maximum number of digits retained)
#
# Arguments:
# newscale New scale
#
# Result:
# None
#
proc ::Decmath::setScale {newscale} {
variable scale
variable maxmantisse
variable bndmantisse
set scale $newscale
set maxmantisse [expr {10**$scale}]
set bndmantisse [expr {10*$maxmantisse}]
}
# + --
# Add two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# Sum of both (rescaled)
#
proc ::Decmath::+ {a b} {
foreach {ma ea} $a {break}
foreach {mb eb} $b {break}
if { $ea > $eb } {
set ma [expr {$ma * 10 ** ($ea-$eb)}]
set er $eb
} else {
set mb [expr {$mb * 10 ** ($eb-$ea)}]
set er $ea
}
set mr [expr {$ma + $mb}]
return [Rescale $mr $er]
}
# - --
# Subtract two numbers (or unary minus)
#
# Arguments:
# a First operand
# b Second operand (optional)
#
# Result:
# Sum of both (rescaled)
#
proc ::Decmath::- {a {b {}}} {
if { $b == {} } {
return [- {0 0} $a]
} else {
lset b 0 [expr {-[lindex $b 0]}]
return [+ $a $b]
}
}
# * --
# Multiply two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# Product of both (rescaled)
#
proc ::Decmath::* {a b} {
foreach {ma ea} $a {break}
foreach {mb eb} $b {break}
set mr [expr {$ma * $mb}]
set er [expr {$ea + $eb}]
return [Rescale $mr $er]
}
# / --
# Divide two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# Quotient of both (rescaled)
#
proc ::Decmath::/ {a b} {
variable scale
variable maxmantisse
foreach {ma ea} $a {break}
foreach {mb eb} $b {break}
set mr [expr {$ma * $maxmantisse / $mb}]
set er [expr {$ea - $eb - $scale}]
return [Rescale $mr $er]
}
# Rescale --
# Rescale the number (using proper rounding)
#
# Arguments:
# mantisse Mantisse of the number
# exponent Exponent of the number
#
# Result:
# Rescaled number (as a list)
#
proc ::Decmath::Rescale {mantisse exponent} {
variable maxmantisse
variable bndmantisse
if { abs($mantisse) <= $maxmantisse } {
return [list $mantisse $exponent]
}
while { abs($mantisse) > $bndmantisse } {
set mantisse [expr {$mantisse/10}]
incr exponent
}
set rest [expr {$mantisse % 10}]
if { $rest > 5 } {
if { $mantisse > 0 } {
incr mantisse
} else {
incr mantisse -1
}
} elseif { $rest == 5 } {
if { ($mantisse/10) % 2 == 1 } {
incr mantisse
}
}
return [list $mantisse $exponent]
}
# tostr --
# Convert number to string
#
# Arguments:
# number Number to be converted
#
# Result:
# Number in the form of a string
#
proc ::Decmath::tostr {number} {
foreach {mantisse exponent} $number {break}
if { $exponent > 0 } {
set string $mantisse[string repeat 0 $exponent]
} else {
set digits [string length $mantisse]
set exponent [expr {-$exponent}]
if { $digits >= $exponent } {
set string [string range $mantisse 0 [expr {$digits-$exponent-1}]].[string range $mantisse [expr {$digits-$exponent}] end]
} else {
set string 0.[string repeat 0 [expr {$exponent-$digits}]]$mantisse
}
}
return $string
}
# fromstr --
# Convert string to number
#
# Arguments:
# string String to be converted
#
# Result:
# Number in the form of {mantisse exponent}
#
proc ::Decmath::fromstr {string} {
set pos [string first . $string]
if { $pos < 0 } {
set mantisse $string
set exponent 0
} else {
set fraction [string range $string [expr {$pos+1}] end]
set mantisse [string trimleft [string map {. ""} $string] 0]
set exponent [expr {-[string length $fraction]}]
}
return [Rescale $mantisse $exponent]
}
# main --
# Test it
#
namespace import ::Decmath::*
set a [fromstr 1.02]
set b [fromstr 2.1]
set c [fromstr 5.25]
puts "$a + $b = [+ $a $b] -- [tostr [+ $a $b]] -- (expr) [expr {1.02+2.1}]"
puts "$a * $b = [* $a $b] -- [tostr [* $a $b]] -- (expr) [expr {1.02*2.1}]"
puts "$c / $b = [/ $c $b] -- [tostr [/ $c $b]] -- (expr) [expr {5.25/2.1}]"
puts "0.001: [fromstr 0.001] -- [tostr {1 -3}]"
LV So, has anyone considered adding this to tcllib?
AM I have certainly thought about it, but it needs to mature a bit (I concocted my version last night, hardly a mature product :)). It will definitely be useful, but we need to know more about people's needs for this. (Financial computations are one obvious application area, but they pose very strict rules!)