## Incomplete gamma

EKB I've implemented the incomplete gamma function (and submitted it to tcllib). Here's the code:

```package require math
package require math::statistics

# Adapted from Fortran code in the Royal Statistical Society's StatLib
# library (http://lib.stat.cmu.edu/apstat/), algorithm AS 32 (with
# some modifications from AS 239)
#
# Calculate normalized incomplete gamma function
#
#                 1       / x               p-1
#   P(p,x) =  --------   |   dt exp(-t) * t
#             Gamma(p)  / 0
#
# Tested some values against R's pgamma function

proc incompleteGamma {x p {tol 1.0e-9}} {
set overflow 1.0e37

if {\$x < 0} {
error "x must be positive"
return
}
if {\$p <= 0} {
error "p must be greater than or equal to zero"
return
}

# If x is zero, incGamma is zero
if {\$x == 0.0} {
return 0.0
}

# Use normal approx is p > 1000
if {\$p > 1000} {
set pn1 [expr {3.0 * sqrt(\$p) * (pow(1.0 * \$x/\$p, 1.0/3.0) + 1.0/(9.0 * \$p) - 1.0)}]
# pnorm is not robust enough for this calculation (overflows); cdf-normal could also be used
return [::math::statistics::pnorm_quicker \$pn1]
}

# If x is extremely large compared to a (and now know p < 1000), then return 1.0
if {\$x > 1.e8} {
return 1.0
}

set factor [expr {exp(\$p * log(\$x) -\$x - [::math::ln_Gamma \$p])}]

# Use series expansion (first option) or continued fraction
if {\$x <= 1.0 || \$x < \$p} {
set gin 1.0
set term 1.0
set rn \$p
while {1} {
set rn [expr {\$rn + 1.0}]
set term [expr {1.0 * \$term * \$x/\$rn}]
set gin [expr {\$gin + \$term}]
if {\$term < \$tol} {
set gin [expr {1.0 * \$gin * \$factor/\$p}]
break
}
}
} else {
set a [expr {1.0 - \$p}]
set b [expr {\$a + \$x + 1.0}]
set term 0.0
set pn1 1.0
set pn2 \$x
set pn3 [expr {\$x + 1.0}]
set pn4 [expr {\$x * \$b}]
set gin [expr {1.0 * \$pn3/\$pn4}]
while {1} {
set a [expr {\$a + 1.0}]
set b [expr {\$b + 2.0}]
set term [expr {\$term + 1.0}]
set an [expr {\$a * \$term}]
set pn5 [expr {\$b * \$pn3 - \$an * \$pn1}]
set pn6 [expr {\$b * \$pn4 - \$an * \$pn2}]
if {\$pn6 != 0.0} {
set rn [expr {1.0 * \$pn5/\$pn6}]
set dif [expr {abs(\$gin - \$rn)}]
if {\$dif <= \$tol && \$dif <= \$tol * \$rn} {
break
}
set gin \$rn
}
set pn1 \$pn3
set pn2 \$pn4
set pn3 \$pn5
set pn4 \$pn6
# Too big? Rescale
if {abs(\$pn5) >= \$overflow} {
set pn1 [expr {\$pn1 / \$overflow}]
set pn2 [expr {\$pn2 / \$overflow}]
set pn3 [expr {\$pn3 / \$overflow}]
set pn4 [expr {\$pn4 / \$overflow}]
}
}
set gin [expr {1.0 - \$factor * \$gin}]
}

return \$gin

}```

 Category Mathematics