static double tgammaRatio()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java [2215:2267]
37 lines of code
14 McCabe index (conditional complexity)

    static double tgammaRatio(double x, double y) {
        if (x <= 0 || !Double.isFinite(x) || y <= 0 || !Double.isFinite(y)) {
            return Double.NaN;
        }
        if (x <= Double.MIN_NORMAL) {
            // Special case for denorms...Ugh.
            return TWO_POW_53 * tgammaRatio(x * TWO_POW_53, y);
        }

        if (x <= MAX_GAMMA_Z && y <= MAX_GAMMA_Z) {
            // Rather than subtracting values, lets just call the gamma functions directly:
            return tgamma(x) / tgamma(y);
        }
        double prefix = 1;
        if (x < 1) {
            if (y < 2 * MAX_GAMMA_Z) {
                // We need to sidestep on x as well, otherwise we'll underflow
                // before we get to factor in the prefix term:
                prefix /= x;
                x += 1;
                while (y >= MAX_GAMMA_Z) {
                    y -= 1;
                    prefix /= y;
                }
                return prefix * tgamma(x) / tgamma(y);
            }
            //
            // result is almost certainly going to underflow to zero, try logs just in case:
            //
            return Math.exp(lgamma(x) - lgamma(y));
        }
        if (y < 1) {
            if (x < 2 * MAX_GAMMA_Z) {
                // We need to sidestep on y as well, otherwise we'll overflow
                // before we get to factor in the prefix term:
                prefix *= y;
                y += 1;
                while (x >= MAX_GAMMA_Z) {
                    x -= 1;
                    prefix *= x;
                }
                return prefix * tgamma(x) / tgamma(y);
            }
            //
            // Result will almost certainly overflow, try logs just in case:
            //
            return Math.exp(lgamma(x) - lgamma(y));
        }
        //
        // Regular case, x and y both large and similar in magnitude:
        //
        return tgammaDeltaRatio(x, y - x);
    }