function yreturn = sqrt(x,var1,var2)
% SQRT Fixed-point square root function for MATLAB Code Generation

% Copyright 2006-2020 The MathWorks, Inc.
%#codegen

coder.extrinsic('eml_fisqrt_helper');
    
% This function accepts mxArray input argument
eml_allow_mx_inputs;      

if eml_ambiguous_types
    yreturn = coder.ignoreConst(zeros(size(x)));
    return;
end

% Error if incorrect number of inputs
eml_invariant((nargin >= 1 && nargin <= 3), ...
    'fixed:coder:incorrectNumberOfInputs');

% Error if complex, or slope bias or negative
eml_invariant(isreal(x), ...
    'fixed:fi:unsupportedComplexInput', 'SQRT');

if isscaledtype(x)
    % FI objects with Fixed datatype

    % Parse the inputs
    switch nargin
      case 1 % sqrt(x) get Ty and Fy and method = bisection

        if eml_is_const(x)
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,x);
        else % ~eml_is_const(x)
            data = localGetConstX(x);
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,data);
        end
        eml_assert(isempty(errMsg),errMsg); % OK for i18n (translated errMsg)
        
      case 2 % sqrt(x,T), sqrt(x,F) or sqrt(x,method)
        eml_invariant(eml_is_const(var1), ...
            'fixed:coder:inputMustBeConstantOrAFI', ...
            'var1', 'sqrt(x,var1)');
        if eml_is_const(x)
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,x,var1);
        else % ~eml_is_const(x)
            data = localGetConstX(x);
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,data,var1);
        end
        eml_assert(isempty(errMsg),errMsg); % OK for i18n (translated errMsg)
        
      case 3 % sqrt(x,T,F), sqrt(x,T,method) or sqrt(x,F,method)
        eml_invariant(eml_is_const(var1) && eml_is_const(var2),...
            'fixed:coder:inputMustBeConstantOrAFI', ...
            'var1,var2', 'sqrt(x,var1,var2)');
        if eml_is_const(x)
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,x,var1,var2);
        else % ~eml_is_const(x)
            data = localGetConstX(x);
            [Ty,Fy,fimathSpecified,methodName,errMsg] = coder.const(@eml_fisqrt_helper,data,var1,var2);
        end
        eml_assert(isempty(errMsg),errMsg); % OK for i18n (translated errMsg)
        
    end

    % Error if x or Ty is slope-bias scaled
    if (isslopebiasscaled(numerictype(x)) || isslopebiasscaled(Ty))
        eml_invariant(false, ...
            'fixed:coder:sqrtSlopeBiasInputsNotSupported');
    end

    Tx = numerictype(x);
    if isscaleddouble(Tx) || isscaleddouble(Ty)
        Ty = numerictype(Ty,'DataType','scaleddouble');
        y = fi(sqrt(double(x)),Ty,Fy);
    else
        % Initialize output y with Ty and Fy to prevent compiler error
        y = fi(zeros(size(x)),Ty,Fy);
        % Compute sqrt based on the square root method
        switch methodName
            case 0 % bisection
                y = bisectionSqrt(x,y,Ty,Fy);
        end
    end

elseif isfloat(x) 
    % True Double or True Single FI
    fimathSpecified = coder.const(false);
    % Parse the inputs
    switch nargin
      case 1 % sqrt(x)
        Ty = coder.const(eml_typeof(x));
        Fy = coder.const(eml_fimath(x));
      case 2 % sqrt(x,T), sqrt(x,F) or sqrt(x,method)->This syntax is not yet
             % supported, hence error out
        eml_invariant(eml_is_const(var1),'fixed:coder:sqrt_X_Var1_Var1MustBeConstant');
        if isnumerictype(var1)
            Ty = localGetSqrtNumericType(x,var1);
            Fy = coder.const(eml_fimath(x));
        elseif isfimath(var1)
            Ty = coder.const(eml_typeof(x));
            Fy = coder.const(var1);
        else
            eml_invariant(false, ...
                'fixed:coder:fiInvalidMethodSignature','SQRT');
        end
        
      case 3 % sqrt(x,T,F), sqrt(x,T,method) or sqrt(x,F,method)
        eml_invariant(eml_is_const(var1) && eml_is_const(var2),...
            'fixed:coder:inputMustBeConstantOrAFI', ...
            'var1,var2', 'sqrt(x,var1,var2)');
        if ischar(var2) % sqrt(x,T,method) or sqrt(x,F,method)
            if isnumerictype(var1)
                Ty = localGetSqrtNumericType(x,var1);
                Fy = coder.const(eml_fimath(x));
            elseif isfimath(var1)
                Ty = coder.const(eml_typeof(x));
                Fy = coder.const(var1);
            else
                eml_invariant(false, ...
                    'fixed:coder:fiInvalidMethodSignature','SQRT');
            end
        else % sqrt(x,T,F)
            Ty = localGetSqrtNumericType(x,var1);
            Fy = coder.const(var2);
        end
    end

    % Call MATLAB realsqrt as complex inputs are not supported.
    dType  = eml_fi_getDType(x);
    xDType = eml_cast(x,dType);
    yDType = realsqrt(xDType);
    y      = eml_cast(yDType,Ty,Fy);    
    
else
    % FI datatype not supported
    eml_fi_assert_dataTypeNotSupported('SQRT');
end

% If the input fi does not have an attached fimath, then the output doesn't either
% Also if a fimath is specified as one of the input parameters, the output is then fimathless.
if coder.const(~eml_fimathislocal(x)) || fimathSpecified
    yreturn = eml_fimathislocal(y,false);
else
    yreturn = y;
end

%---------------------------------------------------------------------------------
function y = bisectionSqrt(x,yinit,Ty,Fy)

coder.extrinsic('emlGetNTypeForTimes');
coder.extrinsic('eml_iscomplexroundmode_helper');

% Binary search method of finding the square root

% This function accepts mxArray input argument
eml_allow_mx_inputs;      

% Get numel of x
Nx = numel(x);

% Initialize y 
y = yinit;
maxWL = getFixedPointWidthLimit();

% Initialize other temporary variables
Tsquare  = coder.const(emlGetNTypeForTimes(Ty,Ty,Fy,true,true,maxWL));
wlTsquare = coder.const(get(Tsquare,'WordLength'));
flTsquare = coder.const(get(Tsquare,'FractionLength'));
yWL = coder.const(get(Ty,'WordLength'));
yFL = coder.const(get(Ty,'FractionLength'));

Ftemp = fimath(Fy,'ProductMode','SpecifyPrecision',...
               'ProductWordLength',wlTsquare,...
               'ProductFractionLength',flTsquare,...
               'SumMode','SpecifyPrecision',...
               'SumWordLength',wlTsquare,...
               'SumFractionLength',flTsquare,...
               'CastBeforeSum',true);
FtempForPlusLSB = fimath(Fy,'ProductMode','SpecifyPrecision',...
                         'ProductWordLength',wlTsquare,...
                         'ProductFractionLength',flTsquare,...
                         'SumMode','SpecifyPrecision',...
                         'SumWordLength',yWL,...
                         'SumFractionLength',yFL,...
                         'CastBeforeSum',true);
ytemp = fi(0,Ty,FtempForPlusLSB);
lsbY = lsb(ytemp);
ytempSquare = fi(0,Tsquare,Ftemp);
ytempPlusLsb = fi(0,Ty,Ftemp);
ytempPlusLsbSquare = fi(0,Tsquare,Ftemp);
tempx = fi(0,numerictype(x),Ftemp);

% Determine bit location 
wlTy = coder.const(int32(get(Ty,'WordLength')));
if issigned(y)
    wn = wlTy - 1;
else    
    wn = wlTy;
end

% Compute output y
for k = 1:Nx
    % If x is negative or 0, then y is 0, otherwise y is computed.
    % In MATLAB and MEX, we error out for x < 0.
    % In generated C code we return 0 and do not error out because
    % a C run time error is not useful.
    coder.internal.assert(x(k)>=0,'fixed:fi:sqrtNoNegative');
    if (x(k) <= 0)
        y(k) = 0;
    else

        % Find y(k)
        for i=wn:-1:1
            ytemp(1) = bitset(y(k),i);
            ytempSquare(1) = ytemp.*ytemp;
            if (ytempSquare <= x(k))
                y(k) = ytemp;
            end
        end
        % "Polish" the final result only when roundmode
        % is ceiling, round, nearest or convergent
        if (coder.const(eml_iscomplexroundmode_helper(Fy)))

            if (y(k) < realmax(y))

                ytemp(1) = y(k);
                ytempSquare(1) = ytemp.*ytemp;
                ytempPlusLsb(1) = ytemp+lsbY;
                tempx(1) = x(k);
                if coder.const(eml_iscomplexroundmode_helper(Fy,'ceil'))

                    if (ytempSquare < tempx)
                        y(k) = ytempPlusLsb;
                    end
                else % for round, nearest or convergent
                    
                    ytempPlusLsbSquare(1) = ytempPlusLsb*ytempPlusLsb;

                    if ( ytempPlusLsbSquare-tempx < tempx-ytempSquare)
                        y(k) = ytempPlusLsb;
                    end
                
                end
            end

        end
    end
end

%--------------------------------------------------------------------------------
function data = localGetConstX(x)
% Local function to return a const fi with same type and fimath as x
data = eml_cast(0,numerictype(x),fimath(x));

%--------------------------------------------------------------------------------
function Ty = localGetSqrtNumericType(x,T)
% Returns numerictype of output of sqrt function. This local function should
% be used only when input x is a float fi. The output, Ty is const.
eml_transient;
if ( coder.const(isnumerictype(T)) && coder.const(isfloat(T)) ) 
    % NumericType is float (double/single)
    Ty = coder.const(T);
else
    Ty = coder.const(eml_typeof(x));
end

%--------------------------------------------------------------------------------
