[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 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)


----
'''[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.

   

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