Gamma Function

In mathematics, the gamma function is an extension of the factorial function, with its argument shifted down by 1, to real and complex numbers.

Gamma FunctionThe gamma function is a very common special function, which finds application in such diverse areas as quantum physics, astrophysics and fluid dynamics.

The gamma distribution, which is formulated in terms of the gamma function, is used in statistics to model a wide range of processes.

The source code presents an accurate estimation of the gamma function and factorials both on the absolute and the logarithmic scale.

The arguments can be all real numbers except negative integers. The source code uses Euler’s reflection formula to estimate gamma for negative numbers.

#COMPILE EXE
#DIM ALL

FUNCTION GammLn(BYVAL x AS EXT) AS EXT

    ' Returns Ln(Gamma()) or 0 on error
    ' Based on Numerical Recipes gamma.h
    ' Lanczos, C. 1964, "A Precision Approximation
    ' of the Gamma Function," SIAM Journal on Numerical
    ' Analysis, ser. B, vol. 1, pp. 86-96.

    LOCAL j AS LONG, tmp, y, ser AS EXT
    DIM cof(0 TO 13) AS LOCAL EXT
    cof(0) = 57.1562356658629235
    cof(1) = -59.5979603554754912
    cof(2) = 14.1360979747417471
    cof(3) = -0.491913816097620199
    cof(4) = 0.339946499848118887e-4
    cof(5) = 0.465236289270485756e-4
    cof(6) = -0.983744753048795646e-4
    cof(7) = 0.158088703224912494e-3
    cof(8) = -0.210264441724104883e-3
    cof(9) = 0.217439618115212643e-3
    cof(10) = -0.164318106536763890e-3
    cof(11) = 0.844182239838527433e-4
    cof(12) = -0.261908384015814087e-4
    cof(13) = 0.368991826595316234e-5
    IF x <= 0.0 THEN FUNCTION = 0.0 : EXIT FUNCTION ' Bad argument
    y = x
    tmp = x + 5.2421875
    tmp = (x + 0.5) * LOG(tmp) - tmp
    ser = 0.999999999999997092
    FOR j = 0 TO 13
        y = y + 1
        ser = ser + cof(j)/y
    NEXT j
    FUNCTION = tmp + LOG(2.5066282746310005 * ser / x)
END FUNCTION
'
FUNCTION PBMAIN () AS LONG
    LOCAL s1, s2, s3, s4, X, X1, PI AS EXT, t, t2 AS STRING
    PI = 4.0 * ATN(1.0)
    t = INPUTBOX$("Enter argument X", "Get Gamma Results",)
    X = VAL(t)
    IF X < -1.0 AND ROUND(X, 0) <> X THEN ' X < -1 and non-integer
        X1 = ABS(X) + 1.0
        s1 = GammLn(X1)
        s3 = PI / (EXP(s1) * SIN(PI * X1)) ' Gamma
        IF s3 > 0.0 THEN s1 = LOG(s3) ELSE s1 = 1000.0 ' Ln of Gamma
        X1 = X1 - 1.0
        s2 = GammLn(X1)
        s4 = PI / (EXP(s2) * SIN(PI * X1)) ' Factorial
        IF s4 > 0.0 THEN s2 = LOG(s4) ELSE s2 = 1000.0 ' Ln of Factorial
    END IF
    IF X > -1.0 AND X < 0.0 THEN
        X1 = ABS(X) + 1.0
        s1 = GammLn(X1)
        s3 = PI / (EXP(s1) * SIN(PI * X1)) ' Gamma
        IF s3 > 0.0 THEN s1 = LOG(s3) ELSE s1 = 1000.0 ' Ln of Gamma
        X1 = X1 + 1.0
        s2 = GammLn(X1)                    ' Ln of Factorial
        s4 = EXP(s2)                       ' Factorial
    END IF
    IF X > 0.0 THEN
        s1 = GammLn(X)                    ' Ln of Gamma
        s2 = GammLn(X + 1.0)              ' Ln of Factorial
        IF VAL(t) <= 1753 AND VAL(t) >= -1753 THEN ' max factorial limit
            s3 = EXP(s1)
            s4 = EXP(s2)
        END IF
    END IF
    IF s3 <> 0.0 AND s4 <> 0.0 THEN
        t2 = "Gamma of " + t + " =" + _
        num(s3) + $CRLF + $CRLF
        t2 = t2 + "Factorial of " + t + " =" + _
        num(s4) + $CRLF + $CRLF
    ELSE
        t2 = t2 + "Gamma and Factorial: overflow." + _
        $CRLF + $CRLF
    END IF
    IF s1 <> 1000.0 THEN t2 = t2 + _
        "Ln of Gamma of " + t + " =" + num(s1) + $CRLF + $CRLF
    IF s2 <> 1000.0 THEN t2 = t2 + "Ln of Factorial of " + _
        t + " =" + num(s2)
    IF X <=-1.0 AND ROUND(X, 0) = X THEN ' X <= -1 and integer
        IF X = -1.0 THEN
            t2 = "Gamma not defined for X = -1" + _
            $CRLF + $CRLF + "Factorial of -1 = 1"
        ELSE
            t2 = "Gamma and Factorial not defined for X = " + t
        END IF
    END IF
    MSGBOX t2, , "Gamma Function Results:"
END FUNCTION
'
FUNCTION num(BYVAL v AS EXT) AS STRING
    IF ABS(v) < 10000.0 AND ABS(v) > 0.0001 THEN
        FUNCTION = FORMAT$(v, "#####.################")
    ELSE
        FUNCTION = FORMAT$(v," 0.################E-####")
    END IF
END FUNCTION

Gamma Function Program

Gamma Function CalculationEstimate Gamma Function for negative argument Gamma Results for negative argument