combinatorics.tcl

Go to the documentation of this file.

/* ----------------------------------------------------------------------*/
/* */
/*  math/combinatorics.tcl --*/
/* */
/*  This file contains definitions of mathematical functions*/
/*  useful in combinatorial problems.  */
/* */
/*  Copyright (c) 2001, by Kevin B. Kenny.  All rights reserved.*/
/* */
/*  See the file "license.terms" for information on usage and redistribution*/
/*  of this file, and for a DISCLAIMER OF ALL WARRANTIES.*/
/*  */
/*  RCS: @(#) $Id: combinatorics.tcl,v 1.5 2004/02/09 19:31:54 hobbs Exp $*/
/* */
/* ----------------------------------------------------------------------*/

package require Tcl 8.0

namespace ::math {

    /*  Commonly used combinatorial functions*/

    /*  ln_Gamma is spelt thus because it's a capital gamma (\u0393)*/

    namespace export ln_Gamma;      /*  Logarithm of the Gamma function*/
    namespace export factorial;     /*  Factorial*/
    namespace export choose;        /*  Binomial coefficient*/

    /*  Note that Beta is spelt thus because it's conventionally a*/
    /*  capital beta (\u0392).  It is exported from the package even*/
    /*  though its name is capitalized.*/

    namespace export Beta;      /*  Beta function*/

}

/* ----------------------------------------------------------------------*/
/* */
/*  ::math::InitializeFactorial --*/
/* */
/*  Initialize a table of factorials for small integer arguments.*/
/* */
/*  Parameters:*/
/*  None.*/
/* */
/*  Results:*/
/*  None.*/
/* */
/*  Side effects:*/
/*  The variable, ::math::factorialList, is initialized to hold*/
/*  a table of factorial n for 0 <= n <= 170.*/
/* */
/*  This procedure is called once when the 'factorial' procedure is*/
/*  being loaded.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::InitializeFactorial () {

    variable factorialList

    set factorialList [list 1]
    set f 1
    for { set i 1 } { $i < 171 } { incr i } {
    if { $i > 12. } {
        set f [expr { $f * double($i)}]
    } else {
        set f [expr { $f * $i }]
    }
    lappend factorialList $f
    }

}

/* ----------------------------------------------------------------------*/
/* */
/*  ::math::InitializePascal --*/
/* */
/*  Precompute the first few rows of Pascal's triangle and store*/
/*  them in the variable ::math::pascal*/
/* */
/*  Parameters:*/
/*  None.*/
/* */
/*  Results:*/
/*  None.*/
/* */
/*  Side effects:*/
/*  ::math::pascal is initialized to a flat list containing */
/*  the first 34 rows of Pascal's triangle.  C(n,k) is to be found*/
/*  at [lindex $pascal $i] where i = n * ( n + 1 ) + k.  No attempt*/
/*  is made to exploit symmetry.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::InitializePascal () {

    variable pascal

    set pascal [list 1]
    for { set n 1 } { $n < 34 } { incr n } {
    lappend pascal 1
    set l2 [list 1]
    for { set k 1 } { $k < $n } { incr k } {
        set km1 [expr { $k - 1 }]
        set c [expr { [lindex $l $km1] + [lindex $l $k] }]
        lappend pascal $c
        lappend l2 $c
    }
    lappend pascal 1
    lappend l2 1
    set l $l2
    }

}

/* ----------------------------------------------------------------------*/
/* */
/*  ::math::ln_Gamma --*/
/* */
/*  Returns ln(Gamma(x)), where x >= 0*/
/* */
/*  Parameters:*/
/*  x - Argument to the Gamma function.*/
/* */
/*  Results:*/
/*  Returns the natural logarithm of Gamma(x).*/
/* */
/*  Side effects:*/
/*  None.*/
/* */
/*  Gamma(x) is defined as:*/
/* */
/*        +inf*/
/*          _*/
/*         |    x-1  -t*/
/*  Gamma(x)= _|   t    e   dt*/
/* */
/*         0*/
/* */
/*  The approximation used here is from Lanczos, SIAM J. Numerical Analysis,*/
/*  series B, volume 1, p. 86.  For x > 1, the absolute error of the*/
/*  result is claimed to be smaller than 5.5 * 10**-10 -- that is, the*/
/*  resulting value of Gamma when exp( ln_Gamma( x ) ) is computed is*/
/*  expected to be precise to better than nine significant figures.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::ln_Gamma ( type x ) {

    # Handle the common case of a real argument that's within the
    # permissible range.

    if { [string is double -strict $x]
     && ( $x > 0 )
     && ( $x <= 2.5563481638716906e+305 )
     } {
    set x [expr { $x - 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 ) }]
    } 

    # Handle the error cases.

    if { ![string is double -strict $x] } {
    return -code error [expectDouble $x]
    }

    if { $x <= 0.0 } {
    set proc [lindex [info level 0] 0]
    return -code error \
        -errorcode [list ARITH DOMAIN \
            "argument to $proc must be positive"] \
        "argument to $proc must be positive"
    }

    return -code error \
    -errorcode [list ARITH OVERFLOW \
            "floating-point value too large to represent"] \
    "floating-point value too large to represent"
    
}

/* ----------------------------------------------------------------------*/
/* */
/*  math::factorial --*/
/* */
/*  Returns the factorial of the argument x.*/
/* */
/*  Parameters:*/
/*  x -- Number whose factorial is to be computed.*/
/* */
/*  Results:*/
/*  Returns x!, the factorial of x.*/
/* */
/*  Side effects:*/
/*  None.*/
/* */
/*  For integer x, 0 <= x <= 12, an exact integer result is returned.*/
/* */
/*  For integer x, 13 <= x <= 21, an exact floating-point result is returned*/
/*  on machines with IEEE floating point.*/
/* */
/*  For integer x, 22 <= x <= 170, the result is exact to 1 ULP.*/
/* */
/*  For real x, x >= 0, the result is approximated by computing*/
/*  Gamma(x+1) using the ::math::ln_Gamma function, and the result is*/
/*  expected to be precise to better than nine significant figures.*/
/* */
/*  It is an error to present x <= -1 or x > 170, or a value of x that*/
/*  is not numeric.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::factorial ( type x ) {

    variable factorialList

    # Common case: factorial of a small integer

    if { [string is integer -strict $x]
     && $x >= 0
     && $x < [llength $factorialList] } {
    return [lindex $factorialList $x]
    }

    # Error case: not a number

    if { ![string is double -strict $x] } {
    return -code error [expectDouble $x]
    }

    # Error case: gamma in the left half plane

    if { $x <= -1.0 } {
    set proc [lindex [info level 0] 0]
    set message "argument to $proc must be greater than -1.0"
    return -code error -errorcode [list ARITH DOMAIN $message] $message
    }

    # Error case - gamma fails

    if { [catch { expr {exp( [ln_Gamma [expr { $x + 1 }]] )} } result] } {
    return -code error -errorcode $::errorCode $result
    }

    # Success - computed factorial n as Gamma(n+1)

    return $result

}

/* ----------------------------------------------------------------------*/
/* */
/*  ::math::choose --*/
/* */
/*  Returns the binomial coefficient C(n,k) = n!/k!(n-k)!*/
/* */
/*  Parameters:*/
/*  n -- Number of objects in the sampling pool*/
/*  k -- Number of objects to be chosen.*/
/* */
/*  Results:*/
/*  Returns C(n,k).  */
/* */
/*  Side effects:*/
/*  None.*/
/* */
/*  Results are expected to be accurate to ten significant figures.*/
/*  If both parameters are integers and the result fits in 32 bits, */
/*  the result is rounded to an integer.*/
/* */
/*  Integer results are exact up to at least n = 34.*/
/*  Floating point results are precise to better than nine significant */
/*  figures.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::choose ( type n , type k ) {

    variable pascal

    # Use a precomputed table for small integer args

    if { [string is integer -strict $n]
     && $n >= 0 && $n < 34
     && [string is integer -strict $k]
     && $k >= 0 && $k <= $n } {

    set i [expr { ( ( $n * ($n + 1) ) / 2 ) + $k }]

    return [lindex $pascal $i]

    }

    # Test bogus arguments

    if { ![string is double -strict $n] } {
    return -code error [expectDouble $n]
    }
    if { ![string is double -strict $k] } {
    return -code error [expectDouble $k]
    }

    # Forbid negative n

    if { $n < 0. } {
    set proc [lindex [info level 0] 0]
    set msg "first argument to $proc must be non-negative"
    return -code error -errorcode [list ARITH DOMAIN $msg] $msg
    }

    # Handle k out of range

    if { [string is integer -strict $k] && [string is integer -strict $n]
     && ( $k < 0 || $k > $n ) } {
    return 0
    }

    if { $k < 0. } {
    set proc [lindex [info level 0] 0]
    set msg "second argument to $proc must be non-negative,\
                 or both must be integers"
    return -code error -errorcode [list ARITH DOMAIN $msg] $msg
    }

    # Compute the logarithm of the desired binomial coefficient.

    if { [catch { expr { [ln_Gamma [expr { $n + 1 }]]
             - [ln_Gamma [expr { $k + 1 }]]
             - [ln_Gamma [expr { $n - $k + 1 }]] } } r] } {
    return -code error -errorcode $::errorCode $r
    }

    # Compute the binomial coefficient itself

    if { [catch { expr { exp( $r ) } } r] } {
    return -code error -errorcode $::errorCode $r
    }

    # Round to integer if both args are integers and the result fits

    if { $r <= 2147483647.5 
           && [string is integer -strict $n]
           && [string is integer -strict $k] } {
    return [expr { round( $r ) }]
    }

    return $r

}

/* ----------------------------------------------------------------------*/
/* */
/*  ::math::Beta --*/
/* */
/*  Return the value of the Beta function of parameters z and w.*/
/* */
/*  Parameters:*/
/*  z, w : Two real parameters to the Beta function*/
/* */
/*  Results:*/
/*  Returns the value of the Beta function.*/
/* */
/*  Side effects:*/
/*  None.*/
/* */
/*  Beta( w, z ) is defined as:*/
/* */
/*                1_*/
/*                |  (z-1)     (w-1)*/
/*  Beta( w, z ) = Beta( z, w ) =     | t     (1-t)      dt*/
/*               _|*/
/*                0*/
/* */
/*         = Gamma( z ) Gamma( w ) / Gamma( z + w )*/
/* */
/*  Results are returned as a floating point number precise to better*/
/*  than nine significant figures for w, z > 1.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::Beta ( type z , type w ) {

    # Check form of both args so that domain check can be made

    if { ![string is double -strict $z] } {
    return -code error [expectDouble $z]
    }
    if { ![string is double -strict $w] } {
    return -code error [expectDouble $w]
    }

    # Check sign of both args

    if { $z <= 0.0 } {
    set proc [lindex [info level 0] 0]
    set msg "first argument to $proc must be positive"
    return -code error -errorcode [list ARITH DOMAIN $msg] $msg
    }
    if { $w <= 0.0 } {
    set proc [lindex [info level 0] 0]
    set msg "second argument to $proc must be positive"
    return -code error -errorcode [list ARITH DOMAIN $msg] $msg
    }

    # Compute beta using gamma function, keeping stack trace clean.

    if { [catch { expr { exp( [ln_Gamma $z] + [ln_Gamma $w]
                  - [ln_Gamma [ expr { $z + $w }]] ) } } beta] } {

    return -code error -errorcode $::errorCode $beta

    } 

    return $beta

}

/* ----------------------------------------------------------------------*/
/* */
/*  Initialization of this file:*/
/* */
/*  Initialize the precomputed tables of factorials and binomial*/
/*  coefficients.*/
/* */
/* ----------------------------------------------------------------------*/

namespace ::math {
    InitializeFactorial
    InitializePascal
}