private static double ibetaSeries()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostBeta.java [835:913]
52 lines of code
8 McCabe index (conditional complexity)

    private static double ibetaSeries(double a, double b, double x, double s0, boolean normalised, Policy pol) {
        double result;

        if (normalised) {
            final double c = a + b;

            // incomplete beta power term, combined with the Lanczos approximation:
            final double agh = a + BoostGamma.Lanczos.GMH;
            final double bgh = b + BoostGamma.Lanczos.GMH;
            final double cgh = c + BoostGamma.Lanczos.GMH;
            result = BoostGamma.Lanczos.lanczosSumExpGScaled(c) /
                    (BoostGamma.Lanczos.lanczosSumExpGScaled(a) * BoostGamma.Lanczos.lanczosSumExpGScaled(b));

            final double l1 = Math.log(cgh / bgh) * (b - 0.5f);
            final double l2 = Math.log(x * cgh / agh) * a;
            //
            // Check for over/underflow in the power terms:
            //
            if (l1 > LOG_MIN_VALUE && l1 < LOG_MAX_VALUE && l2 > LOG_MIN_VALUE && l2 < LOG_MAX_VALUE) {
                if (a * b < bgh * 10) {
                    result *= Math.exp((b - 0.5f) * Math.log1p(a / bgh));
                } else {
                    result *= Math.pow(cgh / bgh, b - 0.5f);
                }
                result *= Math.pow(x * cgh / agh, a);
                result *= Math.sqrt(agh / Math.E);
            } else {
                //
                // Oh dear, we need logs, and this *will* cancel:
                //
                result = Math.log(result) + l1 + l2 + (Math.log(agh) - 1) / 2;
                result = Math.exp(result);
            }
        } else {
            // Non-normalised, just compute the power:
            result = Math.pow(x, a);
        }
        double rescale = 1.0;
        if (result < Double.MIN_NORMAL) {
            // Safeguard: series can't cope with denorms.

            // Update:
            // The entire series is only based on the magnitude of 'result'.
            // If the first term can be added to s0 (e.g. if s0 == 0) then
            // scale s0 and result, compute the series and then rescale.

            // Intentional floating-point equality check.
            if (s0 + result / a == s0) {
                return s0;
            }
            s0 *= TWO_POW_53;
            result *= TWO_POW_53;
            rescale = TWO_POW_M53;
        }

        final double eps = pol.getEps();
        final int maxIterations = pol.getMaxIterations();

        // Create effectively final 'result' for initialisation
        final double result1 = result;
        final DoubleSupplier gen = new DoubleSupplier() {
            /** Next result term. */
            private double result = result1;
            /** Pochhammer term. */
            private final double poch = -b;
            /** Iteration. */
            private int n;

            @Override
            public double getAsDouble() {
                final double r = result / (a + n);
                n++;
                result *= (n + poch) * x / n;
                return r;
            }
        };

        return BoostTools.sumSeries(gen, eps, maxIterations, s0) * rescale;
    }