[Arjen Markus] (10 july 2003) A simple alternative to ''a priori'' analysis of numerical algorithms and [interval arithmic] is that you keep track of numerical errors while doing the calculations. The method is described in detail by Nicholas Higham. 

The script below is a basic implementation of the technique (just two operations) to illustrate it. 

----
 # numerr.tcl --
 #    Running error analysis
 #
 # Bibliography:
 #    Nicholas J. Higham:
 #       Accuracy and stability of Numerical Algorithms
 #    SIAM, 1996, ISBN 0-89871-355-2
 #
 # The idea:
 #    Floating-point operations are not exact, but they
 #    can be modelled as:
 #       fl( x op y ) = (1+D) (x op y)
 #    where |D| <= u, u is the unit roundoff or machine precision
 #
 # The procedures here return:
 # - The result fl( x op y )
 # - The accumulated error estimate
 # in the form of a two-elements list
 #
 namespace eval ::math::numerr {
    namespace export + - * / \

    variable roundoff 1.0e-16
 }

 # + --
 #    Add two numbers and return the result with error estimate
 #
 # Arguments:
 #    operand1     First operand (number with error estimate or regular number)
 #    operand2     Second operand (ditto)
 # Result:
 #    Sum and error estimate
 #
 proc ::math::numerr::+ {operand1 operand2} {
    variable roundoff
    set add1 [lindex $operand1 0]
    set err1 [lindex $operand1 1]
    set add2 [lindex $operand2 0]
    set err2 [lindex $operand2 1]
    if { $err1 == "" } { set err1 0.0 }
    if { $err2 == "" } { set err2 0.0 }

    return [list [expr {$add1+$add2}] \
                 [expr {abs($add1+$add2)*$roundoff+$err1+$err2}]]
 }

 # * --
 #    Multiply two numbers and return the result with error estimate
 #
 # Arguments:
 #    operand1     First operand (number with error estimate or regular number)
 #    operand2     Second operand (ditto)
 # Result:
 #    Porduct and error estimate
 #
 proc ::math::numerr::* {operand1 operand2} {
    variable roundoff
    set mult1 [lindex $operand1 0]
    set err1  [lindex $operand1 1]
    set mult2 [lindex $operand2 0]
    set err2  [lindex $operand2 1]
    if { $err1 == "" } { set err1 0.0 }
    if { $err2 == "" } { set err2 0.0 }

    return [list [expr {$mult1*$mult2}] \
       [expr {abs($mult1*$mult2)*$roundoff + \
              abs($mult2)*$err1 + abs($mult1)*$err2}]]
 }

 # Test code
 #
 namespace import ::math::numerr::*

 set sum  0.0
 set val  100.0
 for { set i 0 } { $i < 10 } { incr i } {
    set sum [+ $sum $val]
    puts "$sum -- [expr {[lindex $sum 1]/[lindex $sum 0]}]"
    set val [* $val 1.1]
 }

----
When run, it produces this output:

 100.0 1e-14 -- 1e-16
 210.0 4.2e-14 -- 2e-16
 331.0 9.93e-14 -- 3e-16
 464.1 1.8564e-13 -- 4e-16
 610.51 3.05255e-13 -- 5e-16
 771.561 4.629366e-13 -- 6e-16
 948.7171 6.6410197e-13 -- 7e-16
 1143.58881 9.14871048e-13 -- 8e-16
 1357.947691 1.2221529219e-12 -- 9e-16
 1593.7424601 1.5937424601e-12 -- 1e-15

As you can see, the relative error increases with each iteration (the regularity is a bit surprising, but it seems to be correct :)

----
[[
[Category Mathematics]
]]