static double beta()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostBeta.java [108:158]
38 lines of code
14 McCabe index (conditional complexity)

    static double beta(double p, double q) {
        if (!(p > 0 && q > 0)) {
            // Domain error
            return Double.NaN;
        }

        final double c = p + q;

        // Special cases:
        if (c == p && q < EPSILON) {
            return 1 / q;
        } else if (c == q && p < EPSILON) {
            return 1 / p;
        }
        if (q == 1) {
            return 1 / p;
        } else if (p == 1) {
            return 1 / q;
        } else if (c < EPSILON) {
            return (c / p) / q;
        }

        // Create input a > b
        final double a = p < q ? q : p;
        final double b = p < q ? p : q;

        // Lanczos calculation:
        final double agh = a + BoostGamma.Lanczos.GMH;
        final double bgh = b + BoostGamma.Lanczos.GMH;
        final double cgh = c + BoostGamma.Lanczos.GMH;
        double result = BoostGamma.Lanczos.lanczosSumExpGScaled(a) *
                (BoostGamma.Lanczos.lanczosSumExpGScaled(b) / BoostGamma.Lanczos.lanczosSumExpGScaled(c));
        final double ambh = a - 0.5f - b;
        if (Math.abs(b * ambh) < cgh * 100 && a > 100) {
            // Special case where the base of the power term is close to 1
            // compute (1+x)^y instead:
            result *= Math.exp(ambh * Math.log1p(-b / cgh));
        } else {
            result *= Math.pow(agh / cgh, ambh);
        }

        if (cgh > 1e10f) {
            // this avoids possible overflow, but appears to be marginally less accurate:
            result *= Math.pow((agh / cgh) * (bgh / cgh), b);
        } else {
            result *= Math.pow((agh * bgh) / (cgh * cgh), b);
        }
        result *= Math.sqrt(Math.E / bgh);

        return result;
    }