beernutmark 2011-04-29 I have worked up an implementation of the Decimal Arithmetic standard as seen http://speleotrove.com/decimal/ for tcl 8.5 which passes nearly all the tests found on the aforementioned site (minus a few weird examples, NaN's and exponent overflow.
You can find it here Decimal Arithmetic Package for tcl 8.5
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]
}
# The following line was added by [JMN] but it is completely broken. $significand is not defined anywhere. [beernutmark]
# set rest [expr {$significand % 10}] ;#for case where $maxmantisse < abs($significand) <= $bndmantisse
while { abs($mantisse) > $bndmantisse } {
set rest [expr {$mantisse % 10}]
set mantisse [expr {$mantisse/10}]
incr exponent
}
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!)
This site seems to give useful information: http://speleotrove.com/decimal/
JMN 2010-04-04 There appears to be a bug in Decmath::Rescale.
e.g
/ [fromstr 100000] [fromstr 1.23] 81300813009 -6
As the result is of the form '813008130081300... (some exponent)' repeater (as can be seen if you choose larger scales), there should be no circumstance under which it rounds to contain a 9.
I suspect the line:
set rest [expr {$mantisse % 10}]should be moved up to be the 1st line within the while loop.
AM (6 april 2010) Hm, something to look into - thanks for reporting it. Yep, it should be moved into the while loop. Fixed it in the above code. (In the original code it was looking at the last digit that was kept, but it should look at the last one to be thrown away)
JMN 2010-04-07 I think there was still a range of values not handled - so I've added the line in again above the while loop with comment. Feel free to remove all this commentary about the bug once you've reviewed it.
beernutmark 2011-04-22 I have worked up a more complete implementation at Decimal Arithmetic Package for tcl 8.5.
arjen - 2010-04-07 02:51:26
One interesting feature of (most) implementations of decimal arithmetic that is missing from (IEEE) binary arithmetic is the concept of significant digits: 1.2 has 2 significant digits, 1.200 has four. In ordinary binary arithmetic there is no difference between the two, but with either of the above implementations, there is. (1.2 is usually interpreted as having an implicit error or uncertainty of 0.05, so meaning a number between 1.15 and 1.25; 1.200 has an implicit error of 0.0005) This should be included in the implementations, but it certainly is not in mine, at least not currently.
beernutmark - 2011-04-27
I don't think that the concept of Significant digits is "missing" from decimal arithmetic. It is not there by design. Look at the info on the main resource page on decimal arithmetic: http://speleotrove.com/decimal/decifaq4.html#signif
Specifically it states that every operation is carried out as if an infinitely precise result is obtainable. Just as in regular "school" arithmetic, in decimal arithmetic you keep all the precision. (1.001 + 1.1 = 2.101 not 2.1). In those cases where tracking and using significant digits is important then your application does so after the precise decimal arithmetic. While "significance arithmetic" is important and has its place it most certainly is not in Decimal arithmetic. If my bank decided to start using significance arithmetic when I next deposit 1 * 10^3 dollars into my account that had $232.00 in it I will be angry when it then only has $1000 after the deposit. Or if the point-of-sale package decides that when I sell 2 of an item that costs $5.25 that I should only charge $10 instead of $10.50 I will be equally upset.
jmn -2010-04-07 I happen to have been investigating significant figures a bit over the past days.. so I may as well throw in a little utility proc I've been using to help me get my head around it. I'm sure it can be done more directly and simply with a regexp or two - but iterating through the characters just fit my way of thinking about it. It would be nice to add an entry in the result dictionary for the uncertainty too.
#significant figures 'are a rather crude way of tracking precision'
# see http://scienceblogs.com/goodmath/2009/03/basics_significant_figures.php
proc sigfigs {floatstr} {
#1) all non-zero digits are considered significant
#2) zeros between non-zeros are significant.
#3) leading zeros not significant.
# - for this function - we will assume trailing zeros in numbers with no decimal point are insignificant.
# There is an ambiguity here because it could be that the zeros are significant but the value just happens to be round..
# we adopt the convention that a trailing . on the input will indicate that the previous zeros are significant.
# e.g "100" has 1 sigfigs but "100." has 3 sigfigs.
#rules when performing calculations:
#a) for multiplication and division the final answer should contain only as many sigfigs as the number with the *least* number of sigfigs.
#b) for addition and subtraction - keep the number of sigfigs in the input with the smallest number of *decimal places*
set i 0
set s 0
set in_sigfigs 0
set sigs "" ;#significant digits
set dp 0
set in_dp 0
if {[string first . $floatstr] >= 0} {
set has_radix 1 ;#radix point is a more general term for decimal point.
} else {
set has_radix 0
}
if {[string match "-*" $floatstr]} {
set floatstr [string map {- ""} $floatstr]
incr i
}
foreach c [split $floatstr ""] {
if {$c eq "0"} {
if {$in_sigfigs} {
append sigs "0"
}
} else {
if {$c ne "."} {
if {!$in_sigfigs} {
set in_sigfigs 1
set s $i
}
append sigs $c
} else {
set in_dp 1
}
}
if {$c ne "."} {
if {$in_dp} {
incr dp
}
}
incr i
}
if {!$has_radix} {
#no radix. We will consider trailing zeros to be insignificant.
set sigs [string trimright $sigs "0"]
}
return [list count [string length $sigs] startindex $s digits $sigs dp $dp]
}