This page is under development. Comments are welcome, but please load any comments in the comments section at the bottom of the page. Please include your wiki MONIKER and date in your comment with the same courtesy that I will give you. Aside from your courtesy, your wiki MONIKER and date as a signature and minimal good faith of any internet post are the rules of this TCL-WIKI. Its very hard to reply reasonably without some background of the correspondent on his WIKI bio page. Thanks, gold 12Dec2018
gold Here is some TCL starter code for Kahan compensated summation algorithm and Neumaier variant summation algorithm.
The Kahan summation algorithm is used to counter round-off error with some effectiveness. Note: One won't see these (single precision) errors in normal TCL 8.6 use. The TCL script below uses <format %7.1f $inputs> to force single precision (sp) on expr and math operator statements. Expr from TCL 8.6 uses double precision inside its calculations and throws round off errors beyond normal sp use. However, this example of the kahan summation algorithm seems to be working correctly and results can be compared to the gimmicked summation procs below.
Kahan, Neumaier, Babuska, and others developed balanced or compensated summation algorithms. Some balanced algorithms subtract the compensations from the summing terms as the terms accumulate. Other balanced algorithms model total compensation with a first, second, or nth order curve and subtract the total compensation from the total terms at the end. Second order balanced algorithms can require 50 percent more time, but could be worth it for more accuracy ( eg. function tables).
# using pseudocode for Kahan summation algorithm
trying to "cure" or slow growth of tail digits
from round off errors,
# kahan summation algorithm in pseudocode
a(i) = sequence of floating point numbers, a(i)>=3
c = keeper result
sum of floating point numbers
loop following
y = a(i) -c
t = sum + y
c = (t-sum)-y
sum = t
10000.0, 3.14159, 2.71828= (a + b) + c
=> single precision fp errors => 10005.8.
Kahan function a, b, c results in 10005.9
ref normal TCL precision 12 or 17,
expr is double precision inside calculations
and not show these errors in sp.
have to reproduce rd. error in single precision
have to add fp format statements
correct for inputs of large angles and accumularted round off errors
check_answer (a + b) + c =? a + ( b + c ) (yes/no)
set answers and printout with resulting values+----------------------------------------------------------------------------------+ | 1) THE ROUND-OFF ERROR PROBLEM | | Why naive summation loses precision | | | | Input: a=10000.0 b=3.14159 c=2.71828 | | | | Naive (single precision): | | step 1: 10000.0 + 3.14159 --> 10003.1 (tail digits dropped) | | step 2: 10003.1 + 2.71828 --> 10005.8 (error accumulates) | | | | Correct result (double precision): 10005.9 | | | | Problem: (a + b) + c =/= a + (b + c) in single precision fp | | Each addition can drop low-order bits of the running sum | +----------------------------------------------------------------------------------+
+----------------------------------------------------------------------------------+ | 2) KAHAN SUMMATION ALGORITHM (kahan_summation) | | Tracks round-off in a compensation variable c | | | | Initialize: sum = 0, c = 0 (compensation keeper) | | | | For each input a(i) in sequence: | | | | sum, c | | | | | v | | y = a(i) - c <-- subtract accumulated error from term | | | | | v | | t = sum + y <-- add corrected term to running sum | | | | | v | | c = (t - sum) - y <-- capture new round-off into c | | | | | v | | sum = t <-- update running sum | | | | | v (next a(i)) | | | | Result: sum (compensated, closer to true value) | | | | Example: a=10000.0, b=3.14159, c=2.71828 | | Naive sp result: 10005.8 | | Kahan result: 10005.9 <-- matches double precision | +----------------------------------------------------------------------------------+
+----------------------------------------------------------------------------------+ | 3) NEUMAIER VARIANT SUMMATION ALGORITHM (neumaier_summation) | | Improvement over Kahan: handles case where |new term| > |sum| | | | | Initialize: sum = 0, keeperx = 0 (total compensation accumulator) | | | | For each input number in sequence: | | | | sum, number | | | | | v | | t = sum + number | | | | | v | | if |sum| >= |number|: | | keeper = (sum - t) + number <-- sum is larger, correct for its error | | else: | | keeper = (number - t) + sum <-- term is larger, correct for its error | | | | | v | | sum = t <-- update running sum | | keeperx += keeper <-- accumulate all compensations | | | | | v (next number) | | | | Result: sum + keeperx <-- add total compensation at end | +----------------------------------------------------------------------------------+
+----------------------------------------------------------------------------------+
| 4) ALGORITHM COMPARISON |
| |
| Test case A: {10000.0, 3.14159, 2.71828} |
| +--------------------+-------------+-------------+ |
| | Method | Result | Correct? | |
| +--------------------+-------------+-------------+ |
| | Naive sp (+ op) | 10005.8 | no | |
| | Naive sp (expr) | 10005.8 | no | |
| | Kahan | 10005.9 | yes | |
| | Neumaier | 10005.9 | yes | |
| | Normal TCL dp expr | 10005.9 | yes | |
| +--------------------+-------------+-------------+ |
| |
| Test case B: {1.0, +1.E16, 1.0, -1.E16} (large magnitude swings) |
| +--------------------+-------------+-------------+ |
| | Method | Result | Correct? | |
| +--------------------+-------------+-------------+ |
| | Normal TCL dp expr | 2.0 | yes | |
| | Neumaier | 2.0 | yes | |
| +--------------------+-------------+-------------+ |
| Note: Kahan can fail on case B. Neumaier handles both cases correctly. |
+----------------------------------------------------------------------------------++----------------------------------------------------------------------------------+ | 5) WHAT THE COMPENSATION TERM CAPTURES | | | | True value of a(i): [ integer part ][ tail digits ] | | | | After t = sum + y: [ integer part ][000000000000000] <-- tail lost | | | | Captured in c: [ tail digits ] <-- recovered | | | | c = (t - sum) - y | | = what actually got added - what we tried to add | | = the round-off residual, sign-corrected | | | | Next iteration subtracts c from the next term before adding, | | effectively injecting the lost tail back into the sum. | +----------------------------------------------------------------------------------+
+----------------------------------------------------------------------------------+
| 6) TCL SCRIPT ORGANIZATION |
| |
| SetRandomSeed / package requires / namespace path |
| | |
| v |
| +------------------------------------+ |
| | Define test input lists | |
| | kipper = {10000.0 3.14159 2.71828}| |
| +------------------+-----------------+ |
| | |
| +------------+------------+ |
| | | | |
| v v v |
| summer_w_ summer_w_ kahan_ or |
| rd_error rd_expr neumaier_summation |
| (simulate (simulate (compensated |
| sp errors) sp via expr) algorithm) |
| | | | |
| v v v |
| puts results to console (format %7.1f) |
| |
| TCL 8.6 note: expr uses double precision internally. |
| format %7.1f forces single-precision-like display to show rd. errors. |
+----------------------------------------------------------------------------------+ # console program for kahan summation algorithm
# pretty print from autoindent and ased editor
# kahan summation algorithm into TCL
# written on Windows XP on TCL
# working under TCL version 8.6
# gold on TCL WIKI , 8jan2017, kahan_baby3
package require Tk
package require math::numtheory
package require math::geometry
package require math::constants
package require math::bigfloat
namespace path {::tcl::mathop ::tcl::mathfunc math::numtheory math::geometry math::constants math::bigfloat}
set tclprecision 17
wm title . "kahan summation algorithm"
console show
proc summer_w_rd_error { lister } {
set sum 0
foreach number $lister {
set sum [+ $sum [ format %7.1f $number ] ]
}
# statement to recreate single precision r. errors
set sum [ format %7.1f $sum ]
return $sum}
proc summer_w_rd_expr { lister } {
set sum 0
foreach number $lister {
set number [ format %7.1f $number ]
set sum [ expr { $sum + $number } ]
}
# statement to recreate single precision r. errors
set sum [ format %7.1f $sum ]
return $sum}
proc kahan_summation {lister} {
set a {10000.0 3.14159 2.71828 }
set b {9879879.88 6585497.99 875870989.54 765864865479.32 }
set c { 1. 2. 3. }
set compensation_keeper2 0
set sum 0
set counter 0
set yaada2 0
set term2 0
foreach number $lister {
set yaada2 [- $number $compensation_keeper2 ]
set term2 [+ $sum $yaada2 ]
set compensation_keeper2 [- [- $term2 $sum ] $yaada2 ]
set sum $term2
incr $counter}
# statement to recreate single precision r. errors
set sum [ format %7.1f $sum ]
return $sum }
set kipper { 10000.0 3.14159 2.71828 }
puts " normal expr summation, dp inside expr,format %7.1f [format %7.1f [ expr { 10000.0 + 3.14159 + 2.71828 } ] ]"
puts " modify summer_w_rd_error, format %7.1f [ summer_w_rd_error $kipper ] "
puts " modify summer_w_rd_expr, format %7.1f [ summer_w_rd_expr $kipper ] "
puts " trial kahan summation, format %7.1f [ kahan_summation $kipper ] "output
# normal expr summation, dp inside expr,format %7.1f 10005.9
# modify summer_w_rd_error, format %7.1f 10005.8
# modify summer_w_rd_expr, format %7.1f 10005.8
# trial kahan summation, format %7.1f 10005.9 # console program for neumaier summation algorithm
# pretty print from autoindent and ased editor
# neumaier algorithm summation algorithm into TCL
# written on Windows XP on TCL
# working under TCL version 8.6
# gold on TCL WIKI , 8jan2017,neumaier_baby5
package require Tk
package require math::numtheory
package require math::geometry
package require math::constants
package require math::bigfloat
namespace path {::tcl::mathop ::tcl::mathfunc math::numtheory math::geometry math::constants math::bigfloat}
set tclprecision 17
wm title . " neumaier summation algorithm"
console show
#10000.0, 3.14159, 2.71828= (a + b) + c
# (a + b) + c =? a + (b + c)
# => single precision fp errors => 10005.8.
# neumaier algorithm a, b, c results in 10005.9
proc neumaier_summation {lister} {
set a {10000.0 3.14159 2.71828 }
set b {9879879.88 6585497.99 875870989.54 765864865479.32 }
set c { 1. 2. 3. }
set k { 1.0 +1.E100 1.0 -1.E100 }
set keeper 0
set sum 0
set counter 0
set term2 0
set keeperx 0
foreach number $lister {
set term2 [+ $sum $number ]
if { [abs $sum ] >= [abs $number ] } {
set keeper [+ [- $sum $term2 ] $number ] } else {
set keeper [+ [- $number $term2 ] $sum ] }
set sum $term2
set keeper2 $keeper
set keeperx [+ $keeperx $keeper2 ]
incr $counter}
# statement to recreate single precision r. errors
set sum [ format %7.1f $sum ]
set keeper [ format %7.1f $keeper ]
# puts " sum $sum keeper $keeper sum+keeper [+ $sum $keeper ] $sum $keeperx [+ $sum $keeperx ]"
return [ format %7.1f [+ $sum $keeperx ]] }
set kipper { 10000.0 3.14159 2.71828 }
set kipper { 1.0 +1.E16 1.0 -1.E16 }
set kipper { 1. 2. 3. }
set kipper { 10000.0 3.14159 2.71828 }
puts " normal expr summation, dp inside expr,format %7.1f sum { 10000.0 3.14159 2.71828 } [format %7.1f [ expr { 10000.0 + 3.14159 + 2.71828 } ] ]"
puts " normal expr summation, dp inside expr,format %7.1f sum { +1.0 + +1.E15 + +1.0 + -1.E15 } [format %7.1f [ expr { +1.0 + +1.E15 + +1.0 + -1.E15 } ] ]"
puts " trial neumaier summation, format %7.1f sum { 10000.0 3.14159 2.71828 } [ neumaier_summation $kipper ] "
set kipper { 1.0 +1.E16 1.0 -1.E16 }
puts " trial neumaier summation, format %7.1f sum { 1.0 +1.E16 1.0 -1.E16 } [ neumaier_summation $kipper ] "
set kipper { 1. 2. 3. }
puts " trial neumaier summation, format %7.1f sum { 1. 2. 3. } = [ neumaier_summation $kipper ] "
# normal expr summation, dp inside expr,format %7.1f sum { 10000.0 3.14159 2.71828 } 10005.9
# normal expr summation, dp inside expr,format %7.1f sum { +1.0 + +1.E15 + +1.0 + -1.E15 } 2.0
# trial neumaier summation, format %7.1f sum { 10000.0 3.14159 2.71828 } 10005.9
# trial neumaier summation, format %7.1f sum { 1.0 +1.E16 1.0 -1.E16 } 2.0
# trial neumaier summation, format %7.1f sum { 1. 2. 3. } = 6.0 Please place any comments here, Thanks.
| Category Numerical Analysis | Category Toys | Category Calculator | Category Mathematics | Category Example | Toys and Games | Category Games | Category Application | Category GUI |