public RealSquareMatrix power()

in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/SquareMatrixSupport.java [130:189]

• 31 lines of code
• 5 McCabe index (conditional complexity)

        public RealSquareMatrix power(int n) {
            checkExponent(n);
            if (n == 0) {
                return identity();
            }
            if (n == 1) {
                return this;
            }

            // Here at least 1 multiplication occurs.
            // Compute the power by repeat squaring and multiplication:
            // 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
            // 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 we require 2 matrix array allocations and a column array.

            // Working arrays
            final double[] col = new double[dim];
            double[] b = new double[data.length];
            double[] tmp;

            // Initialise result as A^1.
            final double[] a = data;
            final int ea = exp;
            double[] r = a.clone();
            int er = ea;

            // 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;

            // Process remaining bits below highest set bit.
            for (int i = 32 - shift; i != 0; i--, bits <<= 1) {
                // Square the result
                er = multiply(r, er, r, er, col, b);
                // Recycle working array
                tmp = b;
                b = r;
                r = tmp;
                if (bits < 0) {
                    // Multiply by A
                    er = multiply(r, er, a, ea, col, b);
                    // Recycle working array
                    tmp = b;
                    b = r;
                    r = tmp;
                }
            }

            return new ArrayRealSquareMatrix(dim, r, er);
        }