[MS]: moved here the discussion/proposals from [Additional math functions] '''Binomial coefficient:''' proc over {n k} { # binomial coefficient: e.g. [over 5 3] = (5*4*3)/(1*2*3) = 10 set k [expr {(($n-$k) < $k) ? $n-$k : $k}] if {$k==0} {return 1} set num $n set den $k while {$k>1} { incr n -1 incr k -1 set num [expr {$num*$n}] set den [expr {$den*$k}] } expr {$num/$den} } ;#RS Note that combinatoric functions are notorious for rewarding algorithmic "trickiness". In this case, the multiple "expr" invocations will toggle to floating point when n is sufficiently large. This is both slow and inaccurate, in general, although certainly quite adequate for many common cases of interest. Higher precision is available with algorithms which exploit common factors repeated in the numerator and denominator (and there will be many of them for most calculations). [CL] will make a point some day of returning here to offer an [[over]] written in terms of [[primefactors]] (see below). - [RS]: Here's a variant which uses mostly lists: collect factors for ''num'' and ''den'', remove from ''num'' all that are in ''den'', finally computes one integer expression: proc over2 {n k} { # binomial coefficient: e.g. [over 5 3] = (5*4*3)/(1*2*3) = 10 set k [expr {(($n-$k) < $k) ? $n-$k : $k}] if {$k==0} {return 1} set num $n set den $k while {$k>1} { lappend num [incr n -1] lappend den [incr k -1] } set t 1 ;# neutral element of multiplication foreach i $num { if {[set pos [lsearch $den $i]]<0} { lappend t $i } else { set den [lreplace $den $pos $pos] } } expr [join $t *]/([join $den *]) } Recursive binomial calculations can be surprisingly well-conditioned (someone compare performance on these): proc binom {m n} { set diff [expr $m - $n] if {$diff < $n} { return [binom $m $diff] } switch $n { 0 { return 1 } 1 { return $m } default { # This could use a bit of common-factor elimination. return [expr [binom [expr $m - 1] [expr $n - 1]] * $m / $n] } } } ''MS'': This looks like too much of a good thing, you do not need to recurse over both variables. Recursing over just ''m'' gives the following algorithm. proc binom2 {m n} { set n [expr {(($m-$n) > $n) ? $m-$n : $n}] set k [expr {$m - $n}] if {$k < 0} {return 0} switch $k { 0 {return 1} 1 {return $m} default { return [expr {($m*[binom [expr {$m - 1}] $n])/$k}] } } } ;#MS Both recursive algorithms are stable in the sense of not producing unnecessarily large intermediate results. If you are interested in performance, avoid recursion and [[switch]]. An iterative implementation of ''binom2'' is: proc binom3 {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 } ;#MS Remark that this is very much the same as ''over'', with the loop being traversed in reverse order and avoiding the accumulation in ''den'' and ''num'' which causes an earlier overflow. --- ''MS'': So, which is best? While trying to compute all coefficients in the expansion of (a+b)^20, both [[over]] and [[over2]] fail to produce correct results as some intermediate computations overflow. The timing of the recursive/iterative solutions gives: % time {for {set i 0} {$i <= 20} {incr i} {set x [binom 20 $i]}} 51716 microseconds per iteration % time {for {set i 0} {$i <= 20} {incr i} {set x [binom2 20 $i]}} 47071 microseconds per iteration % time {for {set i 0} {$i <= 20} {incr i} {set x [binom3 20 $i]}} 3490 microseconds per iteration --- % time {for {set i 0} {$i <= 20} {incr i} {set x [binom4 20 $i]}} 36949 microseconds per iteration The best stability (but low speed) is obtained by using set operations on the lists of prime factors of numerator and denominator: proc binom4 {m n} { set n [expr {(($m-$n) > $n) ? $m-$n : $n}] if {$n > $m} {return 0} if {$n == $m} {return 1} # Find the list of prime factors in the numerator and denominator set num {} set den {} set d 0 while {$n < $m} { set num [appendPrimeFactors [incr n] $num] set den [appendPrimeFactors [incr d] $den] } # We now need to multiply all factors in $num that are not in $den set num [lsort -integer $num] set den [lsort -integer $den] set res 1 set i 0 ;# position counter for num foreach f $den { while {[set g [lindex $num $i]] != $f} { set res [expr {$res * $g}] incr i } incr i } foreach g [lrange $num $i end] { set res [expr {$res * $g}] } set res } ;#MS proc appendPrimeFactors {n res} { # a number x is prime if [llength [primefactors $x]]==1 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 There should however exist ''very few'' cases where binom3 overflows and binom4 doesn't ... But, if you're worried, read on: ---- '''Log gamma -- or, factorials the hard way.''' There was a little discussion on the chat about how to do factorials inside [[expr]]. [Larry Virden] pointed the questioner to several pages with classical recursive or iterative approaches. [Kevin Kenny] saw a challenge in the discussion, and decided to do factorials the hard way -- by computing the gamma function. This is actually useful for avoiding overflows in computing binomial coefficients. It can also be used as an obscure way of computing pi. # Returns log(gamma(x)), where x >= 0 # # +inf # _ # | z-1 -t # gamma(x)= _| t e dt # # 0 proc lgamma { xx } { set x [expr { $xx - 1.0 }] set tmp [expr { $x + 5.5 }] set tmp [ expr { ( $x + 0.5 ) * log( $tmp ) - $tmp }] set ser 1.0 foreach cof { 76.18009173 -86.50532033 24.01409822 -1.231739516 .00120858003 -5.36382e-6 } { set x [expr { $x + 1.0 }] set ser [expr { $ser + $cof / $x }] } return [expr { $tmp + log( 2.50662827465 * $ser ) }] } # Return x! given x. proc fac { x } { return [expr { exp( [lgamma [expr { $x + 1 }]] ) }] } # Demonstrate by printing a table of factorials for { set i 0 } { $i <= 10 } { incr i } { puts "$i! == [expr { round( [fac $i] ) }]" } # Compute binomial coefficients. proc choose { n k } { set r [expr { [lgamma [expr { $n + 1 }]] - [lgamma [expr { $k + 1 }]] - [lgamma [expr { $n - $k + 1 }]] }] return [expr { exp( $r ) }] } # Demonstrate binomial coefficients by printing Pascal's triangle. puts "Pascal's triangle:" for { set n 0 } { $n < 15 } { incr n } { set line {} for { set k 0 } { $k <= $n } { incr k } { set r [expr { round( [choose $n $k] ) }] append line [format %5d $r] } puts $line } # Compute pi, the hard way! set sqrtpi [expr { exp( [lgamma 0.5] ) }] puts "pi is approximately [expr { $sqrtpi * $sqrtpi }]" ---- [Category Mathematics]