Version 48 of Additional math functions

Updated 2012-11-13 07:49:51 by pooryorick

See Also

Integer Check

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

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

Maximum and minimum

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.

Average

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

Random Numbers

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.

Factorial

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

Square Mean and Standard Deviation

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

Binomial coefficient

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
}

Prime factors of an integer

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

Linear regression and correlation coefficient

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

Sign of a number

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.

Traditional degrees

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)

Cross-sum of non-negative integers

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

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

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

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

Logarithm to any base

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

A faster logarithm to base two

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)

Epsilon

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

Numerical functions for [expr]

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

** Defining New Math Functions 8*

AM On the c.l.t. the other day [is this as of May 2003?], Martin Russell asked about how to define new math functions. If you want to do it without the help of DKF's extension [??] and CrtLib [ critcl's critlib?], then here is a receipe provided by Pat Thoyts:

Something along these lines.

static Tcl_MathProc ArbLogProc;

static int
ArbLogProc(clientData, interp, args, resultPtr)
  ClientData clientData;
  Tcl_Interp *interp;     /* current interpreter */
  Tcl_Value  *args;       /* input arguments */
  Tcl_Value  *resultPtr;  /* where to store the result */
{
   double b, n, d;
   b = args[0].doubleValue;
   n = args[1].doubleValue;

   /* do your maths and assign d to the result */
   d = 1.0;

   resultPtr->type = TCL_DOUBLE;
   resultPtr->doubleValue = d;
   return TCL_OK;
}

in your package initialisation...

 Tcl_ValueType arblogArgs[2] = { TCL_DOUBLE, TCL_DOUBLE };
 Tcl_CreateMathFunc(interp, "arblog", 2, arblogArgs, ArbLogProc,
                    (ClientData)NULL);

" In Tcl 8.5, math functions are all located in a namespace, 'tcl::mathfunc' which is resolved relative to the current namespace (so either ::tcl::matfunc::f or [namespace current]::tcl::mathfunc::f can resolve f($x)). Log to an arbitrary base can therefore be done with:

proc tcl::mathfunc::logbase {x b} {
    expr {log($x) / log($b)}
}

without any C hackery being needed.

Fibonacci numbers:

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