Additional math functions

Average edit

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

Binomial coefficient edit

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
}

Cross-sum of non-negative integers edit

proc crosssum {x} {expr [join [split $x ""] +]}

Note that this expression may not be braced. (RS)

Epsilon edit

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

Factorial edit

RS 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)
# 120

I 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:

It returns exact results for factorial x, where x is an integer and 0<=x<=21.
It returns floating point results for integer x, 22<=x<=170, that are correct to 1 unit in the least significant bit position.
It returns approximate results, precise to nine significant digits, for all other real x, x>=0, by using the identity x! = Gamma( x + 1 ). In particular, this precision has been exhaustively verified for all half-integer arguments that give results within the range of IEEE floating point.
It has companion functions for binomial coefficients, the Gamma function and the Beta distribution that are as precise as it is. Moreover, these functions do not suffer from premature overflow; they perform well with large arguments: [choose 10000 100] doesn't give the function heartburn.

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 
} ;#-)

Fibonacci numbers edit

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))))} 
} ;# RS

Lars 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).

Integer Check edit

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

Integer maximum edit

(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}
}

Linear regression and correlation coefficient edit

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
} ;#RS

Logarithm to any base edit

proc log {base x} {
    expr {log($x)/log($base)}
} ;# RS

Logarithm to Base Two edit

proc 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)

Maximum and minimum edit

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.

Means of a number list: arithmetic, geometric, quadratic, harmonic edit

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}] } ;# DKF

JPS: 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}]
  }
}

Mid edit

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 [1] 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).

Numerical functions for [expr] edit

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 [2] - JCW

Prime factors of an integer edit

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
} ;#RS

Random Numbers edit

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

Square Mean and Standard Deviation edit

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)}
} ;# RS

Sign of a number edit

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 again

Actually,

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

Traditional degrees edit

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)

Category Mathematics

Arts and Crafts of Tcl-Tk Programming

See Also edit