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 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.
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.2 - 0.2} --> 8.0 instead of expr {8.2 - 0.2} --> 7.999999999999999
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.
package require Tcl 8.5
package provide 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 divideint remainder
#
# 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 ::Decimal {
variable precision 20
variable maxExponent 999
variable minExponent -998
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 setScale
}
# setVariable
# Set the desired variable
#
# Arguments:
# variable setting
#
# Result:
# None
#
proc ::Decimal::setVariable {variable setting} {
variable rounding
variable precision
variable extended
variable maxExponent
variable minExponent
switch -- $variable {
rounding {set rounding $setting}
precision {set precision $setting}
extended {set extended $setting}
maxExponent {set maxExponent $setting}
minExponent {set minExponent $setting}
default {}
}
}
# add or +
# Add two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# Sum of both (rescaled)
#
proc ::Decimal::add {a b {rescale 1}} {
return [+ $a $b $rescale]
}
proc ::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 $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
#
#
# Arguments:
# a operand
#
# Result:
# a with sign flipped
#
proc ::Decimal::negate { a } {
return [copy-negate $a]
}
proc ::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 ::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 ::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 ::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 ::Decimal::subtract {a {b {}} {rescale 1}} {
return [- $a $b]
}
proc ::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 ::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]]
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 ::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
} else {
return [list 1 0 0]
}
}
if { $ma eq "Inf" && $mb eq "Inf" } {
if { $sa == $sb } {
return [list $sa "Inf" 0]
} else {
return [list 1 "Inf" 0]
}
}
set comparison [- [list $sa $ma $ea] [list $sb $mb $eb]]
if { [lindex $comparison 0] && [lindex $comparison 1] != 0 } {
return $a
} elseif { [lindex $comparison 1] == 0 } {
if { $sa == $sb } {
if { $ea < $eb } {
return $a
} else {
return $b
}
} elseif { $sa < $sb } {
return $b
} else {
return $a
}
} else {
return $b
}
}
# max
# Return the larger of two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# larger of a or b
#
proc ::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
} else {
return [list 0 0 0]
}
}
if { $ma eq "Inf" && $mb eq "Inf" } {
if { $sa == $sb } {
return [list $sa "Inf" 0]
} else {
return [list 0 "Inf" 0]
}
}
set comparison [- [list $sa $ma $ea] [list $sb $mb $eb]]
if { [lindex $comparison 0] && [lindex $comparison 1] != 0 } {
return $b
} elseif { [lindex $comparison 1] == 0 } {
if { $sa == $sb } {
if { $ea < $eb } {
return $b
} else {
return $a
}
} elseif { $sa < $sb } {
return $a
} else {
return $b
}
} else {
return $a
}
}
# 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 ::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 $a
} elseif { $comparison == -1 } {
return $b
} else {
if { $sa > $sb } {
return $b
} else {
return $a
}
}
}
# 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 ::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 $b
} elseif { $comparison == -1 } {
return $a
} else {
if { $sa > $sb } {
return $a
} else {
return $b
}
}
}
# 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 ::Decimal::fma {a b c} {
return [+ $c [::Decimal::* $a $b 0]]
}
# multiply or *
# Multiply two numbers
#
# Arguments:
# a First operand
# b Second operand
#
# Result:
# Product of both (rescaled)
#
proc ::Decimal::multiply {a b {rescale 1}} {
return [::Decimal::* $a $b $rescale]
}
proc ::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 $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 ::Decimal::divide {a b {rescale 1}} {
return [/ $a $b]
}
proc ::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 >= [expr {$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] > [expr {$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 $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 ::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 ::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 ::Decimal::abs {a} {
lset a 0 0
return $a
}
# 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 ::Decimal::Rescale {sign mantisse exponent} {
variable precision
variable rounding
variable maxExponent
variable minExponent
set man_length [string length $mantisse]
set adjusted_exponent [expr {$exponent + $man_length}]
# puts "FAILED adjusted: $adjusted_exponent"
# puts "FAILED min: $minExponent"
# puts "FAILED max: $maxExponent"
if { $adjusted_exponent < $minExponent } {
return [list $sign 0 $minExponent]
} elseif { $adjusted_exponent > $maxExponent } {
return [list $sign "Inf" 0]
}
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
}
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 ::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 ::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}]
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 ::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 { [expr {$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 ::Decimal::fromstr {string} {
variable extended
set string [string trim $string "'"]
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 {
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 $sign $mantisse $exponent]
}
}
# ipart --
# Return the integer part of a Decimal Number
#
# Arguments:
# Number in the form of {sign mantisse exponent}
#
#
# Result:
# Integer
#
proc ::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] <= [expr {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 ::Decimal::round_05up {a digits} {
foreach {sa ma ea} $a {break}
if { -$ea== $digits } {
return $a
} elseif { [expr {$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 ::Decimal::round_half_up {a digits} {
foreach {sa ma ea} $a {break}
if {[expr {$digits + $ea}] == 0 } {
return $a
} elseif { [expr {$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 ::Decimal::round_half_even {a digits} {
foreach {sa ma ea} $a {break}
if {[expr {$digits + $ea}] == 0 } {
return $a
} elseif { [expr {$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 { [expr {$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 ::Decimal::round_half_down {a digits} {
foreach {sa ma ea} $a {break}
if {[expr {$digits + $ea}] == 0 } {
return $a
} elseif { [expr {$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 ::Decimal::round_down {a digits} {
foreach {sa ma ea} $a {break}
if { -$ea== $digits } {
return $a
} elseif { [expr {$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 ::Decimal::round_floor {a digits} {
foreach {sa ma ea} $a {break}
if { -$ea== $digits } {
return $a
} elseif { [expr {$digits + $ea}] > 0 } {
set mantissa [expr { $ma * 10**($digits+$ea) }]
} else {
set round_exponent [expr {$digits + $ea}]
if { !$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 ::Decimal::round_up {a digits} {
foreach {sa ma ea} $a {break}
if { -$ea== $digits } {
return $a
} elseif { [expr {$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 ::Decimal::round_ceiling {a digits} {
foreach {sa ma ea} $a {break}
if { -$ea== $digits } {
return $a
} elseif { [expr {$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 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 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 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 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 is-zero { a } {
set mantisse [lindex $a 1]
if { $mantisse == 0 } {
return 1
} else {
return 0
}
}
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.
#!/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 sNaN's
#
# 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
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"
::Decimal::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
# sNaN's and propigating NaN's.
#
# I am not sure how to implement them or if anyone cares.
if { $operand1 eq "sNaN" || $operand2 eq "sNaN" } {
puts "Skipping test $test: sNaN's 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 [::Decimal::fromstr $operand1]
if { $operand2 != "" } {
set b [::Decimal::fromstr $operand2]
}
if { $operand3 != "" } {
set c [::Decimal::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 -nonewline "test: $test"
switch -- $operation {
add -
subtract {
if { [string first "NaN" $line] != -1 } {
puts "Skipping test $test: NaN's not fully implemented"
incr skipped
continue
}
if { $operand2 == "" } {
set result [::Decimal::tostr [::Decimal::$operation $a]]
} else {
set result [::Decimal::tostr [::Decimal::$operation $a $b]]
}
}
min -
minmag -
max -
maxmag -
copysign {
if { [string first "NaN" $line] != -1 } {
puts "Skipping test $test: NaN's not fully implemented"
incr skipped
continue
}
set result [::Decimal::tostr [::Decimal::$operation $a $b]]
}
fma {
if { [string first "NaN" $line] != -1 } {
puts "Skipping test $test: NaN's not fully implemented"
incr skipped
continue
}
set result [::Decimal::tostr [::Decimal::$operation $a $b $c]]
}
divide -
multiply {
if { [string first "NaN" $line] != -1 } {
puts "Skipping test $test: NaN's not fully implemented"
incr skipped
continue
}
set result [::Decimal::tostr [::Decimal::$operation $a $b]]
}
compare {
set result [::Decimal::$operation $a $b]
}
is-infinite -
is-NaN -
is-zero -
is-signed -
copynegate {
set result [::Decimal::tostr [::Decimal::$operation $a]]
}
divideint -
remainder {
set result [::Decimal::$operation $a $b]
}
plus -
minus -
abs {
# We are going to skip sNaN tests for now.
# Not sure how to implement them.
if { [string first "NaN" $line] != -1 } {
puts "Skipping test $test: NaN's not fully implemented"
incr skipped
continue
}
set result [::Decimal::tostr [::Decimal::$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"
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."