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:
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
} ;#-)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)}
} ;# RSBinomial 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
} ;#RSLinear 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
} ;#RSSign of a number:
proc sgn {a} {expr {$a>0 ? 1 : $a<0 ? -1 : 0}} ;# rmaxActually,
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)
Should this function move to the Stats page? Means of a number list: (arithmetic, geometric, quadratic, harmonic)
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}] } ;# DKFLogarithm to any base:
proc log {base x} {expr {log($x)/log($base)}} ;# RSA 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-019CritLib (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
AM On the c.l.t. the other day, 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, 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 assing d to the result */
d = 1.0;
resultPtr->type = TCL_DOUBLE;
resultPtr->doubleValue = d;
return TCL_OK;
}in you package initialisation...
Tcl_ValueType argblogArgs[2] = { TCL_DOUBLE, TCL_DOUBLE };
Tcl_CreateMathFunc(interp, "arblog", 2, arblogArgs, ArgLogProc,
(ClientData)NULL);"
Mathematically oriented extensions - Arts and crafts of Tcl-Tk programming - Category Mathematics