[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
  
----
[AM] (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!)

----
!!!!!!
%| [Category Mathematics] |%
!!!!!!