picomath

# Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.

# Note that the functions Gamma and LogGamma are mutually dependent.

proc gamma {x} {
    if {$x <= 0} {
        error "Invalid input" $x
    }

    # Split the function domain into three intervals:
    # (0, 0.001), [0.001, 12), and (12, infinity)

    ###########################################################################
    # First interval: (0, 0.001)
    #
    # For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
    # So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
    # The relative error over this interval is less than 6e-7.

    # Euler's gamma constant
    set gamma 0.577215664901532860606512090

    if {$x < 0.001} {
        return [expr 1.0/($x*(1.0 + $gamma*$x))]
    }

    ###########################################################################
    # Second interval: [0.001, 12)

    if {$x < 12.0} {
        # The algorithm directly approximates gamma over (1,2) and uses
        # reduction identities to reduce other arguments to this interval.
        
        set y $x
        set n 0
        set arg_was_less_than_one [expr $y < 1.0]

        # Add or subtract integers as necessary to bring y into (1,2)
        # Will correct for this below
        if {$arg_was_less_than_one} {
            set y [expr $y + 1.0]
        } else {
            # will use n later
            set n [expr entier($y) - 1]
            set y [expr $y - $n]
        }

        # numerator coefficients for approximation over the interval (1,2)
        set p [list \
            -1.71618513886549492533811E+0 \
             2.47656508055759199108314E+1 \
            -3.79804256470945635097577E+2 \
             6.29331155312818442661052E+2 \
             8.66966202790413211295064E+2 \
            -3.14512729688483675254357E+4 \
            -3.61444134186911729807069E+4 \
             6.64561438202405440627855E+4 \
        ]

        # denominator coefficients for approximation over the interval (1,2)
        set q [list \
            -3.08402300119738975254353E+1 \
             3.15350626979604161529144E+2 \
            -1.01515636749021914166146E+3 \
            -3.10777167157231109440444E+3 \
             2.25381184209801510330112E+4 \
             4.75584627752788110767815E+3 \
            -1.34659959864969306392456E+5 \
            -1.15132259675553483497211E+5 \
        ]

        set num 0.0
        set den 1.0

        set z [expr $y - 1]
        for {set i 0} {$i < 8} {incr i} {
            set num [expr ($num + [lindex $p $i])*$z]
            set den [expr $den*$z + [lindex $q $i]]
        }
        set result [expr $num/$den + 1.0]

        # Apply correction if argument was not initially in (1,2)
        if {$arg_was_less_than_one} {
            # Use identity gamma(z) = gamma(z+1)/z
            # The variable "result" now holds gamma of the original y + 1
            # Thus we use y-1 to get back the orginal y.
            set result [expr $result / ($y-1.0)]
        } else {
            # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
            for {set i 0} {$i < $n} {incr i} {
                set result [expr $result * $y]
                set y [expr $y + 1]
            }
        }

        return $result
    }

    ###########################################################################
    # Third interval: [12, infinity)

    if {$x > 171.624} {
        # Correct answer too large to display. 
        return 1.0/0 # float infinity
    }

    return [expr exp([log_gamma $x])]
}

proc log_gamma {x} {
    if {$x <= 0} {
        error "Invalid input"
    }

    if {$x < 12.0} {
        return [expr log(abs([gamma $x]))]
    }

    # Abramowitz and Stegun 6.1.41
    # Asymptotic series should be good to at least 11 or 12 figures
    # For error analysis, see Whittiker and Watson
    # A Course in Modern Analysis (1927), page 252

    set c [list \
         1.0/12.0 \
        -1.0/360.0 \
         1.0/1260.0 \
        -1.0/1680.0 \
         1.0/1188.0 \
        -691.0/360360.0 \
         1.0/156.0 \
        -3617.0/122400.0 \
    ]
    set z [expr 1.0/($x*$x)]
    set sum [lindex $c 7]
    for {set i 6} {$i >= 0} {incr i -1} {
        set sum [expr $sum * $z]
        set sum [expr $sum + [lindex $c $i]]
    }
    set series [expr $sum/$x]

    set halfLogTwoPi 0.91893853320467274178032973640562
    set logGamma [expr ($x - 0.5)*log($x) - $x + $halfLogTwoPi + $series]
    return $logGamma
}