function y_out = mpower(a,kin)
%MATLAB Code Generation Library Function

%   Copyright 2007-2015 The MathWorks, Inc.
%#codegen
coder.extrinsic('emlGetNTypeForMpower');
coder.extrinsic('strcmpi');
coder.internal.prefer_const(a, kin);
eml_invariant(nargin == 2,'MATLAB:narginchk:notEnoughInputs');

eml_invariant(isfi(a),'fixed:fi:firstInputNotFi');

eml_invariant(~(isfi(a)&&isslopebiasscaled(numerictype(a))),...
              'fixed:fi:binaryPointOnlyMath','mpower');

eml_invariant(ismatrix(a)&&eml_is_const(ndims(a))&&isequal(size(a,1),size(a,2)),...
              'MATLAB:mpower:notScalarAndSquareMatrix');

if eml_is_const(kin)
    eml_invariant(isnumeric(kin)&&isreal(kin)&&isequal(kin, floor(kin))&&isscalar(kin)&&isfinite(kin)&&(kin >= 0),...
                  'fixed:fi:invalidExponent','mpower');
else
    if (isscaledtype(a))
        if isscalar(a)&&eml_is_const(isscalar(a))
            isProdModeSpecPrec = coder.const(strcmpi(a.ProductMode,'SpecifyPrecision'));
            eml_invariant(isProdModeSpecPrec,...
                          'fixed:coder:mpowerProductModeForNonConstExponent');       
        else 
            isSumModeSpecPrec = coder.const(strcmpi(a.SumMode,'SpecifyPrecision'));
            eml_invariant(isSumModeSpecPrec,...
                          'fixed:coder:mpowerSumModeForNonConstExponent');
        end
    end
    
    % check at run time that kin is a non-negative integer if kin is not const
    eml_invariant(isnumeric(kin)&&isreal(kin)&&isscalar(kin),...
                  'fixed:fi:invalidExponent', 'mpower');
    if ~(isfinite(kin)&&isequal(kin,floor(kin))&&(kin>=0))
        coder.internal.error('fixed:fi:invalidExponent','mpower'); % run-time error
    end    
end

if ~isfloat(a)&&~eml_is_const(size(a))
    sMode = coder.const(get(a, 'SumMode'));

    eml_invariant((strcmpi(sMode,'SpecifyPrecision')||strcmpi(sMode,'KeepLSB')), ...
                  'fixed:coder:sumModeRestrictedForVarS','mpower');

end

% Cast to int32 after error checking that it's an integer
k = int32(kin);

if isfloat(a)
    
    atemp = eml_cast(a, eml_fi_getDType(a));
    y = eml_cast(mpower(atemp, eml_cast(k,'double')), numerictype(a));
elseif isscalar(a)&&eml_is_const(isscalar(a))
    
    y = power(a, k);
elseif ( eml_is_const(k) && (k == 0) )
    % Special case a ^ 0
    ta = eml_typeof(a);
    tEye = numerictype(numerictype(false,1,0),'DataType',ta.DataType);
    if eml_fimathislocal(a)
        y = fi(eye(size(a,1)), tEye, coder.const(eml_fimath(a)));
    else
        y = fi(eye(size(a,1)), tEye);
    end
elseif ( eml_is_const(k) && (k == 1) )
    % Special case a ^ 1
    y = a;
elseif ( eml_is_const(k) && k == 2 )
    % Special case a ^ 2
    y = a * a;
elseif ( eml_is_const(k) && k == 3 )
    % Special case a ^ 3
    y = a * a * a;
elseif ( eml_is_const(k) && k == 4 )
    % Special case a ^ 4
    b = a * a;
    y = b * b;
elseif coder.const(strcmpi(a.ProductMode,'FullPrecision')) && ...
        coder.const(strcmpi(a.SumMode,'FullPrecision'))
    % a^k with repeated squaring if both sum and product are full precision
    y = compute_mpower_repeated_squaring(a, k);
else
    
    fma = eml_fimath(a);
    pmode = coder.const(get(fma,'ProductMode'));
    smode = coder.const(get(fma,'SumMode'));

    isanymodesp = coder.const(strcmpi(pmode,'SpecifyPrecision'))||coder.const(strcmpi(smode,'SpecifyPrecision'));
    arebothmodeskmsb = coder.const(strcmpi(pmode,'KeepMSB'))&&coder.const(strcmpi(smode,'KeepMSB'));
    isanymodekmsb = coder.const(strcmpi(pmode,'KeepMSB'))||coder.const(strcmpi(smode,'KeepMSB'));

    if coder.const(isanymodesp)
        
        y = compute_mpower_typec(a, k);
    elseif coder.const(isanymodekmsb)
        
        ispmodekmsb = coder.const(strcmpi(pmode,'KeepMSB'));
        issmodeklsb = coder.const(strcmpi(smode,'KeepLSB'));
        ispmodeklsb = coder.const(strcmpi(pmode,'KeepLSB'));
        ispmodefp = coder.const(strcmpi(pmode,'FullPrecision'));
        issmodefp = coder.const(strcmpi(smode,'FullPrecision'));
        isanymodefp = coder.const(ispmodefp||issmodefp);
        awl = get(a, 'WordLength');
        swl = get(fma, 'SumWordLength');
        pwl = get(fma, 'ProductWordLength');
        if coder.const(ispmodekmsb&&issmodeklsb)
            
            n = swl - pwl;
        else
            
            if coder.const(isreal(a))
                
                nsz = coder.const(ceil(log2(size(a,1))));
            else
                
                nsz = coder.const(ceil(log2(size(a,1)+1)));
            end
            if (isanymodefp||arebothmodeskmsb)
                
                n = nsz;
            elseif ispmodeklsb
                n = (pwl - swl) - awl + nsz;
            end
        end
        y = compute_mpower_typeb(a, k, n);
    else
        y = compute_mpower_typea(a, k);
    end
    
end
% Propagate input fimath to output
if isfimathlocal(a)
    y_out = setfimath(y, a.fimath);
else
    y_out = removefimath(y);
end


function y_out = compute_mpower_repeated_squaring(a_in, k)
% a^k with repeated squaring
coder.internal.prefer_const(k);
if isscaleddouble(a_in)
    % Use builtin a.^k for scaled double
    [~,ty] = init_mpower(a_in, k);
    y_out = eml_cast(double(a_in)^double(k), ty);
else
    [~, ty] = init_mpower(a_in, k);
    % Compute in pure integer (stripscaling) so fraction length
    % does not change during the iterations.  Set product and sum
    % word lengths to the output type so you don't exceed the
    % 128-bit limit during multiply.
    F = fimath('RoundingMethod', a_in.RoundingMethod, ...
        'OverflowAction', a_in.OverflowAction, ...
        'ProductMode', 'SpecifyPrecision', ...
        'ProductWordLength', ty.WordLength, ...
        'ProductFractionLength', 0, ...
        'SumMode', 'SpecifyPrecision', ...
        'SumWordLength', ty.WordLength, ...
        'SumFractionLength', 0);

    if isreal(a_in)
        y_prototype = fi([],ty.SignednessBool,ty.WordLength,0,'DataType',a_in.DataType,'fimath',F);
    else
        y_prototype = complex(fi(0,ty.SignednessBool,ty.WordLength,0,'DataType',a_in.DataType,'fimath',F));
    end
    % a.^k with repeated squaring
    y = zeros(size(a_in), 'like', y_prototype);
    a = zeros(size(a_in), 'like', y_prototype);
    a(:) = stripscaling(a_in);

    k = uint32(k);
    one = uint32(1);
    initialized = false;
    if coder.const(strcmp(coder.target(),'hdl'))
        % HDL throws an "unbounded loop structure" error for while
        % loops, so bounding this loop with a for loop. It also
        % throws the error if there is a break in the for loop, so
        % you have to let it go to completion.
        N = k;
        for n = 1:N
            % while k>0
            if k ~=0 
                if bitand(k, one)
                    % k is odd
                    if initialized
                        y(:) = y * a;
                    else
                        y(:) = a;
                        initialized = true;
                    end
                end
                k = bitsrl(k,1);
                if k ~= 0
                    a(:) = a * a;
                end
            end
        end
    else
        % C codegen can use the while loop
        while k>0
            if bitand(k, one)
                % k is odd
                if initialized
                    y(:) = y * a;
                else
                    y(:) = a;
                    initialized = true;
                end
            end
            k = bitsrl(k,1);
            if k ~= 0
                a(:) = a * a;
            end
        end
    end
    y_out = reinterpretcast(y,ty);
end

function y = compute_mpower_typea(a, k)%#codegen
coder.internal.prefer_const(k);
[fmA, ty, ega] = init_mpower(a, k);
ty1 = numerictype(ty, 'FractionLength', 0);
bzero = eml_cast(ega, ty1, fmA);
b = coder.nullcopy(eml_expand(bzero, size(a)));
if isscaleddouble(a)
    b(:) = a*a;
    for idx = 3:1:k
        b(:) = b*a;
    end
    y = eml_cast(b, ty);
else
    as = stripscaling(a);
    b(:) = as*as;
    for idx = 3:1:k
        b(:) = b*as;
    end
    y = reinterpretcast(b, ty);
end


function y = compute_mpower_typeb(a, k, n)%#codegen
coder.internal.prefer_const(k);
[fmA, ty, ega] = init_mpower(a, k);
ty1 = numerictype(ty, 'FractionLength', 0);
bzero = eml_cast(ega, ty1, fmA);
b = coder.nullcopy(eml_expand(bzero, size(a)));
if isscaleddouble(a)
    b(:) = a*a;
    for idx = 3:1:k
        b(:) = b*a;
    end
    y = eml_cast(b, ty);
else    
    as = stripscaling(a);
    b(:) = stripscaling(as*as);
    ta = numerictype(a);
    ta1 = numerictype(ta,'FractionLength',ta.WordLength+n);
    as1 = reinterpretcast(as,ta1);
    for idx = 3:1:k
        b(:) = b*as1;
    end %for
    y = reinterpretcast(b, ty);
end

function y = compute_mpower_typec(a, k)%#codegen
% either 'ProductMode' or 'SumMode' is 'SpecifyPrecision'
coder.internal.prefer_const(k);
if eml_is_const(k)
    b1 = a*a;
    if (k == 2)
        y = b1;
    else
        [fmA, ty, ega] = init_mpower(a, k);
        yzero = eml_cast(ega, ty, fmA);
        y = coder.nullcopy(eml_expand(yzero, size(a)));
        y(:) = b1;
        
        for pwridx = 3:1:k
            y(:) = y*a;
        end
    end
else
    % the SumMode' must be 'SpecifyPrecision'
    [fmA, ty, ega] = init_mpower(a, k);
    yzero = eml_cast(ega, ty, fmA);
    y = coder.nullcopy(eml_expand(yzero, size(a)));
    
    if k == 0
        y(:) = fi(eye(size(a,1)));
    elseif k == 1
        y(:) = a;
    else % k >= 2
        y(:) = a*a;
    end
    for pwridx = 3:1:k
        y(:) = y*a;
    end
end


function [fmA, ty, ega] = init_mpower(a, k)
coder.inline('always');
coder.extrinsic('emlGetNTypeForMpower');
coder.internal.prefer_const(k);
maxWL = getFixedPointWidthLimit();
fmA = eml_fimath(a);
if eml_is_const(size(a))

    szAc = coder.const(size(a));
else

    szAc = coder.const([2 2]);
end
ega = coder.internal.scalarEg(a);
if eml_is_const(k)
    ty = coder.const(emlGetNTypeForMpower(ega,szAc,k,fmA,maxWL));
else
    if isscaledtype(a) 
        % fixed-point or scaled double data type
        isSumModeSpecPrec = coder.const(strcmpi(a.SumMode,'SpecifyPrecision'));
        eml_invariant(isSumModeSpecPrec,...
                      'fixed:coder:mpowerSumModeForNonConstExponent');
        ty = numerictype(numerictype(a), 'WordLength', a.SumWordLength,...
                         'FractionLength',a.SumFractionLength);
    else
        ty = numerictype(a);
    end
end
