Decimal Arithmetic Package for tcl 8.5

The latest version of this code is part of tcllib.

beernutmark 2011-04-29

This is a Tcl implementation of General Decimal Arithmetic as defined by the IEEE 754 standard as given on http://speleotrove.com/decimal/ .

Decimal numbers are defined as a Tcl list of sign, mantissa, and exponent.

After the source for the package, there is a test routine that is used with the language-independent test suites: http://speleotrove.com/decimal/dectest.html for the extended implementation and http://speleotrove.com/decimal/dectest0.html for the simplified implementation.

See the source header for the operations included in the package. I will continue work on it later, but for now it suits all my needs, and I figured I had better get it posted somewhere before I forget.


Examples

% namespace import ::math::decimal::*
% set a 8.20
% set b .20
% # First we must convert our numbers to decimal format.
% set a [fromstr $a]
0 820 -2
% # a is converted to decimal sign = 0, mantissa = 820, exponent = -2.
% set b [fromstr $b]
0 20 -2
% set c [* $a $b]
0 16400 -4
% # To convert back to a "regular" number to tostr.
%  puts [tostr $c]
1.6400
% # Notice that we kept the correct number of digits to show full precision. Just like you learned in elementary school.
% # $a had two digits after the decimal point and so did $b. So our result will have 2 + 2 = 4 digits after the decimal.
% set d [/ $a $b]
0 41 0
% puts [tostr $d]
41
% # In division, we want our exact result when possible.  $b went into $a exactly 41 times not 41.0 or 40.9999999.
% set e [- $a $b]
0 800 -2
% puts [tostr $e]
8.00
% # Try that using expr and you will get 7.99999999999999
% set f [tostr [* $a [+ $b $c]]]
15.088000 
% set f [/ $a $c]
0 5 0
% set f [tostr [round_half_even [* $a [+ $b $c]] 2]]
15.09
% # We rounded our initial result 15.088000 to two digits. There are 8 rounding modes to choose from.
% # Rounding will take place automatically if the digits in our mantissa go past precision. 
% # Otherwise no rounding will ever take place.

beernutmark 2011-04-30

I am now working on making this package work cleanly with Creating your own expr command written by Arjen Markus. (Perhaps not coincidentally the author of the Decmath package which was the initial inspiration and outline for this package).

To do this I am rewriting the procedures in the package to accept operands in either decimal format [list sign mantissa exponent] or string format and return the result in the same format as input. This will allow the creation of a dexpr command which will drop in replace expr when decimal arithmetic is desired.

This way one could simply use [dexpr {8.20 - 0.20}] --> 8.00 instead of [expr {8.20 - 0.20}] --> 7.999999999999999

Currently you have to write that expression as

[tostr [- [fromstr 8.2] [fromstr 0.2]]]

Not nearly as easy to read.

With this change tcl code could then be easily changed to decimal arithmetic with a simple sed: -i 's/expr/dexpr/g'

Where speed is of the essence you could still do the conversion to decimal numbers once, do all your decimal calculations, and then convert back to string at the end.


Code

package require Tcl 8.5
package provide math::decimal 1.0.2
#
# Copyright 2011 Mark Alston. All rights reserved.
#
# Redistribution and use in source and binary forms, with or 
# without modification, are permitted provided that the following 
# conditions are met:
#
#   1. Redistributions of source code must retain the above copyright 
#      notice, this list of conditions and the following disclaimer.
#
#   2. Redistributions in binary form must reproduce the above copyright 
#      notice, this list of conditions and the following disclaimer in 
#      the documentation and/or other materials provided with the distribution.
#
#   THIS SOFTWARE IS PROVIDED BY Mark Alston ``AS IS'' AND ANY EXPRESS 
#   OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
#   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
#   ARE DISCLAIMED. IN NO EVENT SHALL Mark Alston OR CONTRIBUTORS 
#   BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
#   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
#   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
#   OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
#   WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 
#   OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
#   EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# decimal.tcl --
#
#     Tcl implementation of a General Decimal Arithmetic as defined 
#     by the IEEE 754 standard as given on http:://speleotrove.com/decimal
#
#     Decimal numbers are defined as a list of sign mantissa exponent
#
#     The following operations are current implemented:
#
#       fromstr tostr  -- for converting to and from decimal numbers.
#
#       add subtract divide multiply abs compare  -- basic operations
#       max min plus minus copynegate copysign is-zero is-signed
#       is-NaN is-infinite is-finite
#
#       round_half_even round_half_up round_half_down   -- rounding methods
#       round_down round_up round_floor round_ceiling 
#       round_05up
#          
#     By setting the extended variable to 0 you get the behavior of the decimal 
#     subset arithmetic X3.274 as defined on 
#     http://speleotrove.com/decimal/dax3274.html#x3274
#
#     This package passes all tests in test suites:
#           http://speleotrove.com/decimal/dectest.html
#      and  http://speleotrove.com/decimal/dectest0.html
#
#      with the following exceptions:
#
#     This version fails some tests that require setting the max 
#     or min exponent to force truncation or rounding.  
#  
#     This version fails some tests which require the sign of zero to be set
#     correctly during rounding     
#
#     This version cannot handle sNaN's (Not sure that they are of any use for 
#     tcl programmers anyway.
#
#     If you find errors in this code please let me know at 
#         mark at beernut dot com
#
# Decimal --
#     Namespace for the decimal arithmetic procedures
#
namespace eval ::math::decimal {
    variable precision 20
    variable maxExponent 999
    variable minExponent -998
    variable tinyExponent [expr {$minExponent - ($precision - 1)}]
    variable rounding half_up
    variable extended 1

    # Some useful variables to set.
    variable zero [list 0 0 0]
    variable one [list 0 1 0]
    variable ten [list 0 1 1]
    variable onehundred [list 0 1 2]
    variable minusone [list 1 1 0]

    namespace export tostr fromstr setVariable getVariable\
                     add + subtract - divide / multiply * \
                     divide-int  remainder \
                     fma fused-multiply-add \
                     plus minus copynegate negate copysign \
                     abs compare max min \
                     is-zero is-signed is-NaN is-infinite is-finite \
                     round_half_even round_half_up round_half_down \
                     round_down round_up round_floor round_ceiling round_05up

}

# setVariable
#     Set the desired variable
#
# Arguments:
#     variable setting
#
# Result:
#     None
#
proc ::math::decimal::setVariable {variable setting} {
    variable rounding
    variable precision
    variable extended
    variable maxExponent
    variable minExponent
    variable tinyExponent

    switch -nocase -- $variable {
        rounding {set rounding $setting}
        precision {set precision $setting}
        extended {set extended $setting}
        maxExponent {set maxExponent $setting}
        minExponent {
            set minExponent $setting
            set tinyExponent [expr {$minExponent - ($precision - 1)}]
        }
        default {}
    }
}

# setVariable
#     Set the desired variable
#
# Arguments:
#     variable setting
#
# Result:
#     None
#
proc ::math::decimal::getVariable {variable} {
    variable rounding
    variable precision
    variable extended
    variable maxExponent
    variable minExponent

    switch -- $variable {
        rounding {return $rounding}
        precision {return $precision}
        extended {return $extended}
        maxExponent {return $maxExponent}
        minExponent {return $minExponent}
        default {}
    }
}

# add or +
#     Add two numbers
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     Sum of both (rescaled)
#
proc ::math::decimal::add {a b {rescale 1}} {
    return [+ $a $b $rescale]
}

proc ::math::decimal::+ {a b {rescale 1}} {
    variable extended
    variable rounding
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if {!$extended} {
        if {$ma == 0 } {
            return $b
        } 
        if {$mb == 0 } {
            return $a
        }
    }

    if { $ma eq "NaN" || $mb eq "NaN" } {
        return [list 0 "NaN" 0]
    } 

    if { $ma eq "Inf" || $mb eq "Inf" } {
        if { $ma ne "Inf" } {
            return $b
        } elseif { $mb ne "Inf" } {
            return $a
        } elseif { $sb != $sa } {
            return [list 0 "NaN" 0]
        } else {
            return $a
        }
    }
        
    if { $ea > $eb } {
        set ma [expr {$ma * 10 ** ($ea-$eb)}]
        set er $eb
    } else {
        set mb [expr {$mb * 10 ** ($eb-$ea)}]
        set er $ea
    }
    if { $sa == $sb } {
        # Both are either postive or negative
        # Sign remains the same.
        set mr [expr {$ma + $mb}]
        set sr $sa
    } else {
        # one is negative and one is positive.
        # Set sign to the same as the larger number
        # and subract the smaller from the larger.
        if { $ma > $mb } {
            set sr $sa
            set mr [expr {$ma - $mb}]
        } elseif { $mb > $ma } {
            set sr $sb
            set mr [expr {$mb - $ma}]
        } else {
            if { $rounding == "floor" } {
                set sr 1
            } else {
                set sr 0
            }
            set mr 0
        } 
    }
    if { $rescale } {
        return [Rescale [list $sr $mr $er]]
    } else {
        return [list $sr $mr $er]
    }
}

# copynegate --
#     Takes one operand and returns a copy with the sign inverted.
#     In this implementation it works nearly the same as minus
#     but is probably much faster. The main difference is that no 
#     rescaling is done.
#
# 
# Arguments:
#     a          operand
#
# Result:
#     a with sign flipped
#
proc ::math::decimal::negate { a } {
    return [copynegate $a]
}

proc ::math::decimal::copynegate { a } {
    lset a 0 [expr {![lindex $a 0]}]
    return $a
}

# copysign --
#     Takes two operands and returns a copy of the first with the
#     sign set to the sign of the second.
#
# 
# Arguments:
#     a          operand
#     b          operand
#
# Result:
#     b with a's sign
#
proc ::math::decimal::copysign { a b } {
    lset a 0 [lindex $b 0]
    return $a
}

# minus --
#     subtract 0 $a
#
#     Note: does not pass all tests on extended mode.
# 
# Arguments:
#     a          operand
#
# Result:
#     0 - $a
#
proc ::math::decimal::minus { a } {
    return [- [list 0 0 0] $a]
}

# plus --
#     add 0 $a
#
#    Note: does not pass all tests on extended mode.
# 
# Arguments:
#     a          operand
#
# Result:
#     0 + $a
#
proc ::math::decimal::plus {a} {
    return [+ [list 0 0 0] $a]
}



# subtract or -
#     Subtract two numbers (or unary minus)
#
# Arguments:
#     a          First operand
#     b          Second operand (optional)
#
# Result:
#     Sum of both (rescaled)
#
proc ::math::decimal::subtract {a {b {}} {rescale 1}} {
    return [- $a $b]
}

proc ::math::decimal::- {a {b {}} {rescale 1}} {
    variable extended

    if {!$extended} {
        foreach {sa ma ea} $a {break}
        foreach {sb mb eb} $b {break}
        if {$ma == 0 } {
            lset b 0 [expr {![lindex $b 0]}]
            return $b
        } 
        if {$mb == 0 } {
            return $a
        }
    }

    if { $b == {} } {    
        lset a 0 [expr {![lindex $a 0]}]
        return $a
    } else {
        lset b 0 [expr {![lindex $b 0]}]
        return [+ $a $b $rescale]
    }
}


# compare
#     Compare two numbers.
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     1 if a is larger than b
#     0 if a is equal to b
#    -1 if a is smaller than b.
#
proc ::math::decimal::compare {a b} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $sa != $sb } {
        if {$ma != 0 } { 
            set ma 1
            set ea 0
        } elseif { $mb != 0 } {
            set mb 1
            set eb 0
        } else {
            return 0
        }
    }
    if { $ma eq "Inf" && $mb eq "Inf" } {
        if { $sa == $sb } {
            return 0
        } elseif { $sa > $sb } {
            return -1
        } else {
            return 1
        }
    }

    set comparison [- [list $sa $ma $ea] [list $sb $mb $eb] 0] 

    if { [lindex $comparison 0] && [lindex $comparison 1] != 0 } {
        return -1
    } elseif { [lindex $comparison 1] == 0 } {
        return 0
    } else {
        return 1
    }
}

# min
#     Return the smaller of two numbers
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     smaller of a or b
#
proc ::math::decimal::min {a b} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $sa != $sb } {
        if {$ma != 0 } { 
            set ma 1
            set ea 0
        } elseif { $mb != 0 } {
            set mb 1
            set eb 0
        } 
    }
    if { $ma eq "Inf" && $mb eq "Inf" } {
        if { $sa == $sb } {
            return [list $sa "Inf" 0]
        } else {
            return [list 1 "Inf" 0]
        }
    }

    set comparison [compare [list $sa $ma $ea] [list $sb $mb $eb]] 

    if { $comparison == 1 } {
        return [Rescale $b]
    } elseif { $comparison == -1 } {
        return [Rescale $a]
    } elseif { $sb != $sa } {
        if { $sa } {
            return [Rescale $a]
        } else {
            return [Rescale $b]
        }
    } elseif { $sb && $eb > $ea } {
        # Both are negative and the same numerically. So return the one with the largest exponent.
        return [Rescale $b]
    } elseif { $sb }  {
        # Negative with $eb < $ea now.
        return [Rescale $a]
    } elseif { $ea > $eb } {
        # Both are positive so return the one with the smaller
        return [Rescale $b]
    } else {
        return [Rescale $a]
    }
}

# max
#     Return the larger of two numbers
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     larger of a or b
#
proc ::math::decimal::max {a b} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $sa != $sb } {
        if {$ma != 0 } { 
            set ma 1
            set ea 0
        } elseif { $mb != 0 } {
            set mb 1
            set eb 0
        } 
    }
    if { $ma eq "Inf" && $mb eq "Inf" } {
        if { $sa == $sb } {
            return [list $sa "Inf" 0]
        } else {
            return [list 0 "Inf" 0]
        }
    }

    set comparison [compare [list $sa $ma $ea] [list $sb $mb $eb]] 

    if { $comparison == 1 } {
        return [Rescale $a]
    } elseif { $comparison == -1 } {
        return [Rescale $b]
    } elseif { $sb != $sa } {
        if { $sa } {
            return [Rescale $b]
        } else {
            return [Rescale $a]
        }
    } elseif { $sb && $eb > $ea } {
        # Both are negative and the same numerically. So return the one with the smallest exponent.
        return [Rescale $a]
    } elseif { $sb }  {
        # Negative with $eb < $ea now.
        return [Rescale $b]
    } elseif { $ea > $eb } {
        # Both are positive so return the one with the larger exponent
        return [Rescale $a]
    } else {
        return [Rescale $b]
    }
}

# maxmag -- max-magnitude
#     Return the larger of two numbers ignoring their signs.
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     larger of a or b ignoring their signs.
#
proc ::math::decimal::maxmag {a b} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}


    if { $ma eq "Inf" && $mb eq "Inf" } {
        if { $sa == 0 || $sb == 0 } {
            return [list 0 "Inf" 0]
        } else {
            return [list 1 "Inf" 0]
        }
    }

    set comparison [compare [list 0 $ma $ea] [list 0 $mb $eb]] 

    if { $comparison == 1 } {
        return [Rescale $a]
    } elseif { $comparison == -1 } {
        return [Rescale $b]
    } elseif { $sb != $sa } {
        if { $sa } {
            return [Rescale $b]
        } else {
            return [Rescale $a]
        }
    } elseif { $sb && $eb > $ea } {
        # Both are negative and the same numerically. So return the one with the smallest exponent.
        return [Rescale $a]
    } elseif { $sb }  {
        # Negative with $eb < $ea now.
        return [Rescale $b]
    } elseif { $ea > $eb } {
        # Both are positive so return the one with the larger exponent
        return [Rescale $a]
    } else {
        return [Rescale $b]
    }
}

# minmag -- min-magnitude
#     Return the smaller of two numbers ignoring their signs.
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     smaller  of a or b ignoring their signs.
#
proc ::math::decimal::minmag {a b} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $ma eq "Inf" && $mb eq "Inf" } {
        if { $sa == 1 || $sb == 1 } {
            return [list 1 "Inf" 0]
        } else { 
            return [list 0 "Inf" 0]
        }
    }

    set comparison [compare [list 0 $ma $ea] [list 0 $mb $eb]] 

    if { $comparison == 1 } {
        return [Rescale $b]
    } elseif { $comparison == -1 } {
        return [Rescale $a]
    } else {
        # They compared the same so now we use a normal comparison including the signs. This is per the specs.
        if { $sa > $sb } {
            return [Rescale $a]
        } elseif { $sb > $sa } {
            return [Rescale $b]
        } elseif { $sb && $eb > $ea } {
            # Both are negative and the same numerically. So return the one with the largest exponent.
            return [Rescale $b]
        } elseif { $sb }  {
            # Negative with $eb < $ea now.
            return [Rescale $a]
        } elseif { $ea > $eb } {
            return [Rescale $b]
        } else {
            return [Rescale $a]
        }
    }
}

# fma - fused-multiply-add
#     Takes three operands. Multiplies the first two and then adds the third.
#     Only one rounding (Rescaling) takes place at the end instead of after
#     both the multiplication and again after the addition.
#
# Arguments:
#     a          First operand
#     b          Second operand
#     c          Third operand
#
# Result:
#     (a*b)+c
#
proc ::math::decimal::fused-multiply-add {a b c} {
    return [fma $a $b $c]
}

proc ::math::decimal::fma {a b c} {
    return [+ $c [* $a $b 0]]
}

# multiply or *
#     Multiply two numbers
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     Product of both (rescaled)
#
proc ::math::decimal::multiply {a b {rescale 1}} {
    return [* $a $b $rescale]
}

proc ::math::decimal::* {a b {rescale 1}} {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}
    
    if { $ma eq "NaN" || $mb eq "NaN" } {
        return [list 0 "NaN" 0]
    } 

    set sr [expr {$sa^$sb}]

    if { $ma eq "Inf" || $mb eq "Inf" } {
        if { $ma == 0 || $mb == 0 } {
            return [list 0 "NaN" 0]
        } else {
            return [list $sr "Inf" 0]
        }
    }

    set mr [expr {$ma * $mb}]
    set er [expr {$ea + $eb}]


    if { $rescale } {
        return [Rescale [list $sr $mr $er]]
    } else {
        return [list $sr $mr $er]
    }
}

# divide or /
#     Divide two numbers
#
# Arguments:
#     a          First operand
#     b          Second operand
#
# Result:
#     Quotient of both (rescaled)
#
proc ::math::decimal::divide {a b {rescale 1}} {
    return [/ $a $b]
}

proc ::math::decimal::/ {a b {rescale 1}} {
    variable precision

    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $ma eq "NaN" || $mb eq "NaN" } {
        return [list 0 "NaN" 0]
    } 

    set sr [expr {$sa^$sb}]

    if { $ma eq "Inf" } {
        if { $mb ne "Inf"} {
            return [list $sr "Inf" 0]
        } else {
            return [list 0 "NaN" 0]
        }
    }

    if { $mb eq "Inf" } {
        if { $ma ne "Inf"} {
            return [list $sr 0 0]
        } else {
            return [list 0 "NaN" 0]
        }
    }

    if { $mb == 0 } {
        if { $ma == 0 } {
            return [list 0 "NaN" 0]
        } else {
            return [list $sr "Inf" 0]
        }
    }
    set adjust 0
    set mr 0


    if { $ma == 0 } {
        set er [expr {$ea - $eb}]
        return [list $sr 0 $er]
    }
    if { $ma < $mb } {
        while { $ma < $mb } {
            set ma [expr {$ma * 10}]
            incr adjust
        }
    } elseif { $ma >= $mb * 10 } {
        while { $ma >= [expr {$mb * 10}] } {
            set mb [expr {$mb * 10}]
            incr adjust -1
        }
    }

    while { 1 } {
        while { $mb <= $ma } {
            set ma [expr {$ma - $mb}]
            incr mr
        }
        if { ( $ma == 0 && $adjust >= 0 ) || [string length $mr] > $precision + 1 } {
            break
        } else {
            set ma [expr {$ma * 10}]
            set mr [expr {$mr * 10}]
            incr adjust
        }
    }
    
    set er [expr {$ea - ($eb + $adjust)}]

    if { $rescale } {
        return [Rescale [list $sr $mr $er]]
    } else {
        return [list $sr $mr $er]
    }
}

# divideint -- Divide integer
#     Divide a by b and return the integer part of the division.
#
#  Basically, if we send a and b to the divideint (which returns i) 
#  and remainder function (which returns r) then the following is true:
#      a = i*b + r
#
# Arguments:
#     a          First operand
#     b          Second operand
#
#
proc ::math::decimal::divideint { a b } {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}
    set sr [expr {$sa^$sb}]


        
    if { $sr == 1 } {
        set sign_string "-"
    } else {
        set sign_string ""
    }

    if { ($ma eq "NaN" || $mb eq "NaN") || ($ma == 0 && $mb == 0 ) } {
        return "NaN"
    }

    if { $ma eq "Inf" || $mb eq "Inf" } {
        if { $ma eq $mb } {
            return "NaN"
        } elseif { $mb eq "Inf" } {
            return "${sign_string}0"
        } else {
            return "${sign_string}Inf"
        }
    }

    if { $mb == 0 } {
        return "${sign_string}Inf"
    }
    if { $mb == "Inf" } {
        return "${sign_string}0"
    }
    set adjust [expr {abs($ea - $eb)}]
    if { $ea < $eb } {
        set a_adjust 0
        set b_adjust $adjust
    } elseif { $ea > $eb } {
        set b_adjust 0
        set a_adjust $adjust
    } else {
        set a_adjust 0
        set b_adjust 0
    }

    set integer [expr {($ma*10**$a_adjust)/($mb*10**$b_adjust)}]
    return $sign_string$integer
}

# remainder -- Remainder from integer division.
#     Divide a by b and return the remainder part of the division.
#
#  Basically, if we send a and b to the divideint (which returns i) 
#  and remainder function (which returns r) then the following is true:
#      a = i*b + r
#
# Arguments:
#     a          First operand
#     b          Second operand
#
#
proc ::math::decimal::remainder { a b } {
    foreach {sa ma ea} $a {break}
    foreach {sb mb eb} $b {break}

    if { $sa == 1 } {
        set sign_string "-"
    } else {
        set sign_string ""
    }

    if { ($ma eq "NaN" || $mb eq "NaN") || ($ma == 0 && $mb == 0 ) } {
        if { $mb eq "NaN" && $mb ne $ma } {
            if { $sb == 1 } {
                set sign_string "-"
            } else {
                set sign_string ""
            }
            return "${sign_string}NaN"
        } elseif { $ma eq "NaN" } {
            return "${sign_string}NaN"
        } else {
            return "NaN"
        }
    } elseif { $mb == 0 } {
        return "NaN"
    }

    if { $ma eq "Inf" || $mb eq "Inf" } {
        if { $ma eq $mb } {
            return "NaN"
        } elseif { $mb eq "Inf" } {
            return [tostr $a]
        } else {
            return "NaN"
        }
    }

    if { $mb == 0 } {
        return "${sign_string}Inf"
    }
    if { $mb == "Inf" } {
        return "${sign_string}0"
    }

    lset a 0 0
    lset b 0 0
    if { $mb == 0 } {
        return "${sign_string}Inf"
    }
    if { $mb == "Inf" } {
        return "${sign_string}0"
    }

    set adjust [expr {abs($ea - $eb)}]
    if { $ea < $eb } {
        set a_adjust 0
        set b_adjust $adjust
    } elseif { $ea > $eb } {
        set b_adjust 0
        set a_adjust $adjust
    } else {
        set a_adjust 0
        set b_adjust 0
    }

    set integer [expr {($ma*10**$a_adjust)/($mb*10**$b_adjust)}]

    set remainder [tostr [- $a [* [fromstr $integer] $b 0]]]
    return $sign_string$remainder
}


# abs --
#     Returns the Absolute Value of a number
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#
# Result:
#     Absolute value (as a list)
#
 proc ::math::decimal::abs {a} {
     lset a 0 0
     return [Rescale $a]
 }


# Rescale --
#     Rescale the number (using proper rounding)
#
# Arguments:
#     a Number in decimal format 
#
# Result:
#     Rescaled number 
#
proc ::math::decimal::Rescale { a } {



    variable precision
    variable rounding
    variable maxExponent
    variable minExponent
    variable tinyExponent

    foreach {sign mantisse exponent} $a {break}

    set man_length [string length $mantisse]

    set adjusted_exponent [expr {$exponent + ($man_length -1)}]

    if { $adjusted_exponent < $tinyExponent } {
        set mantisse [lindex [round_$rounding [list $sign $mantisse [expr {abs($tinyExponent) - abs($adjusted_exponent)}]] 0] 1]
        return [list $sign $mantisse $tinyExponent]
    } elseif { $adjusted_exponent > $maxExponent } {
        if { $mantisse  == 0 } {
            return [list $sign 0 $maxExponent]
        } else {
            switch -- $rounding {
                half_even -
                half_up { return [list $sign "Inf" 0] }
                down -
                05up { 
                    return [list $sign [string repeat 9 $precision] $maxExponent]
                }
                ceiling { 
                    if { $sign } {
                        return [list $sign [string repeat 9 $precision] $maxExponent]
                    } else {
                        return [list 0 "Inf" 0] 
                    }
                }
                floor {
                    if { !$sign } {
                        return [list $sign [string repeat 9 $precision] $maxExponent]
                    } else {
                        return [list 1 "Inf" 0] 
                    }
                }
                default { }
            }
        }
    }

    if { $man_length <= $precision } {
        return [list $sign $mantisse $exponent]
    } 

    set  mantisse [lindex [round_$rounding [list $sign $mantisse [expr {$precision - $man_length}]] 0] 1]
    set exponent [expr {$exponent + ($man_length - $precision)}]
    
    # it is possible now that our rounding gave us a new digit in our mantisse
    # example rounding 999.9 to 1 digits  with precision 3 will give us 
    # 1000 back. 
    # This can only happen by adding a zero on the end of our mantisse however.
    # So we just chomp it off.

    set man_length_now [string length $mantisse]
    if { $man_length_now > $precision } {
        set mantisse [string range $mantisse 0 end-1]
        incr exponent
        # Check again to see if we have overflowed 
        # we change our test to >= because we have incremented exponent.
        if { $adjusted_exponent >= $maxExponent } {
            switch -- $rounding {
                half_even -
                half_up { return [list $sign "Inf" 0] }
                down -
                05up { 
                    return [list $sign [string repeat 9 $precision] $maxExponent]
                }
                ceiling { 
                    if { $sign } {
                        return [list $sign [string repeat 9 $precision] $maxExponent]
                    } else {
                        return [list 0 "Inf" 0] 
                    }
                }
                floor {
                    if { !$sign } {
                        return [list $sign [string repeat 9 $precision] $maxExponent]
                    } else {
                        return [list 1 "Inf" 0] 
                    }
                }
                default { }
            }
        }
    }
    return [list $sign $mantisse $exponent]
}

# tostr --
#     Convert number to string using appropriate method depending on extended 
#     attribute setting.
#
# Arguments:
#     number     Number to be converted
#
# Result:
#     Number in the form of a string
#
proc ::math::decimal::tostr { number } {
    variable extended
    switch -- $extended {
        0 { return [tostr_numeric $number] }
        1 { return [tostr_scientific $number] }
    }
}

# tostr_scientific --
#     Convert number to string using scientific notation as called for in 
#     Decmath specifications.
#
# Arguments:
#     number     Number to be converted
#
# Result:
#     Number in the form of a string
#
proc ::math::decimal::tostr_scientific {number} {
    foreach {sign mantisse exponent} $number {break}

    if { $sign } {
        set sign_string "-"
    } else {
        set sign_string ""
    }

    if { $mantisse eq "NaN" } {
        return "NaN"
    } 
    if { $mantisse eq "Inf" } {
        return ${sign_string}${mantisse}
    }


    set digits [string length $mantisse]
    set adjusted_exponent [expr {$exponent + $digits - 1}]

    # Why -6? Go read the specs on the website mentioned in the header. 
    # They choose it, I'm using it. They actually list some good reasons though.
    if { $exponent <= 0 && $adjusted_exponent >= -6 } {
        if { $exponent == 0 } {
            set string $mantisse
        } else {
            set exponent [expr {abs($exponent)}]
            if { $digits > $exponent } {
                set string [string range $mantisse 0 [expr {$digits-$exponent-1}]].[string range $mantisse [expr {$digits-$exponent}] end]            
                set exponent [expr {-$exponent}]            
            } else {
                set string 0.[string repeat 0 [expr {$exponent-$digits}]]$mantisse
            }
        }
    } elseif { $exponent <= 0 && $adjusted_exponent < -6 } {
        if { $digits > 1 } {

            set string [string range $mantisse 0 0].[string range $mantisse 1 end]            

            set exponent [expr {$exponent + $digits - 1}]            
            set string "${string}E${exponent}"
        }  else {
            set string "${mantisse}E${exponent}"
        }
    } else {
        if { $adjusted_exponent >= 0 } {
            set adjusted_exponent "+$adjusted_exponent"
        }
        if { $digits > 1 } {
            set string "[string range $mantisse 0 0].[string range $mantisse 1 end]E$adjusted_exponent"
        } else {
            set string "${mantisse}E$adjusted_exponent"
        }
    } 
    return $sign_string$string
}

# tostr_numeric --
#     Convert number to string using the simplified number set conversion
#     from the X3.274 subset of Decimal Arithmetic specifications.
#
# Arguments:
#     number     Number to be converted
#
# Result:
#     Number in the form of a string
#
proc ::math::decimal::tostr_numeric {number} {
    variable precision
    foreach {sign mantisse exponent} $number {break}

    if { $sign } {
        set sign_string "-"
    } else {
        set sign_string ""
    }

    if { $mantisse eq "NaN" } {
        return "NaN"
    } 
    if { $mantisse eq "Inf" } {
        return ${sign_string}${mantisse}
    }

    set digits [string length $mantisse]
    set adjusted_exponent [expr {$exponent + $digits - 1}]
    
    if { $mantisse == 0 } {
        set string 0
        set sign_string "" 
    } elseif { $exponent <= 0 && $adjusted_exponent >= -6 } {
        if { $exponent == 0 } {
            set string $mantisse
        } else {
            set exponent [expr {abs($exponent)}]
            if { $digits > $exponent } {
                set string [string range $mantisse 0 [expr {$digits-$exponent-1}]]
                set decimal_part [string range $mantisse [expr {$digits-$exponent}] end]
                set string ${string}.${decimal_part}
                set exponent [expr {-$exponent}]            
            } else {
                set string 0.[string repeat 0 [expr {$exponent-$digits}]]$mantisse
            }
        }
    } elseif { $exponent <= 0 && $adjusted_exponent < -6 } {
        if { $digits > 1 } {
            set string [string range $mantisse 0 0].[string range $mantisse 1 end]            
            set exponent [expr {$exponent + $digits - 1}]            
            set string "${string}E${exponent}"
        }  else {
            set string "${mantisse}E${exponent}"
        }
    } else {
        if { $adjusted_exponent >= 0 } {
            set adjusted_exponent "+$adjusted_exponent"
        }
        if { $digits > 1 && $adjusted_exponent >= $precision } {
            set string "[string range $mantisse 0 0].[string range $mantisse 1 end]E$adjusted_exponent"
        } elseif { $digits + $exponent <= $precision } {
            set string ${mantisse}[string repeat 0 [expr {$exponent}]]
        } else {
            set string "${mantisse}E$adjusted_exponent"
        }
    } 
    return $sign_string$string
}

# fromstr --
#     Convert string to number
#
# Arguments:
#     string      String to be converted
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::fromstr {string} {
    variable extended

    set string [string trim $string "'\""]

    if { [string range $string 0 0] == "-" } {
        set sign 1
        set string [string trimleft $string -]
        incr pos -1
    } else  {
        set sign 0
    }

    if { $string eq "Inf" || $string eq "NaN" } {
        if {!$extended} {
            # we don't allow these strings in the subset arithmetic.
            # throw error.
            error "Infinities and NaN's not allowed in simplified decimal arithmetic"
        } else {
            return [list $sign $string 0]
        }
    }

    set string [string trimleft $string "+-"]
    set echeck [string first "E" [string toupper $string]]
    set epart 0
    if { $echeck >= 0 } {
        set epart [string range $string [expr {$echeck+1}] end]
        set string [string range $string 0 [expr {$echeck -1}]]
    }

    set pos [string first . $string]    
    
    if { $pos < 0 } {
        if { $string == 0 } {
            set mantisse 0
            if { !$extended } {
                set sign 0
            }
        } else {
            set mantisse $string
        }
        set exponent 0
    } else {
        if { $string == "" } {
            return [list 0 0 0]
        } else {
            #stripping the leading zeros here is required to avoid some octal issues.  
            #However, it causes us to fail some tests with numbers like 0.00 and 0.0 
            #which test differently but we can't deal with now.
            set mantisse [string trimleft [string map {. ""} $string] 0]
            if { $mantisse == "" } {
                set mantisse 0
                if {!$extended} {
                    set sign 0
                }
            }
            set fraction [string range $string [expr {$pos+1}] end]
            set exponent [expr {-[string length $fraction]}]
        }
    }
    set exponent [expr {$exponent + $epart}]

    if { $extended } {
        return [list $sign $mantisse $exponent]
    } else {
        return [Rescale [list $sign $mantisse $exponent]]
    }
}

# ipart --
#     Return the integer part of a Decimal Number
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     
#
# Result:
#     Integer
#
proc ::math::decimal::ipart { a } {

    foreach {sa ma ea} $a {break}
 
    if { $ea == 0 } {
        if { $sa } {
            return -$ma
        } else {
            return $ma
        }
    } elseif { $ea > 0 } {
        if { $sa } {
            return [expr {-1 * $ma * 10**$ea}]
        } else {
            return [expr {$ma * 10**$ea}]
        }
    } else {
        if { [string length $ma] <= abs($ea) } {
            return 0
        } else {
            if { $sa } {
                set string_sign "-"
            } else {
                set string_sign ""
            }
            set ea [expr {abs($ea)}]
            return "${string_sign}[string range $ma 0 end-$ea]"
        }
    }
}

# round_05_up --
#     Round zero or five away from 0.
#     The same as round-up, except that rounding up only occurs 
#     if the digit to be rounded up is 0 or 5, and after overflow 
#     the result is the same as for round-down.
#
#     Bias: away from zero
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_05up {a digits} {
    foreach {sa ma ea} $a {break}

    if { -$ea== $digits } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr { $ma * 10**($digits+$ea) }]
    } else {
        set round_exponent [expr {$digits + $ea}]
        if { [string length $ma] <= $round_exponent } {
            if { $ma != 0 } {
                set mantissa 1
            } else {
                set mantissa 0
            }
            set exponent 0
        } else {
            set integer_part [ipart [list 0 $ma $round_exponent]]

            if { [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}0 -1]] == 0 } {
                # We are rounding something with fractional part .0 
                set mantissa  $integer_part
            } elseif { [string index $integer_part end] eq 0 || [string index $integer_part end] eq 5 } {
                set mantissa [expr {$integer_part + 1}]
            } else {
                set mantissa  $integer_part
            }
            set exponent [expr {-1 * $digits}]
        }
    }
    return [list $sa $mantissa $exponent]
}

# round_half_up --
#     
#     Round to the nearest. If equidistant, round up.
#     
#
#     Bias: away from zero
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_half_up {a digits} {
    foreach {sa ma ea} $a {break}
    
    if { $digits + $ea == 0 } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr {$ma *10 **($digits+$ea)}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer_part [ipart [list 0 $ma $round_exponent]]

        switch -- [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}5 -1]] {
            0 { 
                # We are rounding something with fractional part .5
                set mantissa [expr {$integer_part + 1}]
            }
            -1 { 
                set mantissa $integer_part
            }
            1 {
                set mantissa [expr {$integer_part + 1}]
            }
            
        }
    }
    set exponent [expr {-1 * $digits}]
    return [list $sa $mantissa $exponent]
}

# round_half_even --
#     Round to the nearest. If equidistant, round so the final digit is even.
#     Bias: none
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_half_even {a digits} {

    foreach {sa ma ea} $a {break}

    if { $digits + $ea == 0 } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr {$ma * 10**($digits+$ea)}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer_part [ipart [list 0 $ma $round_exponent]]

        switch -- [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}5 -1]] {
            0 { 
                # We are rounding something with fractional part .5
                if { $integer_part % 2 } {
                    # We are odd so round up
                    set mantissa [expr {$integer_part + 1}]
                } else {
                    # We are even so round down
                    set mantissa $integer_part
                }
            }
            -1 { 
                set mantissa $integer_part
            }
            1 {
                set mantissa [expr {$integer_part + 1}]
            }
        }
    }
    set exponent [expr {-1 * $digits}]
    return [list $sa $mantissa $exponent]
}

# round_half_down --
#
#     Round to the nearest. If equidistant, round down.
#     
#     Bias: towards zero
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_half_down {a digits} {
    foreach {sa ma ea} $a {break}
    
    if { $digits + $ea == 0 } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr {$ma * 10**($digits+$ea)}]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set integer_part [ipart [list 0 $ma $round_exponent]]
        switch -- [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}5 -1]] {
            0 {
                # We are rounding something with fractional part .5
                # The rule is to round half down.
                set mantissa $integer_part
            } 
            -1 {
                set mantissa $integer_part
            }
            1 {
                set mantissa [expr {$integer_part + 1}]
            }
        }
    }
    set exponent [expr {-1 * $digits}]
    return [list $sa $mantissa $exponent]
}

# round_down --
#
#     Round toward 0.  (Truncate)
#     
#     Bias: towards zero
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_down {a digits} {
    foreach {sa ma ea} $a {break}


    if { -$ea== $digits } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr { $ma * 10**($digits+$ea) }]
    } else {
        set round_exponent [expr {$digits + $ea}]
        set mantissa [ipart [list 0 $ma $round_exponent]]
    }

    set exponent [expr {-1 * $digits}]
    return [list $sa $mantissa $exponent]
}

# round_floor --
#
#     Round toward -Infinity.
#     
#     Bias: down toward -Inf.
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_floor {a digits} {
    foreach {sa ma ea} $a {break}

    if { -$ea== $digits } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr { $ma * 10**($digits+$ea) }]
    } else {
        set round_exponent [expr {$digits + $ea}]
        if { $ma == 0 } {
            set mantissa 0
        } elseif { !$sa } {
            set mantissa [ipart [list 0 $ma $round_exponent]]
        } else {
            set mantissa [expr {[ipart [list 0 $ma $round_exponent]] + 1}]
        }
    }
    set exponent [expr {-1 * $digits}]
    return [list $sa $mantissa $exponent]
}        

# round_up --
#
#     Round away from 0
#     
#     Bias: away from 0
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_up {a digits} {
    foreach {sa ma ea} $a {break}


    if { -$ea== $digits } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr { $ma * 10**($digits+$ea) }]
    } else {
        set round_exponent [expr {$digits + $ea}]
        if { [string length $ma] <= $round_exponent } {
            if { $ma != 0 } {
                set mantissa 1
            } else {
                set mantissa 0
            }
            set exponent 0
        } else {
            set integer_part [ipart [list 0 $ma $round_exponent]]
            switch -- [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}0 -1]] {
                0 {
                    # We are rounding something with fractional part .0
                    set mantissa $integer_part
                } 
                default {
                    set mantissa [expr {$integer_part + 1}]
                }
            }
            set exponent [expr {-1 * $digits}]
        }
    }
    return [list $sa $mantissa $exponent]
}

# round_ceiling --
#
#     Round toward Infinity
#     
#     Bias: up toward Inf.
#
# Arguments:
#     Number in the form of {sign mantisse exponent}
#     Number of decimal points to round to.
#
# Result:
#     Number in the form of {sign mantisse exponent}
#
proc ::math::decimal::round_ceiling {a digits} {
    foreach {sa ma ea} $a {break}
    if { -$ea== $digits } {
        return $a
    } elseif { $digits + $ea > 0 } {
        set mantissa [expr { $ma * 10**($digits+$ea) }]
    } else {
        set round_exponent [expr {$digits + $ea}]
        if { [string length $ma] <= $round_exponent } {
            if { $ma != 0 } {
                set mantissa 1
            } else {
                set mantissa 0
            }
            set exponent 0
        } else {
            set integer_part [ipart [list 0 $ma $round_exponent]]
            switch -- [compare [list 0 $ma $round_exponent] [list 0 ${integer_part}0 -1]] {
                0 {
                    # We are rounding something with fractional part .0
                    set mantissa $integer_part
                } 
                default {
                    if { $sa } {
                        set mantissa [expr {$integer_part}]
                    } else {
                        set mantissa [expr {$integer_part + 1}]
                    }
                }
            }
            set exponent [expr {-1 * $digits}]
        }
    }
    return [list $sa $mantissa $exponent]
}        

# is-finite 
#
#     Takes one operand and returns: 1 if neither Inf or Nan otherwise 0.
#     
#
# Arguments:
#     a - decimal number
#
# Returns:
#    
proc ::math::decimal::is-finite { a } {
    set mantissa [lindex $a 1]
    if { $mantissa == "Inf" || $mantissa == "NaN" } {
        return 0
    } else {
        return 1
    }
}

# is-infinite 
#
#     Takes one operand and returns: 1 if Inf otherwise 0.
#     
#
# Arguments:
#     a - decimal number
#
# Returns:
#    
proc ::math::decimal::is-infinite { a } {
    set mantissa [lindex $a 1]
    if { $mantissa == "Inf" } {
        return 1
    } else {
        return 0
    }
}

# is-NaN
#
#     Takes one operand and returns: 1 if NaN otherwise 0.
#     
#
# Arguments:
#     a - decimal number
#
# Returns:
#    
proc ::math::decimal::is-NaN { a } {
    set mantissa [lindex $a 1]
    if { $mantissa == "NaN" } {
        return 1
    } else {
        return 0
    }
}

# is-signed
#
#     Takes one operand and returns: 1 if sign is 1 (negative).
#     
#
# Arguments:
#     a - decimal number
#
# Returns:
#    
proc ::math::decimal::is-signed { a } {
    set sign [lindex $a 0]
    if { $sign } {
        return 1
    } else {
        return 0
    }
}

# is-zero
#
#     Takes one operand and returns: 1 if operand is zero otherwise 0.
#     
#
# Arguments:
#     a - decimal number
#
# Returns:
#    
proc ::math::decimal::is-zero { a } {
    set mantisse [lindex $a 1]
    if { $mantisse == 0 } {
        return 1
    } else {
        return 0
    }
}

Test Program

Here is the tcl program to run the test suites. Run the test program with a single option, which is the location of the test suite to run. As noted in the source, the tests relying on NaN's and some minor precision issues (very small zeros with fixed maxExponents smaller than the results) may fail. You will need to download the tests from http://speleotrove.com/decimal/dectest.html and http://speleotrove.com/decimal/dectest0.html

#!/usr/bin/tclsh
#
# Copyright 2011 Mark Alston. All rights reserved.
#
# Redistribution and use in source and binary forms, with or 
# without modification, are permitted provided that the following 
# conditions are met:
#
#   1. Redistributions of source code must retain the above copyright 
#      notice, this list of conditions and the following disclaimer.
#
#   2. Redistributions in binary form must reproduce the above copyright 
#      notice, this list of conditions and the following disclaimer in 
#      the documentation and/or other materials provided with the distribution.
#
#   THIS SOFTWARE IS PROVIDED BY Mark Alston ``AS IS'' AND ANY EXPRESS 
#   OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
#   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
#   ARE DISCLAIMED. IN NO EVENT SHALL Mark Alston OR CONTRIBUTORS 
#   BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
#   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
#   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
#   OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
#   WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 
#   OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
#   EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# decimal_test.tcl --
#
#     Test program for the Decimal package.
#
#     This package passes all tests in test suites:
#           http://speleotrove.com/decimal/dectest.html
#      and  http://speleotrove.com/decimal/dectest0.html
#
#      with the following exceptions:
#
#     The Decimal package fails some tests that require setting the max 
#     or min exponent to force truncation or rounding.  
#  
#     The Decimal package fails some tests which require the sign of 
#     zero to be set correctly during rounding     
#
#     The Decimal package cannot handle sNaNs 
#
#     If you find errors in this code please let me know at 
#         mark at beernut dot com
#
#
#
proc mcsplit "str splitStr {mc {\x00}}" {
    return [split [string map [list $splitStr $mc] $str] $mc]
}
source decimal.tcl
namespace import ::math::decimal::*

set filename $argv 

set file [open $filename r]
set max_command_length 0

set tests 1
set failed 0
set success 0
set skipped 0
foreach line [split [read $file] "\n"] {
    
    if { $line == "" } {
        continue 
    } elseif { [string range [string trim $line] 0 1] eq "--"} {
        #Comment line. Print it out for reference in the output.
        puts "$line"
    } elseif { [string first ":" $line]  >0 } { 
        #We are setting a variable
        # Clear out the comments in the line and get the variable and setting.
        set variable_line [split [lindex [mcsplit $line "--"] 0] ":"]
        set variable [string trim [lindex $variable_line 0]]
        set setting [string trim [lindex $variable_line 1]]
        puts "===== set $variable = $setting"
        setVariable $variable $setting
    } else {
        # We are perfoming a test
        regsub -all " +" $line " " line 
        set command_list [mcsplit $line "->"]

        set command_string [split [string trim [lindex $command_list 0]]]

        set test [lindex $command_string 0]
        set operation [lindex $command_string 1]

        set operand1 [lindex $command_string 2]


        if {[llength $command_string] > 3 } {
            set operand2 [lindex $command_string 3]
        } else {
            set operand2 ""
        }
        if {[llength $command_string] > 4 } {
            set operand3 [lindex $command_string 4]
        } else {
            set operand3 ""
        }
        

        if { [string first "?" $line] != -1 } {
            puts "Skipping test $test: Invalid operations not fully implemented"
            incr skipped
            continue
        }


        # We are going to skip NaN tests for now. 
        # We actually pass all the regular nan tests but not the 
        # sNaNs and propigating NaNs.
        #
        # I am not sure how to implement them or if anyone cares.
        
        if { [regexp {\-*sNaN} $operand1] || [regexp {\-*sNaN} $operand2] } {
            puts "Skipping test $test: sNaNs not implemented"
            incr skipped
            continue
        }

        if { [regexp {NaN[0-9].*} $operand1] || [regexp {NaN[0-9].*} $operand2] } {
            puts "Skipping test $test: propigating NaNs not implemented"
            incr skipped
            continue
        }

        if { $operand1 eq "\#" || $operand2 eq "\#" } {
            puts "Skipping test $test: nulls not implemented"
            incr skipped
            continue
        }

        set result_string [split [string trim [lindex $command_list 1]]]
        set desired_result [string trim [lindex $result_string 0] "'"]



        set a [fromstr $operand1]
        if { $operand2 != "" } {
            set b [fromstr $operand2]
        }

        if { $operand3 != "" } {
            set c [fromstr $operand3]
        }

        #skip operations with numbers with abs(exponent) > 999 above tends to lockup expr
        
        if { [expr {abs([lindex $a 2])}] > 999  } {
            puts "Skipping test $test: Exponent greater than 999"
            incr skipped
            continue
        }

        if { $operand2 != "" } {
            if { [expr {abs([lindex $b 2])}] > 999 } {
                puts "Skipping test $test: Exponent greater than 999"
                incr skipped
                continue
            }
        }

        if { $operand3 != "" } {
            if { [expr {abs([lindex $c 2])}] > 999 } {
                puts "Skipping test $test: Exponent greater than 999"
                incr skipped
                continue
            }
        }
        puts $line
        puts -nonewline "test: $test"        
        switch -- $operation {
            add -
            subtract {
                if { $operand2 == "" } {
                    set result [tostr [$operation $a]]
                } else {
                    set result [tostr [$operation $a $b]]
                }
            }
            min -
            minmag -
            max -
            maxmag -
            copysign {
                set result [tostr [$operation $a $b]]
            }
            fma {
                set result [tostr [$operation $a $b $c]]
            }
            divide -
            multiply {
                set result [tostr [$operation $a $b]]
            }
            compare {
                set result [$operation $a $b]
            }
            is-infinite -
            is-NaN -
            is-zero -
            is-signed - 
            copynegate {
                set result [tostr [$operation $a]]
            }
            divideint -
            remainder {
                set result [$operation $a $b]
            }
            plus -
            minus -
            abs {
                set result [tostr [$operation $a]]
            }
            default {
                puts "Skipping test $test: $operation not implemented"
                incr skipped
                continue
            }                
        }
        if { $result eq "Inf" } {
            set result "Infinity"
        } elseif { $result eq "-Inf" } {
            set result "-Infinity"
        }
        if { $result ne $desired_result } {
            puts "  FAILED $operation $operand1 $operand2 $operand3 $result != $desired_result"
            puts "precision  : [getVariable precision]"
            puts "maxExponent: [getVariable maxExponent]"
            puts "minExponent: [getVariable minExponent]"
            puts "rounding   : [getVariable rounding]"
            incr failed
        } else {
            puts "  Passed $operation $operand1 $operand2 $operand3 $result == $desired_result"
            incr success
        }
        incr tests
    }
}
close $file
puts "Passed $success tests."
puts "Failed $failed tests."
puts "Skipped $skipped tests."

arjen (3 may 2011) Simplified the string trim commands in [fromstr]. This caused some confusion in the syntax highlighting and only one command suffices.

beernutmark 2011-05-03 Thanks. I just cleaned up a bunch of if's to remove the expr. I also changed max, min, maxmag, minmag which now pass 100% of the tests. I made some minor other changes to tostr and fromstr as well to help pass a few more tests. I also changed rescale to now directly take a decimal number since it really didn't fit with the rest of the proc's by taking the separate list elements.

beernutmark 2011-05-05 Changed the namespace to ::decimal::math to help future integration with tcllib. Fixed namespace on a few procs.)