private static double log1pmxSmall()

in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/SpecialMath.java [157:218]


    private static double log1pmxSmall(double x, double a) {
        // Use a direct Taylor series:
        // ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
        // Reverse the summation (small to large) for a marginal increase in precision.
        // To stop the Taylor series the next term must be less than 1 ulp from the
        // answer.
        // x^n/n < |log(1+x)-x| * eps
        // eps = machine epsilon = 2^-53
        // x^n < |log(1+x)-x| * eps
        // n < (log(|log(1+x)-x|) + log(eps)) / log(x)
        // In practice this is a conservative limit.

        final double x2 = x * x;

        if (a < TWO_POW_M53) {
            // Below machine epsilon. Addition of x^3/3 is not possible.
            // Subtract from zero to prevent creating -0.0 for x=0.
            return 0 - x2 / 2;
        }

        final double x4 = x2 * x2;

        // +/-9.5367431640625e-07: log1pmx = -4.547470617660916e-13 :
        // -4.5474764000725028e-13
        // n = 4.69
        if (a < TWO_POW_M20) {
            // n=5
            return x * x4 / 5 -
                       x4 / 4 +
                   x * x2 / 3 -
                       x2 / 2;
        }

        // +/-2.44140625E-4: log1pmx = -2.9797472637290841e-08 : -2.9807173914456693e-08
        // n = 6.49
        if (a < TWO_POW_M12) {
            // n=7
            return x * x2 * x4 / 7 -
                       x2 * x4 / 6 +
                        x * x4 / 5 -
                            x4 / 4 +
                        x * x2 / 3 -
                            x2 / 2;
        }

        // Assume |x| < 2^-6
        // +/-0.015625: log1pmx = -0.00012081346403474586 : -0.00012335696813916864
        // n = 10.9974

        // n=11
        final double x8 = x4 * x4;
        return x * x2 * x8 / 11 -
                   x2 * x8 / 10 +
                    x * x8 /  9 -
                        x8 /  8 +
               x * x2 * x4 /  7 -
                   x2 * x4 /  6 +
                    x * x4 /  5 -
                        x4 /  4 +
                    x * x2 /  3 -
                        x2 /  2;
    }