private static DD computePowScaled()

in commons-numbers-core/src/main/java/org/apache/commons/numbers/core/DD.java [2083:2163]

• 44 lines of code
• 4 McCabe index (conditional complexity)

    private static DD computePowScaled(long b, double x, double xx, int n, long[] exp) {
        // Compute the power by multiplication (keeping track of the scale):
        // 13 = 1101
        // x^13 = x^8 * x^4 * x^1
        //      = ((x^2 * x)^2)^2 * x
        // 21 = 10101
        // x^21 = x^16 * x^4 * x^1
        //      = (((x^2)^2 * x)^2)^2 * x
        // 1. Find highest set bit in n (assume n != 0)
        // 2. Initialise result as x
        // 3. For remaining bits (0 or 1) below the highest set bit:
        //    - square the current result
        //    - if the current bit is 1 then multiply by x
        // In this scheme the factors to multiply by x can be pre-computed.

        // Scale the input in [0.5, 1) to be above 1. Represented as 2^be * b.
        final long be = b - 1;
        final double b0 = x * 2;
        final double b1 = xx * 2;
        // Split b
        final double b0h = highPart(b0);
        final double b0l = b0 - b0h;

        // Initialise the result as x^1. Represented as 2^fe * f.
        long fe = be;
        double f0 = b0;
        double f1 = b1;

        double u;
        double v;
        double w;

        // Shift the highest set bit off the top.
        // Any remaining bits are detected in the sign bit.
        final int shift = Integer.numberOfLeadingZeros(n) + 1;
        int bits = n << shift;

        // Multiplication is done without using DD.multiply as the arguments
        // are always finite and the product will not overflow. The square can be optimised.
        // Process remaining bits below highest set bit.
        for (int i = 32 - shift; i != 0; i--, bits <<= 1) {
            // Square the result
            // Inline multiply(f0, f1, f0, f1, f), adapted for squaring
            fe <<= 1;
            u = f0 * f0;
            v = twoSquareLow(f0, u);
            // Inline fastTwoSum(hi, lo + (2 * f0 * f1), f)
            w = v + (2 * f0 * f1);
            f0 = u + w;
            f1 = fastTwoSumLow(u, w, f0);
            // Rescale
            if (Math.abs(f0) > SAFE_MULTIPLY) {
                // Scale back to the [1, 2) range. As safe multiply is 2^500
                // the exponent should be < 1001 so the twoPow scaling factor is supported.
                final int e = Math.getExponent(f0);
                final double s = twoPow(-e);
                fe += e;
                f0 *= s;
                f1 *= s;
            }
            if (bits < 0) {
                // Multiply by b
                fe += be;
                // Inline multiply(f0, f1, b0, b1, f)
                u = highPart(f0);
                v = f0 - u;
                w = f0 * b0;
                v = twoProductLow(u, v, b0h, b0l, w);
                // Inline fastTwoSum(w, v + (f0 * b1 + f1 * b0), f)
                u = v + (f0 * b1 + f1 * b0);
                f0 = w + u;
                f1 = fastTwoSumLow(w, u, f0);
                // Avoid rescale as x2 is in [1, 2)
            }
        }

        final int[] e = {0};
        final DD f = new DD(f0, f1).frexp(e);
        exp[0] = fe + e[0];
        return f;
    }