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
    
}