private static double betaSmallBLargeASeries()

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

    private static double betaSmallBLargeASeries(double a, double b, double x, double y, double s0, double mult,
            Policy pol, boolean normalised) {
        //
        // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
        //
        // Some values we'll need later, these are Eq 9.1:
        //
        final double bm1 = b - 1;
        final double t = a + bm1 / 2;
        double lx;
        double u;
        if (y < 0.35) {
            lx = Math.log1p(-y);
        } else {
            lx = Math.log(x);
        }
        u = -t * lx;
        // and from from 9.2:
        final double h = BoostGamma.regularisedGammaPrefix(b, u);
        if (h <= Double.MIN_NORMAL) {
            // Update:
            // Boost returns s0.
            // If this is zero then compute an expected sub-normal value
            // using the classic continued fraction representation.
            if (s0 == 0) {
                return ibetaFraction(a, b, x, y, pol, normalised);
            }
            return s0;
        }
        double prefix;
        if (normalised) {
            prefix = h / GammaRatio.delta(a, b);
            prefix /= Math.pow(t, b);
        } else {
            prefix = BoostGamma.fullIgammaPrefix(b, u) / Math.pow(t, b);
        }
        prefix *= mult;
        //
        // now we need the quantity Pn, unfortunately this is computed
        // recursively, and requires a full history of all the previous values
        // so no choice but to declare a big table and hope it's big enough...
        //
        final double[] p = new double[PN_SIZE];
        // see 9.3.
        p[0] = 1;
        //
        // Now an initial value for J, see 9.6:
        //
        double j = BoostGamma.gammaQ(b, u, pol) / h;
        //
        // Now we can start to pull things together and evaluate the sum in Eq 9:
        //
        // Value at N = 0
        double sum = s0 + prefix * j;
        // some variables we'll need:
        // 2*N+1
        int tnp1 = 1;
        double lx2 = lx / 2;
        lx2 *= lx2;
        double lxp = 1;
        final double t4 = 4 * t * t;
        double b2n = b;

        for (int n = 1; n < PN_SIZE; ++n) {
            //
            // begin by evaluating the next Pn from Eq 9.4:
            //
            tnp1 += 2;
            p[n] = 0;
            int tmp1 = 3;
            for (int m = 1; m < n; ++m) {
                final double mbn = m * b - n;
                p[n] += mbn * p[n - m] / BoostGamma.uncheckedFactorial(tmp1);
                tmp1 += 2;
            }
            p[n] /= n;
            p[n] += bm1 / BoostGamma.uncheckedFactorial(tnp1);
            //
            // Now we want Jn from Jn-1 using Eq 9.6:
            //
            j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
            lxp *= lx2;
            b2n += 2;
            //
            // pull it together with Eq 9:
            //
            final double r = prefix * p[n] * j;
            final double previous = sum;
            sum += r;
            // Update:
            // This continues until convergence at machine epsilon
            // |r| < eps * |sum|; r < 1
            // |r| / eps < |sum|; r > 1
            if (sum == previous) {
                break;
            }
        }
        return sum;
    }