[Arjen Markus] (9 july 2003) Having read a book about numerical accuracy and stability, as well as several technical descriptions about interval arithmetic (certain Fortran compilers support it), I wanted to have this facility in Tcl too. 

The idea: a number is not completely accurate, but there is some error margin. Then summing two such numbers becomes a matter of manipulating ''intervals'':

   x = [[a,b]]
   y = [[c,d]]

then:

   x+y = [[a+c,b+d]]

for instance.

I know that it has severe drawbacks (and for these I want to provide some counter measures, like an alternative technique, called "running error analysis"):

   * It assumes that any arithmetic operation works on independent intervals, which is not true for x*x for instance
   * As a consequence intervals will grow indefinitely, to the point that there is no practical information anymore

Well, creating the procs was easy and fast, so we can focus on obtaining experience :)
Here is the script with a simple example.

'''Note''': the functional or prefix style of arithmetic is really useful here - it will allow us to overload the procedures without changing the code everywhere.

----
 # interval.tcl --
 #    Rudimentary support for interval arithmetic in Tcl
 #
 # To do:
 #    Elementary functions - cos, sin, tan, abs, exp, log, pow
 #    General function - function
 #
 namespace eval ::math::interval {
    namespace export + - * / \
                     contains interval begin-end \
                     square abs
 }

 # + --
 #    Add two interval numbers
 #
 # Arguments:
 #    operand1     First operand (interval number or regular number)
 #    operand2     Second operand (interval number or regular number)
 # Result:
 #    Sum of the two intervals
 #
 proc ::math::interval::+ {operand1 operand2} {
    foreach {begin1 end1} [concat $operand1 $operand1] {break}
    foreach {begin2 end2} [concat $operand2 $operand2] {break}

    return [list [expr {$begin1+$begin2}] [expr {$end1+$end2}]]
 }

 # - --
 #    Subtract two interval numbers
 #
 # Arguments:
 #    operand1     First operand (interval number or regular number)
 #    operand2     Second operand (interval number or regular number)
 # Result:
 #    Difference of the two intervals
 #
 proc ::math::interval::- {operand1 operand2} {
    foreach {begin1 end1} [concat $operand1 $operand1] {break}
    foreach {begin2 end2} [concat $operand2 $operand2] {break}

    return [list [expr {$begin1-$begin2}] [expr {$end1-$end2}]]
 }

 # * --
 #    Multiply two interval numbers
 #
 # Arguments:
 #    operand1     First operand (interval number or regular number)
 #    operand2     Second operand (interval number or regular number)
 # Result:
 #    Difference of the two intervals
 #
 proc ::math::interval::* {operand1 operand2} {
    foreach {begin1 end1} [concat $operand1 $operand1] {break}
    foreach {begin2 end2} [concat $operand2 $operand2] {break}

    set v1  [expr {$begin1*$begin2}]
    set v2  [expr {$begin1*$end2  }]
    set v3  [expr {$end1  *$begin2}]
    set v4  [expr {$end1  *$end2  }]
    set min $v1
    set max $v1
    foreach v [list $v2 $v3 $v4] {
       set min [expr { $min<$v? $min : $v }]
       set max [expr { $max>$v? $max : $v }]
    }
    return [list $min $max]
 }

 # / --
 #    Divide two interval numbers
 #
 # Arguments:
 #    operand1     First operand (interval number or regular number)
 #    operand2     Second operand (interval number or regular number)
 # Result:
 #    Difference of the two intervals
 # Note:
 #    The second interval may include 0!
 #
 proc ::math::interval::/ {operand1 operand2} {
    foreach {begin1 end1} [concat $operand1 $operand1] {break}
    foreach {begin2 end2} [concat $operand2 $operand2] {break}

    if { $begin2 <= 0.0 && $end2 >= 0.0 } { error "Division by interval containing zero" }

    *  $operand1 [list [expr {1.0/$begin2}] [expr {1.0/$end2}]]
 }

 # contains --
 #    Check if the interval contains a given value
 #
 # Arguments:
 #    interval     Interval number in question
 #    value        Value that is or is not contained
 # Result:
 #    0 if not, 1 if it is contained
 #
 proc ::math::interval::contains {interval value} {
    foreach {begin end} [concat $interval $interval] {break}

    if { $begin <= $value && $end >= $value } {
       return 1
    } else {
       return 0
    }
 }

 # interval --
 #    Create an interval number from a centre and a tolerance
 #
 # Arguments:
 #    centre       Central number for the interval
 #    tolerance    Tolerance (relative error, defaults to 1.0e-10)
 # Result:
 #    Interval number
 # Note:
 #    For centre = 0.0, the tolerance is chosen to be absolute
 #    - a hack in fact.
 #
 proc ::math::interval::interval { centre {tolerance 1.0e-10} } {
    if { $centre != 0.0 } {
       set begin [expr {(1.0-$tolerance)*$centre}]
       set end   [expr {(1.0+$tolerance)*$centre}]
    } else {
       set begin [expr {-$tolerance}]
       set end   $tolerance
    }

    return [list $begin $end]
 }

 # begin-end --
 #    Create an interval number from two extremes
 #
 # Arguments:
 #    low          Lower bound for the interval
 #    high         Higher bound for the interval
 # Result:
 #    Interval number
 #
 proc ::math::interval::begin-end { low high } {
    if { $low < $high } {
       return [list $low $high]
    } else {
       return [list $high $low]
    }
 }

 # square --
 #    Return the interval containing the square of an interval
 #
 # Arguments:
 #    interval     Interval to be squared
 # Result:
 #    Encompassing interval
 # Note:
 #    The trick is that the square of a number can not be negative!
 #    And the operands are correlated so that a direct multiplication
 #    gives a very pessimistic answer
 #
 proc ::math::interval::square { interval } {
    foreach {begin end} [concat $interval $interval] {break}

    set min2 [expr {$begin*$begin}]
    set max2 [expr {$end*$end}]

    if { $begin < 0.0 && $end > 0.0 } {
       return [list 0.0 [expr {$min2>$max2? $min2 : $max2}]]
    } else {
       return [list [expr {$min2<$max2? $min2 : $max2}] \
                    [expr {$min2>$max2? $min2 : $max2}] ]
    }
 }

 # abs --
 #    Return the interval containing the absolute value of an interval
 #
 # Arguments:
 #    interval     Interval to be treated
 # Result:
 #    Encompassing interval
 #
 proc ::math::interval::abs { interval } {
    foreach {begin end} [concat $interval $interval] {break}

    set min [expr {abs($begin)}]
    set max [expr {abs($end)}]

    if { $begin < 0.0 && $end > 0.0 } {
       return [list 0.0 [expr {$min>$max? $min : $max}]]
    } else {
       return [list [expr {$min<$max? $min : $max}] \
                    [expr {$min>$max? $min : $max}] ]
    }
 }

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

 puts [+ {0.9 1.1} {2.1 2.2}]
 puts [+ {0.9 1.1} 1]
 puts [+ 1 1]
 puts [- {0.9 1.1} {2.1 2.2}]
 puts [* {0.9 1.1} {2.1 2.2}]
 puts [/ {0.9 1.1} {2.1 2.2}]

 puts "Square:"
 puts [square {1 2}]
 puts [square {-1 2}]
 puts [square {-1 -2}]

 puts "Abs:"
 puts [abs {1 2}]
 puts [abs {-1 2}]
 puts [abs {-1 -2}]

 #
 # Now a more practical example
 #
 puts "Integrating a parabolic function f(x) = alpha * x**2"
 puts "Assumes that all steps have independent uncertainties!"

 set alpha   [interval 1 0.001]
 set alphapl 1
 set dx    0.1
 set sum   0.0
 set sumpl 0.0

 for { set i 0 } { $i < 300 } { incr i } {
    set x     [* $i $dx]
    set sum   [+ $sum   [* $alpha   [* $x $x]]]
    set sumpl [+ $sumpl [* $alphapl [* $x $x]]]
 }

 puts "(Interval) estimate of integral: $sum"
 puts "(Plain) estimate of integral:    $sumpl"

----
The output:

 3.0 3.3
 1.9 2.1
 2 2
 -1.2 -1.1
 1.89 2.42
 0.409090909091 0.52380952381
 Square:
 1 4
 0.0 4
 1 4
 Abs:
 1 2
 0.0 2
 1 2
 Integrating a parabolic function f(x) = alpha * x**2
 Assumes that all steps have independent uncertainties!
 (Interval) estimate of integral: 89460.9495 89640.0505
 (Plain) estimate of integral:    89550.5 89550.5

----
[TV] I remember the idea from 1st years physics practicum: propagation of measurement errors. Maybe give the whole idea a standard deviation functional integral interpretation and we're in interesting physicist area..

A main point to think about which is very relevant in digital signal processing as opposed to electronics and other continuous equations based computations and approximation methods is the ''correlation'' measure of the error, also expressable in bits (the information unit), which comes from  rounding errors in the computations themselves, also in the inevitable limited accuracy computations of the intervals, so you make sure they are real outer bounds.

For instance 

 (bin) 4 % set a [expr 1.0/6]
 0.166666666667
 (bin) 5 % expr 1.0/$a
 5.99999999999

is simplistic example of what in senstive boundary computations (as in electrical engineering is well known rule in differential equations related problems: most problems arise in the boundary conditions) can easily make your 'certain' interval just wrong.

I just reread and found you mention correlation yourself. The accumulation of error does however not take all information away, the uncertainty interval becomes large, as does the product, their ratio doesn't change to dramatically I'd think.

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