Arjen Markus (10 august 2004) "Real" numbers on computers can cause a lot of agony because they are implemented (most of the time) using binary rather than decimal arithmetic. Thus "obvious" computations like:
0.1 + 0.1 == 0.2
fail on most computer systems. There are many ways to solve these problems, though they almost always mean that you need another package or library than the default one - for instance arbitrary precision arithmetic packages.
An important class of applications where these problems are really annoying is that of financial computations. I am no expert in this field, but here is my thought: Why not use fixed-point arithmetic ?
The range for the numbers one deals with is fairly limited, the precision is too, but you do need "decimally" predictable results.
I have an implementation of the basic arithmetical operations in mind that should be easy to use. It is Tcl only, will not imply too much overhead as most will be done using integer computations and ought to be flexible enough for most such applications. I have not coded it yet, just thought it over and I want to use this Wiki page to see if there is any interest in it ...
AM (19 august 2004) I am slowly making progress with this little package. I do have a new addition procedure, but I want to make things more modular to save code and avoid mistakes. That is on its way ...
In answer to the question Why not use fixed-point arithmetic ?, why not use rational numbers, where the numerator and denominator are separately maintained? I believe some Scheme implementations can do that -- escargo 19 Aug 2004 (I was always amazed how close the rational approximation for pi, 355/113, came to the correct value.)
Lars H: Rational arithmetic suffers from the "combinatorial explosion"; the number of digits in the representation is typically doubled with every operation. Unless you want to compute an exact value, you probably don't want that. (As for pi having good rational approximations, there is another side also of that coin. The common proofs that e and pi are transcendental numbers are based on a lemma that non-rational algebraic numbers cannot be that easily approximated using rational numbers.)
AM I agree with Lars here. Furthermore: being able to represent 1/3 exactly does not solve the problem of how to represent it in a decimal way .... It is easier to see the magnitude of 0.33 than of 1740/5742.
escargo 23 Aug 2004 - Good answers. I have used fixed-point arithmetic on processors that had integer multiply/divide but no floating point. One of the problems using it was that the binary point was programmer maintained. You were OK doing addition/subtraction, but when doing multiplication and division, you had to normalize results yourself. There are times when doing decimal arithmetic in BCD really does seem like the best way. (While I think rational numbers have their place, I was interested in arguments for and against them, rather than being solely for them.)
KBK 2005-12-23 Hmm, how'd I miss this page the first time around? Anyway, one comment about rationals and the combinatorial explosion; it is real. The simplest example is (1 + sqrt(5))/2 = 1.618. The convergents for the continued fraction representation of this number are ratios of successive Fibonacci numbers:
2/1, 3/2, 5/3, 8/5, 13/8, ...
and numerator and denominator grow exponentially. Any quadratic surd that is not an integer has a similar problem. (Incidentally, approximating surds by Newton-Raphson iteration yields precisely the same convergents.)
The common case for fixed point is in financial calculations, which have well-defined rounding conventions (and numbers of decimal digits to retain) in accounting practice.
Here is a second version of such a package ... (It surprised me how difficult it is to adequately insert the decimal point ;)
(More modular in set-up than the first one, but still it lacks proper rounding - I want to get a look at the source code on paper before moving on to tackle that problem ;) -- AM 6 september 2004
# fixedpoint.tcl --
# Implement fixed-point arithmetic
#
# fixedpoint --
# Namespace for the procedures and variables
#
namespace eval ::math::fixedpoint {
variable std_precision 3
variable scale {1 10 100 1000 10000 100000 1000000
100000000 1000000000 10000000000}
variable scale_current
variable rounding_method
}
# setprecision --
# Set the standard precision for computations
#
# Arguments:
# newvalue Integer between 0 and 9, giving the number
# of decimals to work with (optional)
# Result:
# Returns the standard precision (also if no argument given)
# Side effect:
# The variable std_precision is set along with several others
#
proc ::math::fixedpoint::setprecision { {newvalue {}} } {
variable std_precision
variable scale
variable scale_current
if { $newvalue == {} } {
return $std_precision
}
if { $newvalue < 0 || $newvalue > 9 } {
return -code error "Number of decimals must be between 0 and 9"
}
set scale_current [expr {wide([lindex $scale $std_precision])}]
return $std_precision
}
# setrounding --
# Set the rounding method
#
# Arguments:
# newname Name of method (optional):
# round-up, round-to-zero, floor, ceil, round-to-even
# Result:
# Returns the rouding method (also if name not given)
# Side effect:
# The variable rounding_method is set
#
proc ::math::fixedpoint::setrounding { {newname {}} } {
variable rounding_method
if { $newname == {} } {
return $rounding_method
}
switch -- $newname {
"round-up" -
"round-to-zero" -
"floor" -
"ceil" -
"round-to-even" {
set rounding_method $newname
}
default {
return -code error "Rounding method must be one of:
round-up, round-to-zero, floor, ceil or round-to-even"
}
}
return $rounding_method
}
# fixed --
# Convert a numerical value (interpreted as string!) to a fixed-point
# number using the current precision
#
# Arguments:
# strval String representation of the value
# precision Precision to use (optional; defaults to standard)
#
# Result:
# A fixed-point number
#
proc ::math::fixedpoint::fixed { strval {precision {}}} {
variable std_precision
if { ![string is double -strict $strval] } {
return -code error "Argument must be a valid number"
}
if { $precision == {} } {
set precision $std_precision
}
#
# Find the decimal point (if any) and the exponent (if any)
#
set p [string first . $strval]
set e [string first e $strval]
set E [string first e $strval]
set newval $strval
set expon 0
if { $e != -1 } {
set newval [string range $strval 0 [expr {$e-1}]]
set expon [string range $strval [expr {$e+1}] end]
}
if { $E != -1 } {
set newval [string range $strval 0 [expr {$E-1}]]
set expon [string range $strval [expr {$E+1}] end]
}
if { $p != -1 } {
set whole [string range $newval 0 [expr {$p-1}]]
set fract [string range $newval [expr {$p+1}] end]
}
set expon [expr {$precision+$expon-[string length $fract]}]
while { $expon > 0 } {
set fract "${fract}0"
incr expon -1
}
if { $expon < 0 } {
set fract [string range $fract 0 end$expon]
# TODO: insert proper rounding ...
}
return [list [expr wide($whole$fract)] $precision]
}
# tostring --
# Convert a fixed-point number to a string
#
# Arguments:
# fixedval Fixed-point number
#
# Result:
# A string with the decimal representation
#
proc ::math::fixedpoint::tostring { fixedval } {
set m [lindex $fixedval 0]
if { $m == 0 } {
return "0."
}
set p [lindex $fixedval 1]
set sign ""
if { $m < 0 } {
set sign "-"
set m [expr {-$m}]
}
if { [string length $m] < $p } {
set m "[string repeat 0 [expr {$p-[string length $m]}]]$m"
set sign "${sign}0"
}
return $sign[string range $m 0 end-$p].[string range $m end-[expr {$p-1}] end]
}
# SplitNumber --
# Split the encoded number in an integral and a fractional part
#
# Arguments:
# fixed Fixed-point number
# int Name of variable to hold integral part
# fract Name of variable to hold fractional part
# scale_factor Name of variable to hold scale factor
#
# Result:
# int and fract are set
#
proc ::math::fixedpoint::SplitNumber { fixed int fract scale_factor} {
variable scale
upvar 1 $int Int
upvar 1 $fract Fract
upvar 1 $scale_factor ScaleFactor
set m [lindex $fixed 0]
set p [lindex $scale [lindex $fixed 1]]
set Int [expr {($m>=0)? $m/$p : -(-$m/$p)}]
set Fract [expr {($m>=0)?$m%$p:(-(-$m%$p))}]
set ScaleFactor $p
}
# RoundOff --
# Round off to the given precision
#
# Arguments:
# int Name of variable holding integral part
# fract Name of variable holding fractional part
# scale_factor Scale factor
#
# Result:
# int and fract are set
#
proc ::math::fixedpoint::RoundOff { int fract scale_factor } {
variable rounding_method
variable std_precision
variable scale
upvar 1 $int Int
upvar 1 $fract Fract
if { abs($Fract) > 0.5*$scale_factor } {
if { $Fract >= 0 } {
set Fract [expr {int($Fract/$scale_factor+0.5)}]
} else {
set Fract [expr {int($Fract/$scale_factor-0.5)}]
}
} else {
set Fract [expr {int($Fract/$scale_factor)}]
}
}
# + --
# Add two fixed-point numbers
#
# Arguments:
# op1 First operand
# op2 Second operand
#
# Result:
# The sum of the two
#
proc ::math::fixedpoint::+ { op1 {op2 {}} } {
variable std_precision
variable scale
variable scale_current
if { [llength $op2] == 0 } {
return $op1
}
# TODO: handle different precisions
# TODO: handle overflow
SplitNumber $op1 I1 F1 p1
SplitNumber $op2 I2 F2 p2
if { $p1 > $p2 } {
set F2 [expr {$F2*$p1/$p2}]
}
if { $p1 < $p2 } {
set F1 [expr {$F1*$p2/$p1}]
set p1 $p2
}
set I [expr {$I1+$I2}]
set F [expr {$F1+$F2}]
if { $p1 > $scale_current } {
set reduce [expr {$p1/$scale_current}]
RoundOff I F $reduce
# TODO: rounding
}
if { $p1 < $scale_current } {
set enlarge [expr {$scale_current/$p1}]
set F [expr {$F*$enlarge}]
}
set I [expr {wide($I)*$scale_current+$F}]
return [list $I $std_precision]
}
# - --
# Subtract two fixed-point numbers
#
# Arguments:
# op1 First operand
# op2 Second operand
#
# Result:
# The difference of the two
#
proc ::math::fixedpoint::- { op1 {op2 {}} } {
if { $op2 == {} } {
foreach {m p} $op1 {break}
return [list [expr {-$m}] $p]
} else {
foreach {m p} $op2 {break}
return [+ $op1 [list [expr {-$m}] $p]]
}
}
# * --
# Multiply two fixed-point numbers
#
# Arguments:
# op1 First operand
# op2 Second operand
#
# Result:
# The product of the two
#
proc ::math::fixedpoint::* { op1 op2 } {
variable std_precision
variable scale
variable scale_current
# TODO: handle different precisions
# TODO: handle overflow
SplitNumber $op1 I1 F1 p1
SplitNumber $op2 I2 F2 p2
if { $p1 > $p2 } {
set F2 [expr {$F2*$p1/$p2}]
}
if { $p1 < $p2 } {
set F1 [expr {$F1*$p2/$p1}]
set p1 $p2
}
set I [expr {$I1*$I2}]
set F [expr {$I1*$F2+$I2*$F1+($F1*$F2)/$p1}]
if { $p1 > $scale_current } {
set reduce [expr {$p1/$scale_current}]
RoundOff I F $reduce
# TODO: rounding
}
if { $p1 < $scale_current } {
set enlarge [expr {$scale_current/$p1}]
set F [expr {$F*$enlarge}]
}
set I [expr {wide($I)*$scale_current+$F}]
return [list $I $std_precision]
}
# /--
# Divide two fixed-point numbers
#
# Arguments:
# op1 First operand
# op2 Second operand
#
# Result:
# The quotient of the two
#
proc ::math::fixedpoint::/ { op1 op2 } {
variable std_precision
variable scale
variable scale_current
# TODO: handle different precisions
# TODO: handle overflow
SplitNumber $op1 I1 F1 p1
SplitNumber $op2 I2 F2 p2
if { $p1 > $p2 } {
set F2 [expr {$F2*$p1/$p2}]
}
if { $p1 < $p2 } {
set F1 [expr {$F1*$p2/$p1}]
set p1 $p2
}
set I [expr {($I1*$p1+$F1)/($I2*$p1+$F2)}]
set F [expr {($p1*($I1*$p1+$F1-$I*($I2*$p1+$F2)))/($I2*$p1+$F2)}]
if { $p1 > $scale_current } {
set reduce [expr {$p1/$scale_current}]
RoundOff I F $reduce
# TODO: rounding
}
if { $p1 < $scale_current } {
set enlarge [expr {$scale_current/$p1}]
set F [expr {$F*$enlarge}]
}
set I [expr {wide($I)*$scale_current+$F}]
return [list $I $std_precision]
}
# initialisation --
#
namespace eval ::math::fixedpoint {
setprecision 3
setrounding round-up
}
set tcl_precision 17
puts [::math::fixedpoint::fixed "1.01"]
puts [::math::fixedpoint::fixed "1.011e3"]
puts [::math::fixedpoint::fixed "1.0112222222e3"]
puts [::math::fixedpoint::tostring [::math::fixedpoint::fixed "1.0112222222e3"]]
set op1 [::math::fixedpoint::fixed "1.01"]
set op2 [::math::fixedpoint::fixed "1.021"]
puts [::math::fixedpoint::tostring [::math::fixedpoint::+ $op1 $op2]]
puts [expr {1.01+1.021}]
puts [::math::fixedpoint::tostring [::math::fixedpoint::- $op1 $op2]]
puts [expr {1.01-1.021}]
puts [::math::fixedpoint::tostring [::math::fixedpoint::+ $op1 {10210 2}]]
puts "Sums:"
puts "Fixed -- floating-point"
foreach {x1 x2} {1.01 1.021 1.1 1.2 2.3 4.5 1.99 2.01 20.3 0.001} {
set op1 [::math::fixedpoint::fixed $x1]
set op2 [::math::fixedpoint::fixed $x2]
puts "$x1+$x2: \
[::math::fixedpoint::tostring [::math::fixedpoint::+ $op1 $op2]] \
-- [expr {$x1+$x2}]"
}
puts "Differences:"
puts "Fixed -- floating-point"
foreach {x1 x2} {1.01 1.021 1.1 1.2 2.3 4.5 1.99 2.01 20.3 0.001} {
set op1 [::math::fixedpoint::fixed $x1]
set op2 [::math::fixedpoint::fixed $x2]
puts "$x1-$x2: \
[::math::fixedpoint::tostring [::math::fixedpoint::- $op1 $op2]] \
-- [expr {$x1-$x2}]"
}
puts "Products:"
puts "Fixed -- floating-point"
foreach {x1 x2} {1.01 1.021 1.1 1.2 2.3 4.5 1.99 2.01 20.3 0.001} {
set op1 [::math::fixedpoint::fixed $x1]
set op2 [::math::fixedpoint::fixed $x2]
puts "$x1*$x2: \
[::math::fixedpoint::tostring [::math::fixedpoint::* $op1 $op2]] \
-- [expr {$x1*$x2}]"
}
puts "Quotients:"
puts "Fixed -- floating-point"
foreach {x1 x2} {1.021 1.01 1.1 1.2 2.3 4.5 1.99 2.01 20.03 0.001} {
set op1 [::math::fixedpoint::fixed $x1]
set op2 [::math::fixedpoint::fixed $x2]
puts "$x1/$x2: \
[::math::fixedpoint::tostring [::math::fixedpoint::/ $op1 $op2]] \
-- [expr {$x1/$x2}]"
}
----
[AM] (20 december 2005) Here is a much simpler approach (it is also less complete, but I know how to proceed):
# decimal.tcl --
# Dealing with decimal arithmetic
#
# Note:
# - This is a limited implementation
# - After using "fixed-precision" all existing "decimal numbers"
# become invalid, so use it at initialisation only to control
# the precision (and the range)
#
# namespace --
# Define the DecimalArithmetic namespace
#
package require Tcl 8.4
namespace eval ::DecimalArithmetic {
variable _precision_
variable _maxvalue_
variable _tofrac_
namespace export fixed-precision + < > >= == <= < != \
tostring fromstring
}
# fixed-precision --
# Set the precision and compute the auxiliary parameters
# Arguments:
# prec Number of digits in a fraction
# Result:
# None
# Side effect:
# Auxiliary parameter set - see also the note at the start
#
proc ::DecimalArithmetic::fixed-precision {prec} {
variable _precision_
variable _maxvalue_
variable _tofrac_
set _precision_ $prec
set _maxvalue_ [expr {pow(10.0,18-$prec)}]
set _tofrac_ [expr "wide(1[string repeat 0 $prec])"]
}
# comparisons --
# Define the comparison procedures
#
namespace eval ::DecimalArithmetic {
foreach op {< > >= == <= < !=} {
proc $op {op1 op2} [string map [list OP $op] \
{expr {$op1 OP $op2}}]
}
}
# default precision --
# Define the comparison procedures
#
namespace eval ::DecimalArithmetic {
fixed-precision 8
}
# + --
# Add two decimal numbers
# Arguments:
# op1 First operand
# op2 Second operand
# Result:
# Sum of the two or Inf
# Note:
# Take care not to overflow the range!
#
proc ::DecimalArithmetic::+ {op1 op2} {
variable _maxvalue_
if { $op1 == "Inf" || $op2 == "Inf" } {
return "Inf"
}
if { abs(double($op1)+double($op2)) > $_maxvalue_ } {
return "Inf"
}
return [expr {$op1+$op2}]
}
# tostring --
# Convert a decimal number to a string ("1.001" for instance)
# Arguments:
# number Number to be converted
# Result:
# Decimal representation of the number
#
proc ::DecimalArithmetic::tostring {number} {
variable _precision_
if { $number == "Inf" } {
return "Inf"
}
return "[string range $number 0 end-$_precision_].[string range $number end-[expr {$_precision_-1}] end]"
}
# fromstring --
# Convert a string to a decimal number ("1.001" for instance)
# Arguments:
# string String to be converted
# Result:
# Decimal number equivalent to the string (in current precision)
#
proc ::DecimalArithmetic::fromstring {string} {
variable _precision_
variable _maxvalue_
variable _tofrac_
if { ! [string is double $string] } {
return "Inf"
}
foreach {int frac} [split $string "."] {break}
if { $int > $_maxvalue_ } {
return "Inf"
}
set frac [string range "$frac[string repeat 0 $_precision_]" 0 [expr {$_precision_-1}]]
set frac [string triml $frac 0]
if { $frac == "" } {
set frac "0"
}
set frac [expr {wide($frac)}]
if { $int >= 0 } {
return [expr {wide($int)*$_tofrac_+$frac}]
} else {
return [expr {wide($int)*$_tofrac_-$frac}]
}
}
# main --
# Test the procedures ...
#
namespace import ::DecimalArithmetic::*
set x 0.0
set dx 1.01
set y [fromstring 0.0]
set dy [fromstring 1.01]
set tcl_precision 17
for { set i 0 } { $i < 10 } { incr i } {
set x [expr {$x+$dx}]
set y [+ $y $dy]
puts "$x -- [tostring $y]"
}escargo 20 Dec 2005 - It is the normalizing of the decimal point (or binary point) after multiplication and division that are the true tests of fixed-point arithmetic. And fixed-point does not imply decimal; it can just as easily be a binary point as a decimal point.
AM Hm, the above strategy makes it possible to work with any radix ... Though I can not foresee the usefulness of a radix like 17, it would work ... The main reason for doing it like this is that it is much simpler than my first idea.
escargo - I guess using a DecimalArithmetic namespace made me think that nondecimal radix (or binary radix) was less supported. (I used to use some computers that used a radix 40 character set; it allowed for 26 letters, 10 digits, a space character and three other characters (",", "$", ".", if memory serves) in the character set. This let them pack 3 characters in a 16-bit word.)
Sarnold What if I want 3.2571% of 17000.00$ rounded to 2 decimals? Could the above scripts handle this?
AM It should be possible (rounding is not supported yet, but there is no essential limitation) - as long as you stay within the limits of the precision.
beernutmark 2011-05-02 - Decimal Arithmetic Package for tcl 8.5 handles this easily with 8 rounding modes to choose from.
puts "[::Decimal::tostr [::Decimal::round_half_even [::Decimal::/ [::Decimal::* [::Decimal::fromstr 3.2571] [::Decimal::fromstr 17000.00]] $::Decimal::onehundred] 2]]"
553.71