<> ** See Also ** * [math] module in [tcllib] * [integrate] * [Stats] * [converting between rectangular and polar co-ordinates] * [Fraction math] - [Complex math made simple] (Complex numbers) * [Sample Math Programs] * [Binomial coefficients] - including the gamma function. * [Combinatorial mathematics functions] * [Fast Fourier Transform] * [Multivariate Linear Regression] * [Mathematically oriented extensions] * [Arts and crafts of Tcl-Tk programming] ** 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 [http://www.equi4.com/critlib/mathf.README] - [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. ======c 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... ======c 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). <> Mathematics