static double tgamma()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java [595:690]
45 lines of code
11 McCabe index (conditional complexity)

    static double tgamma(double z) {
        // Handle integers
        if (Math.rint(z) == z) {
            if (z <= 0) {
                // Pole error
                return Double.NaN;
            }
            if (z <= MAX_GAMMA_Z) {
                // Gamma(n) = (n-1)!
                return FACTORIAL[(int) z - 1];
            }
            // Overflow
            return Double.POSITIVE_INFINITY;
        }

        if (Math.abs(z) <= LANCZOS_THRESHOLD) {
            // Small z
            // NSWC Library of Mathematics Subroutines
            // Note:
            // This does not benefit from using extended precision to track the sum (t).
            // Extended precision on the product reduces the error but the majority
            // of error remains in InvGamma1pm1.

            if (z >= 1) {
                /*
                 * From the recurrence relation
                 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
                 * then
                 * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
                 * where t = x - n. This means that t must satisfy
                 * -0.5 <= t - 1 <= 1.5.
                 */
                double prod = 1;
                double t = z;
                while (t > 2.5) {
                    t -= 1;
                    prod *= t;
                }
                return prod / (1 + InvGamma1pm1.value(t - 1));
            }
            /*
             * From the recurrence relation
             * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
             * then
             * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
             * which requires -0.5 <= x + n <= 1.5.
             */
            double prod = z;
            double t = z;
            while (t < -0.5) {
                t += 1;
                prod *= t;
            }
            return 1 / (prod * (1 + InvGamma1pm1.value(t)));
        }

        // Large non-integer z
        // Boost C++ tgamma implementation

        if (z < 0) {
            /*
             * From the reflection formula
             * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
             * and the recurrence relation
             * Gamma(1 - x) = -x * Gamma(-x),
             * it is found
             * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
             */
            return -Math.PI / (sinpx(z) * tgamma(-z));
        } else if (z > MAX_GAMMA_Z + 1) {
            // Addition to the Boost code: Simple overflow detection
            return Double.POSITIVE_INFINITY;
        }

        double result = Lanczos.lanczosSum(z);
        final double zgh = z + Lanczos.GMH;
        final double lzgh = Math.log(zgh);
        if (z * lzgh > LOG_MAX_VALUE) {
            // we're going to overflow unless this is done with care:

            // Updated
            // Check for overflow removed:
            // if (lzgh * z / 2 > LOG_MAX_VALUE) ... overflow
            // This is replaced by checking z > MAX_FACTORIAL + 2

            final double hp = Math.pow(zgh, (z / 2) - 0.25);
            result *= hp / Math.exp(zgh);
            // Check for overflow has been removed:
            // if (Double.MAX_VALUE / hp < result) ... overflow
            result *= hp;
        } else {
            result *= Math.pow(zgh, z - 0.5) / Math.exp(zgh);
        }

        return result;
    }