function Z = mpower(X,Y)
%^   Matrix power.
%   Z = X^y is X to the y power if y is a scalar and X is square. If y
%   is an integer greater than one, the power is computed by repeated
%   squaring. For other values of y the calculation involves
%   eigenvalues and eigenvectors.
%
%   Z = x^Y is x to the Y power if Y is a square matrix and x is a scalar.
%   Computed using eigenvalues and eigenvectors.
%
%   Z = X^Y, where both X and Y are matrices, is an error.
%
%   C = MPOWER(A,B) is called for the syntax 'A ^ B' when A or B is an
%   object.
%
%   See also POWER.

%   Copyright 1984-2018 The MathWorks, Inc.

if isscalar(X) && isscalar(Y)
    Z = X.^Y;
else
    if isinteger(X) || isinteger(Y)
        error(message('MATLAB:mpower:integerNotSupported'));
    end

    if ~ismatrix(X) || ~ismatrix(Y)
        error(message('MATLAB:mpower:inputsMustBe2D'));
    end

    if isa(X,'single') || isa(Y, 'single')
        try
            X = single(X);
            Y = single(Y);
        catch err
            if (strcmp(err.identifier, 'MATLAB:unimplementedSparseType'))
                error(message('MATLAB:mpower:sparseSingleNotSupported'))
            else
                rethrow(err);
            end
        end
    else
        X = double(X);
        Y = double(Y);
    end

    if isscalar(Y) && isreal(Y) && isfinite(Y) && round(Y) == Y ... 
            && size(X,1) == size(X,2)
        Z = integerMpower(X,Y);
    else
        Z = generalMpower(X,Y);
    end
end

function Z = integerMpower(X,y)
if y == 0
    Z = eye(size(X),'like',X([])); %always real
else
    % X and y can be sparse
    % Z = X^y for integer y. Use repeated squaring.
    % For example: A^13 = A^1 * A^4 * A^8
    p = abs(y);
    D = X;
    first = true;
    while p > 0
        if rem(p,2) == 1 %if odd
            if first
                Z = D;  % hit first time. D*I
                first = false;
            else
                Z = D*Z;
            end
        end
        p = fix(p/2);
        if p ~= 0
            D = D*D;
        end
    end
    if y < 0
        Z = inv(Z);
    end
end


function Z = generalMpower(X,Y)
if isscalar(Y) && size(X,1) == size(X,2)
    % If y is a not a real integer-valued scalar, use eigenvalues
    % X and y should never be sparse.
    try
        [V,d] = eig(full(X),'vector');
        Z = (V.*(d.^Y).')/V;
    catch err
        if (strcmp(err.identifier,'MATLAB:eig:matrixWithNaNInf'))
            Z = NaN(size(X),superiorfloat(X,Y));
        else
            rethrow(err);
        end
    end
    if issparse(X)
       Z = sparse(Z);
    end
elseif isscalar(X) && size(Y,1) == size(Y,2)
    % If x is a scalar and Y is a matrix, use eigenvalues.
    % x and Y should never be sparse.
    try
        [V,d] = eig(full(Y),'vector');
        Z = (V.*(X.^d).')/V;
    catch err
        if (strcmp(err.identifier,'MATLAB:eig:matrixWithNaNInf'))
            Z = NaN(size(Y),superiorfloat(X,Y));
        else
            rethrow(err);
        end
    end
    if issparse(X)
       Z = sparse(Z);
    end
else
    Me = MException(message('MATLAB:mpower:notScalarAndSquareMatrix'));
    throwAsCaller(Me);
end
