public static double value()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogBeta.java [155:247]
88 lines of code
19 McCabe index (conditional complexity)

    public static double value(double p,
                               double q) {
        if (Double.isNaN(p) ||
            Double.isNaN(q) ||
            p <= 0 ||
            q <= 0) {
            return Double.NaN;
        }

        final double a = Math.min(p, q);
        final double b = Math.max(p, q);
        if (a >= TEN) {
            final double w = sumDeltaMinusDeltaSum(a, b);
            final double h = a / b;
            final double c = h / (1 + h);
            final double u = -(a - 0.5) * Math.log(c);
            final double v = b * Math.log1p(h);
            if (u <= v) {
                return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
            }
            return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
        } else if (a > TWO) {
            if (b > THOUSAND) {
                final int n = (int) Math.floor(a - 1);
                double prod = 1;
                double ared = a;
                for (int i = 0; i < n; i++) {
                    ared -= 1;
                    prod *= ared / (1 + ared / b);
                }
                return (Math.log(prod) - n * Math.log(b)) +
                        (LogGamma.value(ared) +
                         logGammaMinusLogGammaSum(ared, b));
            }
            double prod1 = 1;
            double ared = a;
            while (ared > TWO) {
                ared -= 1;
                final double h = ared / b;
                prod1 *= h / (1 + h);
            }
            if (b < TEN) {
                double prod2 = 1;
                double bred = b;
                while (bred > TWO) {
                    bred -= 1;
                    prod2 *= bred / (ared + bred);
                }
                return Math.log(prod1) +
                       Math.log(prod2) +
                       (LogGamma.value(ared) +
                       (LogGamma.value(bred) -
                        LogGammaSum.value(ared, bred)));
            }
            return Math.log(prod1) +
                   LogGamma.value(ared) +
                   logGammaMinusLogGammaSum(ared, b);
        } else if (a >= 1) {
            if (b > TWO) {
                if (b < TEN) {
                    double prod = 1;
                    double bred = b;
                    while (bred > TWO) {
                        bred -= 1;
                        prod *= bred / (a + bred);
                    }
                    return Math.log(prod) +
                           (LogGamma.value(a) +
                            (LogGamma.value(bred) -
                             LogGammaSum.value(a, bred)));
                }
                return LogGamma.value(a) +
                    logGammaMinusLogGammaSum(a, b);
            }
            return LogGamma.value(a) +
                   LogGamma.value(b) -
                   LogGammaSum.value(a, b);
        } else {
            if (b >= TEN) {
                return LogGamma.value(a) +
                       logGammaMinusLogGammaSum(a, b);
            }
            // The original NSWC implementation was
            //   LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
            // but the following command turned out to be more accurate.
            // Note: Check for overflow that occurs if a and/or b are tiny.
            final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b);
            if (Double.isFinite(beta)) {
                return Math.log(beta);
            }
            return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
        }
    }