[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.

----
 # 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

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