public RealMatrix power(final int p)
        throws NotPositiveException, NonSquareMatrixException {
        if (p < 0) {
            throw new NotPositiveException(LocalizedFormats.NOT_POSITIVE_EXPONENT, p);
        }

        if (!isSquare()) {
            throw new NonSquareMatrixException(getRowDimension(), getColumnDimension());
        }

        if (p == 0) {
            return MatrixUtils.createRealIdentityMatrix(this.getRowDimension());
        }

        if (p == 1) {
            return this.copy();
        }

        final int power = p - 1;

        /*
         * Only log_2(p) operations is used by doing as follows:
         * 5^214 = 5^128 * 5^64 * 5^16 * 5^4 * 5^2
         *
         * In general, the same approach is used for A^p.
         */

        final char[] binaryRepresentation = Integer.toBinaryString(power).toCharArray();
        final ArrayList<Integer> nonZeroPositions = new ArrayList<Integer>();
        int maxI = -1;

        for (int i = 0; i < binaryRepresentation.length; ++i) {
            if (binaryRepresentation[i] == '1') {
                final int pos = binaryRepresentation.length - i - 1;
                nonZeroPositions.add(pos);

                // The positions are taken in turn, so maxI is only changed once
                if (maxI == -1) {
                    maxI = pos;
                }
            }
        }

        RealMatrix[] results = new RealMatrix[maxI + 1];
        results[0] = this.copy();

        for (int i = 1; i <= maxI; ++i) {
            results[i] = results[i-1].multiply(results[i-1]);
        }

        RealMatrix result = this.copy();

        for (Integer i : nonZeroPositions) {
            result = result.multiply(results[i]);
        }

        return result;
    }