1# math.tcl -- 2# 3# Collection of math functions. 4# 5# Copyright (c) 1998-2000 by Ajuba Solutions. 6# 7# See the file "license.terms" for information on usage and redistribution 8# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 9# 10# RCS: @(#) $Id: misc.tcl,v 1.6 2005/10/10 14:02:47 arjenmarkus Exp $ 11 12package require Tcl 8.2 ;# uses [lindex $l end-$integer] 13namespace eval ::math { 14} 15 16# ::math::cov -- 17# 18# Return the coefficient of variation of three or more values 19# 20# Arguments: 21# val1 first value 22# val2 second value 23# args other values 24# 25# Results: 26# cov coefficient of variation expressed as percent value 27 28proc ::math::cov {val1 val2 args} { 29 set sum [ expr { $val1+$val2 } ] 30 set N [ expr { [ llength $args ] + 2 } ] 31 foreach val $args { 32 set sum [ expr { $sum+$val } ] 33 } 34 set mean [ expr { $sum/$N } ] 35 set sigma_sq 0 36 foreach val [ concat $val1 $val2 $args ] { 37 set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ] 38 } 39 set sigma_sq [ expr { $sigma_sq/($N-1) } ] 40 set sigma [ expr { sqrt($sigma_sq) } ] 41 if { $mean != 0.0 } { 42 set cov [ expr { ($sigma/$mean)*100 } ] 43 } else { 44 return -code error -errorinfo "Cov undefined for data with zero mean" -errorcode {ARITH DOMAIN} 45 } 46 set cov 47} 48 49# ::math::fibonacci -- 50# 51# Return the n'th fibonacci number. 52# 53# Arguments: 54# n The index in the sequence to compute. 55# 56# Results: 57# fib The n'th fibonacci number. 58 59proc ::math::fibonacci {n} { 60 if { $n == 0 } { 61 return 0 62 } else { 63 set prev0 0 64 set prev1 1 65 for {set i 1} {$i < $n} {incr i} { 66 set tmp $prev1 67 incr prev1 $prev0 68 set prev0 $tmp 69 } 70 return $prev1 71 } 72} 73 74# ::math::integrate -- 75# 76# calculate the area under a curve defined by a set of (x,y) data pairs. 77# the x data must increase monotonically throughout the data set for the 78# calculation to be meaningful, therefore the monotonic condition is 79# tested, and an error is thrown if the x value is found to be 80# decreasing. 81# 82# Arguments: 83# xy_pairs list of x y pairs (eg, 0 0 10 10 20 20 ...); at least 5 84# data pairs are required, and if the number of data 85# pairs is even, a padding value of (x0, 0) will be 86# added. 87# 88# Results: 89# result A two-element list consisting of the area and error 90# bound (calculation is "Simpson's rule") 91 92proc ::math::integrate { xy_pairs } { 93 94 set length [ llength $xy_pairs ] 95 96 if { $length < 10 } { 97 return -code error "at least 5 x,y pairs must be given" 98 } 99 100 ;## are we dealing with x,y pairs? 101 if { [ expr {$length % 2} ] } { 102 return -code error "unmatched xy pair in input" 103 } 104 105 ;## are there an even number of pairs? Augment. 106 if { ! [ expr {$length % 4} ] } { 107 set xy_pairs [ concat [ lindex $xy_pairs 0 ] 0 $xy_pairs ] 108 } 109 set x0 [ lindex $xy_pairs 0 ] 110 set x1 [ lindex $xy_pairs 2 ] 111 set xn [ lindex $xy_pairs end-1 ] 112 set xnminus1 [ lindex $xy_pairs end-3 ] 113 114 if { $x1 < $x0 } { 115 return -code error "monotonicity broken by x1" 116 } 117 118 if { $xn < $xnminus1 } { 119 return -code error "monotonicity broken by xn" 120 } 121 122 ;## handle the assymetrical elements 0, n, and n-1. 123 set sum [ expr {[ lindex $xy_pairs 1 ] + [ lindex $xy_pairs end ]} ] 124 set sum [ expr {$sum + (4*[ lindex $xy_pairs end-2 ])} ] 125 126 set data [ lrange $xy_pairs 2 end-4 ] 127 128 set xmax $x1 129 set i 1 130 foreach {x1 y1 x2 y2} $data { 131 incr i 132 if { $x1 < $xmax } { 133 return -code error "monotonicity broken by x$i" 134 } 135 set xmax $x1 136 incr i 137 if { $x2 < $xmax } { 138 return -code error "monotonicity broken by x$i" 139 } 140 set xmax $x2 141 set sum [ expr {$sum + (4*$y1) + (2*$y2)} ] 142 } 143 144 if { $xmax > $xnminus1 } { 145 return -code error "monotonicity broken by xn-1" 146 } 147 148 set h [ expr { ( $xn - $x0 ) / $i } ] 149 set area [ expr { ( $h / 3.0 ) * $sum } ] 150 set err_bound [ expr { ( ( $xn - $x0 ) / 180.0 ) * pow($h,4) * $xn } ] 151 return [ list $area $err_bound ] 152} 153 154# ::math::max -- 155# 156# Return the maximum of two or more values 157# 158# Arguments: 159# val first value 160# args other values 161# 162# Results: 163# max maximum value 164 165proc ::math::max {val args} { 166 set max $val 167 foreach val $args { 168 if { $val > $max } { 169 set max $val 170 } 171 } 172 set max 173} 174 175# ::math::mean -- 176# 177# Return the mean of two or more values 178# 179# Arguments: 180# val first value 181# args other values 182# 183# Results: 184# mean arithmetic mean value 185 186proc ::math::mean {val args} { 187 set sum $val 188 set N [ expr { [ llength $args ] + 1 } ] 189 foreach val $args { 190 set sum [ expr { $sum + $val } ] 191 } 192 set mean [expr { double($sum) / $N }] 193} 194 195# ::math::min -- 196# 197# Return the minimum of two or more values 198# 199# Arguments: 200# val first value 201# args other values 202# 203# Results: 204# min minimum value 205 206proc ::math::min {val args} { 207 set min $val 208 foreach val $args { 209 if { $val < $min } { 210 set min $val 211 } 212 } 213 set min 214} 215 216# ::math::product -- 217# 218# Return the product of one or more values 219# 220# Arguments: 221# val first value 222# args other values 223# 224# Results: 225# prod product of multiplying all values in the list 226 227proc ::math::product {val args} { 228 set prod $val 229 foreach val $args { 230 set prod [ expr { $prod*$val } ] 231 } 232 set prod 233} 234 235# ::math::random -- 236# 237# Return a random number in a given range. 238# 239# Arguments: 240# args optional arguments that specify the range within which to 241# choose a number: 242# (null) choose a number between 0 and 1 243# val choose a number between 0 and val 244# val1 val2 choose a number between val1 and val2 245# 246# Results: 247# num a random number in the range. 248 249proc ::math::random {args} { 250 set num [expr {rand()}] 251 if { [llength $args] == 0 } { 252 return $num 253 } elseif { [llength $args] == 1 } { 254 return [expr {int($num * [lindex $args 0])}] 255 } elseif { [llength $args] == 2 } { 256 foreach {lower upper} $args break 257 set range [expr {$upper - $lower}] 258 return [expr {int($num * $range) + $lower}] 259 } else { 260 set fn [lindex [info level 0] 0] 261 error "wrong # args: should be \"$fn ?value1? ?value2?\"" 262 } 263} 264 265# ::math::sigma -- 266# 267# Return the standard deviation of three or more values 268# 269# Arguments: 270# val1 first value 271# val2 second value 272# args other values 273# 274# Results: 275# sigma population standard deviation value 276 277proc ::math::sigma {val1 val2 args} { 278 set sum [ expr { $val1+$val2 } ] 279 set N [ expr { [ llength $args ] + 2 } ] 280 foreach val $args { 281 set sum [ expr { $sum+$val } ] 282 } 283 set mean [ expr { $sum/$N } ] 284 set sigma_sq 0 285 foreach val [ concat $val1 $val2 $args ] { 286 set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ] 287 } 288 set sigma_sq [ expr { $sigma_sq/($N-1) } ] 289 set sigma [ expr { sqrt($sigma_sq) } ] 290 set sigma 291} 292 293# ::math::stats -- 294# 295# Return the mean, standard deviation, and coefficient of variation as 296# percent, as a list. 297# 298# Arguments: 299# val1 first value 300# val2 first value 301# args all other values 302# 303# Results: 304# {mean stddev coefvar} 305 306proc ::math::stats {val1 val2 args} { 307 set sum [ expr { $val1+$val2 } ] 308 set N [ expr { [ llength $args ] + 2 } ] 309 foreach val $args { 310 set sum [ expr { $sum+$val } ] 311 } 312 set mean [ expr { $sum/double($N) } ] 313 set sigma_sq 0 314 foreach val [ concat $val1 $val2 $args ] { 315 set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ] 316 } 317 set sigma_sq [ expr { $sigma_sq/double($N-1) } ] 318 set sigma [ expr { sqrt($sigma_sq) } ] 319 set cov [ expr { ($sigma/$mean)*100 } ] 320 return [ list $mean $sigma $cov ] 321} 322 323# ::math::sum -- 324# 325# Return the sum of one or more values 326# 327# Arguments: 328# val first value 329# args all other values 330# 331# Results: 332# sum arithmetic sum of all values in args 333 334proc ::math::sum {val args} { 335 set sum $val 336 foreach val $args { 337 set sum [ expr { $sum+$val } ] 338 } 339 set sum 340} 341 342#---------------------------------------------------------------------- 343# 344# ::math::expectDouble -- 345# 346# Format an error message that an argument was expected to be 347# double and wasn't 348# 349# Parameters: 350# arg -- Misformatted argument 351# 352# Results: 353# Returns an appropriate error message 354# 355# Side effects: 356# None. 357# 358#---------------------------------------------------------------------- 359 360proc ::math::expectDouble { arg } { 361 return [format "expected a floating-point number but found \"%.50s\"" $arg] 362} 363 364#---------------------------------------------------------------------- 365# 366# ::math::expectInteger -- 367# 368# Format an error message that an argument was expected to be 369# integer and wasn't 370# 371# Parameters: 372# arg -- Misformatted argument 373# 374# Results: 375# Returns an appropriate error message 376# 377# Side effects: 378# None. 379# 380#---------------------------------------------------------------------- 381 382proc ::math::expectInteger { arg } { 383 return [format "expected an integer but found \"%.50s\"" $arg] 384} 385 386