arithmetic mean of a list of numbers:
proc average L {
expr ([join $L +])/[llength $L].
}Note that empty lists produce a syntax error. The dot behind llength casts it to double (not dangerous here, as llength will always return a non-negative integer) -- RS
Perhaps the best (what criteria?) should move to the Binomial page and just a pointer to the page should be here? (This got too long; I'm keeping the best algorithm here, moving the previous discussion to Binomial Coefficients. This solution is called binom3 in that page.)
proc binom {m n} {
set n [expr {(($m-$n) > $n) ? $m-$n : $n}]
if {$n > $m} {return 0}
if {$n == $m} {return 1}
set res 1
set d 0
while {$n < $m} {
set res [expr {($res*[incr n])/[incr d]}]
}
set res
}proc crosssum {x} {expr [join [split $x ""] +]}Note that this expression may not be braced. (RS)
Comparing two floats x,y for equality is most safely done by testing abs($x-$y)<$eps, where eps is a sufficiently small number. You can find out which eps is good for your machine with the following code:
proc eps {{base 1}} {
set eps 1e-20
while {$base-$eps==$base} {
set eps [expr {$eps+1e-22}]
}
set eps [expr {$eps+1e-22}]
}
% eps 1
5.55112000002e-017 ;# on both my Win2K/P3 and Sun/Solaris
% eps 0.1
6.93889999999e-018
% eps 0.01
8.674e-019
% eps 0.001
1.085e-019RS 2008-01-02: Here's a little example for a user-defined recursive factorial function:
proc tcl::mathfunc::fac x {
expr {$x<2? 1: $x*fac($x-1)}
}
expr fac(5)
# 120I see a factorial function on 3-4 different pages -some not even about math. And yet none in the tcllib math library. Perhaps one should be submitted. How to determine best?
KBK: There is indeed a factorial in ::tcllib::math. It's in some sense 'better' than any of the ones I've seen here on the Wiki:
RS I like this one, compact but recursive:
proc fac n {
expr {$n<2? 1: $n*[fac [expr {$n-1}]]}
}However, this one runs 1/3 faster:
proc fac2 n {
expr $n<2? 1: [join [iota 1 $n] *]+0
}given an index generator iota, e.g. iota 1 5 => {1 2 3 4 5}
proc iota {base n} {
set res {}
for {set i $base} {$i<$n+$base} {incr i} {lappend res $i}
set res
}However, factorials computed in terms of expr are correct only until 12!; above that you get "false positives", negatives, or zeroes.. Of course one could use doubles, which seem to be exact up to 18! (at the maximum tcl_precision 17). But the fastest fac is still tabulated:
proc fac3 n {
lindex {
1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600
479001600.0 87178291200.0 1307674368000.0 20922789888000.0
355687428096000.0 6402373705728000.0
} $n
} ;#-)tcllib::math has an iterative version, but here's the "closed form" if anyone cares:
proc fib n {
expr {round(1/sqrt(5)*(pow((1+sqrt(5))/2,$n) - (pow((1-sqrt(5))/2,$n))))}
} ;# RSLars H: Actually, you don't need to compute the second term, since it always contributes < 1/2 for non-negative n. You can simply do
proc fib2 n {
expr {round(1/sqrt(5)*pow((1+sqrt(5))/2,$n))}
}For negative n it is instead the first term that can be ignored, but one rarely needs those Fibonacci numbers.
BTW, I also changed an "int" to a "round" in RS's proc (if you're unlucky with the numerics, "int" can give you one less than the correct answer).
see whether variable has an integer value
Since Tcl 8.1.1, the built-in string is int does the same for a value.
proc is_int x {
expr {![catch {incr x 0}]}
}
proc is_no_int x {
catch {incr x 0}
}(MAXINT): determine biggest positive signed integer (by Jeffrey Hobbs):
proc largest_int {} {
set int 1
set exp 7; # assume we get at least 8 bits
while {$int > 0} { set int [expr {1 << [incr exp]}] }
expr {$int-1}
}proc reg,cor points {
# linear regression y=ax+b for {{x0 y0} {x1 y1}...}
# returns {a b r}, where r: correlation coefficient
foreach i {N Sx Sy Sxy Sx2 Sy2} {set $i 0.0}
foreach point $points {
foreach {x y} $point break
set Sx [expr {$Sx + $x}]
set Sy [expr {$Sy + $y}]
set Sx2 [expr {$Sx2 + $x*$x}]
set Sy2 [expr {$Sy2 + $y*$y}]
set Sxy [expr {$Sxy + $x*$y}]
incr N
}
set t1 [expr {$N*$Sxy - $Sx*$Sy}]
set t2 [expr {$N*$Sx2 - $Sx*$Sx}]
set a [expr {double($t1)/$t2}]
set b [expr {double($Sy-$a*$Sx)/$N}]
set r [expr {$t1/(sqrt($t2)*sqrt($N*$Sy2-$Sy*$Sy))}]
list $a $b $r
} ;#RSproc log {base x} {
expr {log($x)/log($base)}
} ;# RSproc ld x "expr {log(\$x)/[expr log(2)]}"This is an example of a "live" proc body - the divisor is computed only once, at definition time. With a single backslash escape needed, it's worth the fun ;-) (RS)
proc max {a args} {
foreach i $args {if {$i>$a} {set a $i}};return $a
}
proc min {a args} {
foreach i $args {if {$i<$a} {set a $i}};return $a
}Works with whatever < and > can compare (strings included). Or how about (float numbers only):
proc max args {
lindex [lsort -real $args] end
}
proc min args {
lindex [lsort -real $args] 0
}Or, use -dictionary to handle strings, ints, real.... and also allow to be called with a single list arg (FYI, it's actually a bit faster to use the sort method)
proc min args {
if {[llength $args] == 1} {set args [lindex $args 0]}
lindex [lsort -dict $args] 0
}
proc max args {
if {[llength $args] == 1} {set args [lindex $args 0]}
lindex [lsort -dict $args] end
}RS: ... only that you get lsort results like
{-1 -5 -10 0 5 10}if you use the -dict mode of lsort. Numeric max/min should rather use -integer or -float. Max/min of strings must be left to dedicated procs, if ever needed.
Should this function move to the Stats page?
proc mean L {
expr ([join $L +])/[llength $L].
}
proc gmean L {
expr pow([join $L *],1./[llength $L])
}
proc qmean L {
expr sqrt((pow([join $L ,2)+pow(],2))/[llength $L])
}
proc hmean L {
expr [llength $L]/(1./[join $L +1./])
}where qmean is the best braintwister... For a list of {1 2} the string
sqrt((pow( 1 ,2)+pow( 2 ,2))/ 2)
(blanks added for clarity) is built up and fed to expr, where it makes a perfectly well-formed expression if not braced. (RS)
proc median L {lindex $L [expr {[llength $L]/2}] } ;# DKFJPS: That median assumes the list is already sorted. This one doesn't:
proc median {l} {
if {[set len [llength $l]] % 2} then {
return [lindex [lsort -real $l] [expr {($len - 1) / 2}]]
} else {
return [expr {([lindex [set sl [lsort -real $l]] [expr {($len / 2) - 1}]] \
+ [lindex $sl [expr {$len / 2}]]) / 2.0}]
}
}AMG: Here's a math function I sometimes find useful. It accepts three arguments, and it returns whichever of the three is between the other two. It's mostly useful to clamp a number to a range.
proc ::tcl::mathfunc::mid {a b c} {
lindex [lsort -real [list $a $b $c]] 1
}It can also be implemented as a bunch of [if]s, which is how I do it in C.
Here is one incorrect implementation you should watch out for:
proc ::tcl::mathfunc::mid {a b c} {
expr {max($a, min($b, $c))}
}This is what Allegro (include/allegro/base.h) has used since the dawn of time. :^( I'm reporting it now; hopefully it'll be fixed. If you're curious, see [L1 ] for my writeup.
KPV The folk algorithm for finding the middle number (or second highest in a longer list) is to take the max of the pair-wise mins. To wit:
max(min($a,$b), min($a,$c), min($b,$c))
LV So what is an example of a case in which the second, incorrect, version of the algorithm fails? Answer: "incorrect_mid 1 0 0" returns 1. The problem is it doesn't (always) handle the case where two of the inputs are the same. Doh.
AMG: I thought the problem was that it doesn't handle the case of the first input being greater than the other two. This wasn't a problem for Allegro because everyone used its MID macro thus: MID(minimum_value, value_to_clamp, maximum_value).
CritLib (see the Critcl page) now includes an adapted version of Donal K. Fellows' extension which lets you write numerical functions for "expr" in Tcl. See the "mathf" readme [L2 ] - JCW
proc primefactors n {
# a number x is prime if [llength [primefactors $x]]==1
set res {}
set f 2
while {$f<=$n} {
while {$n%$f==0} {
set n [expr {$n/$f}]
lappend res $f
}
set f [expr {$f+2-($f==2)}]
}
set res
} ;#RSOf course, since 8.0 just say
expr {rand()}Jeffrey Hobbs has this substitute for pre-8.0 Tcl:
set _ran [clock seconds]
proc random {range} {
global _ran
set _ran [expr ($_ran * 9301 + 49297) % 233280]
return [expr int($range * ($_ran / double(233280)))]
}Pass in an int and it returns a number (0..int). Also, the Wiki page on "rand" has more on the subject.
Perhaps this function should move to the Stats page mentioned above? Square mean and standard deviation:
proc mean2 list {
set sum 0
foreach i $list {set sum [expr {$sum+$i*$i}]}
expr {double($sum)/[llength $list]}
}
proc stddev list {
set m [mean $list] ;# see below for [mean]
expr {sqrt([mean2 $list]-$m*$m)}
} ;# RSLWS 19 Feb 2021: This stddev gives different results than ::math::statistics::stdev (the latter of which matches a spreadsheet calculation). I'm not sure why, but thought I would point it out in case anyone was going to use it.
proc sgn {a} {expr {$a>0 ? 1 : $a<0 ? -1 : 0}} ;# rmax
proc sgn x {expr {$x<0? -1: $x>0}} ;# RS
proc sgn x {expr {($x>0)+($x>>31)}} ;# jcw (32-bit arch)
proc sgn x {expr {($x>0)-($x<0)}} ;# rmax againActually,
string compare $a 0
seems to give the correct result for all integer values and floating point values not equal to 0.
0.0 (and 0.00 etc) [string compare 0.0 0] returns 1, however.
sergiol: I was playing codegolf and based on the C answer http://codegolf.stackexchange.com/a/103831/29325 I confirmed it on tcl
puts [expr !!$n|$n>>31]
can be seen on: http://rextester.com/live/BKGZ8868
clock format can be put to un-timely uses. As degrees especially in geography are also subdivided in minutes and seconds, how's this one-liner for formatting decimal degrees:
proc dec2deg x {
concat [expr int($x)] [clock format [expr round($x*3600)] -format "%M' %S\""]
}An additional -gmt 1 switch is needed if you happen to live in a non-integer timezone. (RS)
When dealing with hardware registers, communication protocols or low-level drivers there's often a need to saturate values and convert them to an integer having a particular bit size.
The following math function does that:
# Saturate X and convert the result to an integer value.
#
# N defines the bit size of the integer value:
# - If N >= 0:
# Positive values of N (and zero) define a range from 0 to 2**N-1.
# This is the range of unsigned integer values having N bits.
# For example if N = 8 the range is 0 to 255.
# - If N < 0:
# Negative values of N define a range from -2**abs(N+1) to 2**abs(N+1)-1.
# This is the range of signed integer values having abs(N) bits in two's complement.
# For example if N = -8 the range is -128 to 127.
#
# If X is -inf the result is the maximum negative value in given range.
# If X is +inf the result is the maximum positive value in given range.
# If X is nan an error is raised.
proc ::tcl::mathfunc::satint {x n} {
if {![string is entier -strict $n]} {
error "N must be an integer value"
}
if {$n >= 0} {
# unsigned integer range
set a 0
set z [expr {2**$n - 1}]
} else {
# signed integer range (two's complement)
set a [expr {-2**abs($n + 1)}]
set z [expr {2**abs($n + 1) - 1}]
}
return [expr {entier(min($z, max($a, $x)))}]
}Usage
;# N >= 0: unsigned integer range
% expr {satint(1234, 8)}
255
% expr {satint(-1, 8)}
0
% expr {satint(123.6, 8)}
123
% expr {satint(+inf, 16)}
65535
% expr {satint(1, 8.5)}
N must be an integer value
;# N < 0: signed integer range (two's complement)
% expr {satint(-200, -8)}
-128
% expr {satint(1234, -8)}
127
% expr {satint(123.6, -8)}
123The rounding mode (see Rounding in Tcl) can be changed by combining with math functions like ceil, floor or round:
expr {satint(round(123.6), 8)}
124
expr {satint(floor(123.6), 8)}
123
expr {satint(ceil(123.6), 8)}
124This is a math function computing the median value:
proc ::tcl::mathfunc::median {x args} {
# expand x to allow handling a list argument
set p [lsort -real [list {*}$x {*}$args]]
set n [llength $p]
set i [expr {$n >> 1}]
if {0 != ($n & 1)} {
set y [lindex $p $i]
} else {
# for an even number of elements the median is
# the average of the two center elements
set a [lindex $p $i]
set b [lindex $p [expr {$i - 1}]]
set y [expr {0.5 * ($a + $b)}]
}
return $y
}Note that this implementation can handle list arguments:
expr {median(4,3,1,2)}
% -> 2.5
set x {1 3 2}
expr {median($x)}
% -> 2
if {median($x) == 2} { puts "It works!" }
% -> It works!