[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]
    }

#    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 [JMN] added a bug in Decmath::Rescale.  in the line:

======
set rest [expr {$significand % 10}] ;#for case where $maxmantisse < abs($significand) <= $bndmantisse
======

significand has not been set.  As I am just learning about decimal math I am not sure what it should be set to at this point.

[beernutmark] 2011-04-23 While working on adding proper rounding functions to the above Decmath package I discovered that it has quite a few errors when dealing with negative numbers and 0.  I changed tostr and fromstr and added several rounding functions.  In order to make cutting and pasting easier I have attached the entire corrected version below. I didn't want to replace [AM]'s work until he has a chance to look at the changes.
======
package provide Pos 1.0

# decmath.tcl --
#     Simple package for dealing with decimal arithmetic.
#     Numbers are represented as:
#         mantisse * 10**exponent
#
#

# 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]
}

# Abs --
#     Returns the Absolute Value of a number
#
# Arguments:
#     mantisse   Mantisse of the number
#     exponent   Exponent of the number
#
# Result:
#     Absolute value (as a list)
#
 proc ::Decmath::abs {mantisse exponent} {
     set mantisse [expr {abs($mantisse)}]
     return [list $mantisse $exponent]
 }


# 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]
    } elseif { $exponent == 0 } {
        set string $mantisse
    } else {
        set digits [string length [string trimleft $mantisse "-"]]
        set exponent [expr {-$exponent}]
        if { $digits > $exponent } {
            if { [string range $mantisse 0 0] == "-" } {
                set string [string range $mantisse 0 [expr {$digits-$exponent}]].[string range $mantisse [expr {$digits-$exponent+1}] end]
            } else {
                set string [string range $mantisse 0 [expr {$digits-$exponent-1}]].[string range $mantisse [expr {$digits-$exponent}] end]
            }
        } else {
            if { [string range $mantisse 0 0] == "-" } {
                set string -0.[string repeat 0 [expr {$exponent-$digits}]][string range $mantisse 1 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 } {
        if { abs($string) == 0 } {
            set mantisse 0
        } else {
            set mantisse $string
        }
        set exponent 0
    } else {
        if { [string range $string 0 0] == "-" } {
            set mantisse -[string trimleft [string map {. ""} [string trimleft $string "-"]] 0]
        } else {
            set mantisse [string trimleft [string map {. ""} $string] 0]
        }
        if { $mantisse == "" || $mantisse == "-"} {
            return [list 0 0]
        } else {
            set fraction [string range $string [expr {$pos+1}] end]
            set exponent [expr {-[string length $fraction]}]
        }
    }
    return [Rescale $mantisse $exponent]
}

# ipart --
#     Return the integer part of a Decimal Number
#
# Arguments:
#     Number in the form of {mantisse exponent}
#     
#
# Result:
#     Integer
#
proc ::Decmath::ipart { a } {
    foreach {ma ea} $a {break}
    if { $ea == 0 } {
        return $ma
    } elseif { $ea > 0 } {
        return $ma[string repeat 0 $ea]
    } else {
        if { [string length [string trimleft $ma "-"]] <= -$ea } {
            return 0
        } else {
            return [string range $ma 0 end-[expr {abs($ea)}]]
        }
    }
}

# round_half_up --
#     Round a Decimal Number using round-half-up method. 
#     This is the way that the normal tcl round() function works.
#     Bias: away from zero
#
# Arguments:
#     Number in the form of {mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {mantisse exponent}
#
proc ::Decmath::round_half_up {a digits} {
    foreach {ma ea} $a {break}
    
    if {[expr {$digits + $ea}] == 0 } {
        return $a
    } elseif { [expr {$digits + $ea}] > 0 } {
        set mantissa $ma[string        repeat 0 [expr {$digits+$ea}]]
        set exponent [expr {-1 * $digits}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer [tostr [list $ma $round_exponent]]
        set mantissa [expr {round($integer)}]
        set exponent [expr {-1 * $digits}]
    }
    return [Rescale $mantissa $exponent]
}

# round_half_even --
#     Round a Decimal Number using round-half-even method. 
#     Bias: none
#
# Arguments:
#     Number in the form of {mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {mantisse exponent}
#
proc ::Decmath::round_half_even {a digits} {
    foreach {ma ea} $a {break}
    if {[expr {$digits + $ea}] == 0 } {
        return $a
    } elseif { [expr {$digits + $ea}] > 0 } {
        set mantissa $ma[string        repeat 0 [expr {$digits+$ea}]]
        set exponent [expr {-1 * $digits}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer_part [ipart [list $ma $round_exponent]]
        if { [tostr [- [list [string trimleft $ma "-"] $round_exponent] [list [string trimleft ${integer_part}5 "-"] -1]]] == 0} {
            # We are rounding something with fractional part .5
            if { [expr {$integer_part % 2}] } {
                # We are odd so round up
                if { $integer_part > 0 } {
                    set mantissa [expr {$integer_part + 1}]
                } else {
                    set mantissa [expr {$integer_part - 1}]
                }
            } else {
                # We are even so round down
                set mantissa $integer_part
            }
        } else {
            set integer [tostr [list $ma $round_exponent]]
            set mantissa [expr {round($integer)}]
        }
        set exponent [expr {-1 * $digits}]
    }
    return [Rescale $mantissa $exponent]
}

# round_half_down --
#     Round a Decimal Number using round-half-down method. 
#     Bias: towards zero
#
# Arguments:
#     Number in the form of {mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {mantisse exponent}
#
proc ::Decmath::round_half_down {a digits} {
    foreach {ma ea} $a {break}
    
    if {[expr {$digits + $ea}] == 0 } {
        return $a
    } elseif { [expr {$digits + $ea}] > 0 } {
        set mantissa $ma[string        repeat 0 [expr {$digits+$ea}]]
        set exponent [expr {-1 * $digits}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer_part [ipart [list $ma $round_exponent]]
        if { [tostr [- [list [string trimleft $ma "-"] $round_exponent] [list [string trimleft ${integer_part}5 "-"] -1]]] == 0} {
            # We are rounding something with fractional part .5
            # The rule is to round half down.
            set mantissa $integer_part
        } else {
            set integer [tostr [list $ma $round_exponent]]
            set mantissa [expr {round($integer)}]
        }
        set exponent [expr {-1 * $digits}]
    }
    return [Rescale $mantissa $exponent]
}

proc ::Decmath::floor {a digits} {
    foreach {ma ea} $a {break}
    if { -$ea== $digits } {
        return $a
    } elseif { [expr {$digits + $ea}] > 0 } {
        set mantissa $ma[string        repeat 0 [expr {$digits+$ea}]]
        set exponent [expr {-1 * $digits}]
    } else {
        set round_exponent [expr {abs($digits + $ea)}]
        if { [string length [string trimleft $ma "-"]] <= $round_exponent } {
            set mantissa 0
        } else {
            set mantissa [string range $ma 0 end-$round_exponent]
        }
        set exponent [expr {-1 * $digits}]
    }
    return [Rescale $mantissa $exponent]
}


======
The usefulness of this Decmath package can be seen in the following:
======
% Decmath::floor [Decmath::- [Decmath::fromstr 8.2] [Decmath::fromstr 0.2]] 0
8 0
% expr {floor(8.2 - 0.2)}
7.0
======

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

[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]
    }


<<categories>>Mathematics