#----------------------------------------------------------------------
#
# 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 eval ::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.
#
#----------------------------------------------------------------------

proc ::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.
#
#----------------------------------------------------------------------

proc ::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.
#
#----------------------------------------------------------------------

proc ::math::ln_Gamma { 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.
#
#----------------------------------------------------------------------

proc ::math::factorial { 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.
#
#----------------------------------------------------------------------

proc ::math::choose { n 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.
#
#----------------------------------------------------------------------

proc ::math::Beta { z 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 eval ::math {
    InitializeFactorial
    InitializePascal
}