Version 0 of Binomial Coefficients

Updated 2001-06-13 10:15:31

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