function y = colon(a,d,b)
%MATLAB Code Generation Library Function

%   NOTES:
%   1. Does not accept complex inputs.
%   2. Does not accept logical D.
%   3. Does not accept vector inputs.
%   4. Inputs must be constants.
%   5. Single precision results are filled using single-precision
%   arithmetic.
%   6. Half precision results are filled using double precision arithmetic.

%   Copyright 1984-2023 The MathWorks, Inc.
%#codegen

coder.internal.allowHalfInputs;
coder.internal.userReadableName('');
coder.inline('always');
eml_invariant(nargin >= 2,'MATLAB:minrhs');
coder.internal.prefer_const(a);
coder.internal.prefer_const(d);
eml_invariant(coder.internal.isBuiltInNumeric(a) || ischar(a) || isa(a,'half'), ...
    'Coder:toolbox:unsupportedClass','colon',class(a));
if nargin < 3
    if ischar(a)
        y = colon(a,char(1),d);
    else
        y = colon(a,ones('like',a),d);
    end
    return
end
coder.internal.prefer_const(b);
if ~( ...
        eml_is_const(isscalar(a)) && isscalar(a) && ...
        eml_is_const(isscalar(d)) && isscalar(d) && ...
        eml_is_const(isscalar(b)) && isscalar(b))
    eml_invariant(isscalar(a) && isscalar(d) && isscalar(b), ...
        'Coder:toolbox:colon_operandsMustBeScalars');
    y = colon(a(1),d(1),b(1));
    return
end
eml_invariant(coder.internal.isBuiltInNumeric(d) || ischar(d) || isa(d,'half'), ...
    'Coder:toolbox:unsupportedClass','colon',class(d));
eml_invariant(coder.internal.isBuiltInNumeric(b) || ischar(b) || isa(b,'half'), ...
    'Coder:toolbox:unsupportedClass','colon',class(b));
eml_invariant(coder.const(~coder.internal.eml_option_eq('VariableSizing','DisableInInference')) || ...
    (eml_is_const(a) && eml_is_const(d) && eml_is_const(b)), ...
    'Coder:toolbox:colon_6', ...
    'IfNotConst','Fail');
eml_invariant(isreal(a) && isreal(d) && isreal(b), ...
    'Coder:toolbox:colon_7', ...
    'IfNotConst','Fail');
isda = isa(a,'double');
isdb = isa(b,'double');
isdd = isa(d,'double');
eml_invariant( ...
    (isda || isdd || isa(a,class(d))) && ...
    (isda || isdb || isa(a,class(b))) && ...
    (isdd || isdb || isa(d,class(b))), ...
    'MATLAB:colon:mixedNonDoubleOperands', ...
    'IfNotConst','Fail');
yzero = coder.const(zerosum(a,d,b));
if ischar(a) || ischar(d) || ischar(b)
    eml_invariant(ischar(a) && ischar(b), ...
        'Coder:toolbox:colon_9', ...
        'IfNotConst','Fail');
    eml_invariant(ischar(d) || isa(d,'double'), ...
        'Coder:toolbox:colon_10', ...
        'IfNotConst','Fail');
    eml_invariant( ...
        (~isda || floor(a) == a) && ...
        (~isdd || floor(d) == d) && ...
        (~isdb || floor(b) == b), ...
        'Coder:MATLAB:colon_scalarDoubleMustBeIntegerValued');
    y = char(eml_integer_colon_dispatcher(int8(a),double(d),int8(b)));
elseif isinteger(a) || isinteger(d) || isinteger(b)
    eml_invariant( ...
        (~isda || floor(a) == a) && ...
        (~isdd || floor(d) == d) && ...
        (~isdb || floor(b) == b), ...
        'Coder:MATLAB:colon_scalarDoubleMustBeIntegerValued');
    y = eml_integer_colon_dispatcher(a,d,b);
elseif isnan(a) || isnan(d) || isnan(b)
    y = coder.internal.nan('like',yzero);
elseif (d == 0) || (a < b && d < 0) || (b < a && d > 0)
    y = zeros(1,0,'like',yzero);
elseif (isinf(a) || isinf(b)) && (isinf(d) || ( ...
        ~coder.internal.isConstFalse(isinf(a)) && ...
        ~coder.internal.isConstFalse(isinf(b)) && ...
        a == b))
    % The uses of coder.internal.isConstFalse above help the
    % constant-folder in cases like 1:b, where b is variable and possibly
    % infinite. Otherwise, with current limitations of range analysis, we
    % would generate the equivalent of
    %     elseif isinf(b) && 1 == b
    y = coder.internal.nan('like',yzero);
elseif isinf(d)
    y = cast(a,'like',yzero);
elseif isa(a,'half') || isa(b,'half') || isa(d,'half')
    y = eml_half_colon(a,d,b);
elseif floor(a) == a && floor(d) == d
    cls = class(yzero);
    if ~(isa(a,cls) && isa(b,cls) && isa(d,cls))
        checkrange(a,d,b);
    end
    icls = coder.internal.indexIntClass;
    TWO = coder.internal.indexInt(2);
    IFMAX = coder.const(irdivide(intmax(icls),TWO));
    FMAX = coder.const(cast(IFMAX,cls));
    SIGNED = coder.const(intmin(icls) < 0);
    if isda || isdb || isdd || ~SIGNED || ...
            abs(a) >= FMAX || abs(d) >= FMAX || abs(b) >= FMAX
        y = cast(eml_colonobj(double(a),double(d),double(b)).',cls);
    else
        ia = icast(a,icls);
        ib = icast(fix(b),icls);
        id = icast(d,icls);
        n = iplus(ones(icls),irdivide(iminus(ib,ia),id));
        y = coder.nullcopy(zeros(1,n,cls));
        for k = 1:n
            y(k) = ia + (k - 1)*id;
        end
    end
else
    y = eml_float_colon(a,d,b);
end

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

function y = eml_integer_colon_dispatcher(a,d,b)
% Compute Y = A:D:B for the integer classes.
coder.internal.prefer_const(a,d,b);
if eml_ambiguous_types
    % This will need fixing when we support 64-bit integers.
    y = colon(double(a),double(d),double(b));
    return
end
ZERO = zerosum(a,d,b);
cls = class(ZERO);
checkrange(a,d,b);
ai = icast(a,cls);
bi = icast(b,cls);
if d < intmin(cls) || d > intmax(cls)
    y = eml_integer_colon_d_out_of_range(ai,d,bi);
elseif coder.internal.isaUint(ZERO)
    y = eml_unsigned_integer_colon(ai,icast(d,cls),bi);
else
    y = eml_signed_integer_colon(ai,icast(d,cls),bi);
end

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

function y = eml_unsigned_integer_colon(a,d,b)
% Integer colon.  This function assumes that inputs a, d, and b
% belong to the same integer class.
coder.internal.prefer_const(a,d,b);
cls = class(a);
n = integer_colon_length(a,d,b);
y = coder.nullcopy(zeros(1,n,cls));
for k = 1:n
    y(k) = iplus(a,itimes(icast(coder.internal.indexMinus(k,1),cls),d));
end

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

function y = eml_signed_integer_colon(a,d,b)
% Integer colon.  This function assumes that inputs a, d, and b
% belong to the same integer class.
coder.internal.prefer_const(a,d,b);
cls = class(a);
n = integer_colon_length(a,d,b);
y = coder.nullcopy(zeros(1,n,cls));
if n > 0
    y(1) = a;
    yk = a;
    for k = 2:n
        yk = iplus(yk,d);
        y(k) = yk;
    end
end

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

function y = eml_integer_colon_d_out_of_range(a,d,b)
% Edge case colon expressions where a and b belong to an integer class but
% d is a 'double' that is not in the range of the integer class.
% Consequently, the code below assumes d < intmin(cls) || d > intmax(cls).
coder.internal.prefer_const(a,d,b);
cls = class(a);
if d > intmax(cls)
    % d > intmax for both signed and unsigned integer classes
    if a > b
        y = zeros(1,0,cls);
    elseif a < 0
        y2d = double(a) + d;
        if y2d <= b
            y = [a,icast(y2d,cls)];
        else
            y = a;
        end
    else
        y = a;
    end
elseif coder.internal.isaUint(a)
    % d < intmin for unsigned classes
    if a < b
        y = zeros(1,0,cls);
    elseif -d <= intmax(cls)
        % Step value is in range.
        absd = icast(-d,cls);
        % length(a:d:b) == length(b:-d:a)
        n = integer_colon_length(b,absd,a);
        y = coder.nullcopy(zeros(1,n,cls));
        for k = 1:n
            y(k) = iminus(a,itimes(icast(coder.internal.indexMinus(k,1),cls),absd));
        end
    else
        % Step value is out of range.
        y = a;
    end
else
    % d < intmin for signed classes
    if a < b
        y = zeros(1,0,cls);
    elseif a > 0
        y2d = double(a) + d;
        if y2d >= b
            y = [a,icast(y2d,cls)];
        else
            y = a;
        end
    else
        y = a;
    end
end

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

function n = integer_colon_length(a,d,b)
% Length of (a:d:b).  All inputs must belong to the same integer class.
% The output is coder.internal.indexIntClass.
coder.internal.prefer_const(a,d,b);
coder.inline('always');
if d < 0
    % Change direction. Note that length(a:d:b) == length(b:-d:a).
    % We must handle the edge case d = intmin(cls) because -d wraps.
    if d == intmin(class(a))
        if a < b
            % Return empty in the degenerate case
            n = zeros(coder.internal.indexIntClass);
        elseif a >= 0 && iplus(a,d) >= b
            n = cast(2,coder.internal.indexIntClass);
        else
            n = ones(coder.internal.indexIntClass);
        end
    else
        negd = ineg(d);
        n = integer_colon_length_nonnegd(b,negd,a);
    end
else
    n = integer_colon_length_nonnegd(a,d,b);
end

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

function n = integer_colon_length_nonnegd(a,d,b)
% Length of (a:d:b).  All inputs must belong to the same integer class.
% This function returns zero if d <= 0 || a > b.
coder.internal.prefer_const(a,d,b);
coder.inline('always');
cls = class(a);
ucls = coder.internal.unsignedClass(cls);
IDXMAX = intmax(coder.internal.indexIntClass);
IDX_CLASS_COVERS = intmax(ucls) < IDXMAX;
SIGNED = ~coder.internal.isaUint(a);
if d <= 0 || b < a
    n = zeros(coder.internal.indexIntClass);
else
    if SIGNED && a < 0 && b >= 0
        % Number of points may be greater than intmax(cls).
        a = iplus(a,ones(cls));
        span = iplus(icast(b,ucls),ones(ucls));
        span = iplus(span,icast(ineg(a),ucls));
    else
        % Although a <= b is guaranteed here, the compiler is not able to
        % determine that b - a is always non-negative when these inputs are
        % variable. Considering that ICAST avoids saturation logic with the
        % 'spill' option, we need to indicate to the compiler that b - a is
        % non-negative in order to get good bounds.
        bma = iminus(b,a);
        assert(bma >= 0); %<HINT>
        span = icast(bma,ucls);
    end
    nm1 = irdivide(span,icast(d,ucls));
    assert_pmaxsize(IDX_CLASS_COVERS || nm1 < IDXMAX);
    n = coder.internal.indexPlus(icast(nm1,coder.internal.indexIntClass),1);
end

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

function y = eml_float_colon(a,d,b)
% Compute Y = A:D:B for a possibly non-integer stride D.
% Size-determining computations are performed in double precision floating
% point, but unlike MATLAB, we use single precision for the arithmetic in
% the fill operation if A or D is single precision. For half precision we
% use double precision arithmetic to fill in the values in line with the
% MATLAB implementation.
coder.internal.prefer_const(a,d,b);
ytype = zerosum(a,d,b);
if isa(ytype,'half')
    d1 = double(d);
    EPS = coder.const(feval('eps',ones('like',ytype)));
    [n,a1,b1,n_too_large] = float_colon_length(double(a),d1,double(b),EPS);
else
    d1 = d;
    EPS = coder.const(feval('eps',ones('like',ytype)));
    [n,a1,b1,n_too_large] = float_colon_length(a,d1,b,EPS);
end
assert_pmaxsize(~n_too_large);
y = coder.nullcopy(zeros(1,n,'like',ytype));
if n > 0
    y(1) = a1;
    if n > 1
        y(n) = b1;
        nm1 = n - 1;
        nm1d2 = coder.internal.indexDivide(nm1,2);
        coder.gpu.kernel;
        for k = 1:(nm1d2 - 1)
            kd = cast(k,'like',d1)*d1;
            y(k + 1) = a1 + kd;
            y(n - k) = b1 - kd;
        end
        if 2*nm1d2 == nm1 % N is odd.
            y(nm1d2 + 1) = eml_rdivide(a1 + b1,2);
        else % N is even.
            kd = cast(nm1d2,'like',d1)*d1;
            y(nm1d2 + 1) = a1 + kd;
            y(nm1d2 + 2) = b1 - kd;
        end
    end
end

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

function [n,anew,bnew,n_too_large] = float_colon_length(a,d,b,EPS)
% Calculate the length of the colon expression and adjust B, if necessary,
% so that BNEW = A + (N-1)*D to working precision.  The start index, A, is
% not adjusted unless one or more of the inputs are non-finite, in which
% case ANEW = NaN.
coder.internal.prefer_const(a,d,b);
% This function assumes the following cases have already been handled:
% if isnan(a) || isnan(d) || isnan(b)
% elseif (d == 0) || (a < b && d < 0) || (b < a && d > 0)
% elseif isinf(a) || isinf(b)
% elseif isinf(d)
% end
anew = a;
ndbl = floor(eml_rdivide(double(b)-double(a),double(d)) + 0.5);
apnd = double(a) + ndbl*double(d);
if d > 0
    cdiff = apnd - double(b);
else
    cdiff = double(b) - apnd;
end
thresh = 2*double(EPS)*maxabs(double(a),double(b));
if abs(cdiff) < thresh
    ndbl = ndbl + 1;
    bnew = b;
elseif cdiff > 0 % overshot
    bnew = cast(double(a)+(ndbl-1)*double(d),'like',b);
else % undershot
    ndbl = ndbl + 1;
    bnew = cast(apnd,'like',b);
end
n_too_large = ndbl > intmax(coder.internal.indexIntClass);
% Although ndbl is always non-negative, range analysis needs a hint.
if ndbl >= 0
    n = eml_cast(ndbl,coder.internal.indexIntClass,'to zero','spill');
else
    n = zeros(coder.internal.indexIntClass);
end

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

function checkrange(a,d,b)
% Assert that A and B are in range for class(A:D:B).
coder.internal.prefer_const(a,d,b);
yzero = zerosum(a,d,b);
outcls = class(yzero);
if isa(yzero,'float')
    ubnd = realmax(outcls);
    lbnd = -ubnd;
else
    lbnd = intmin(outcls);
    ubnd = intmax(outcls);
end
eml_invariant( ...
    (isa(a,outcls) || (a <= ubnd && a >= lbnd)) && ...
    (isa(b,outcls) || (b <= ubnd && b >= lbnd)), ...
    'Coder:MATLAB:colon_OutOfRange');

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

function c = maxabs(a,b)
% c = max(abs(a),abs(b)) without going through the machinery of MAX.
coder.internal.prefer_const(a,b);
c = eml_max(abs(a),abs(b));

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

function zero = zerosum(a,d,b)
% Returns the scalar 0 of the combined class of the input arguments.
zero = cast(0,'like',a) + cast(0,'like',d) + cast(0,'like',b);

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

function assert_pmaxsize(p)
coder.internal.prefer_const(p);
eml_invariant(p, ...
    'Coder:MATLAB:pmaxsize');

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

function b = icast(a,cls)
% b = cast(a,cls) without saturation.
coder.inline('always');
coder.internal.prefer_const(a,cls);
b = eml_cast(a,cls,'to zero','spill');

function b = ineg(a)
% b = -a without saturation.
coder.inline('always');
coder.internal.prefer_const(a);
b = eml_uminus(a,class(a),'spill');

function c = iplus(a,b)
% c = a + b in class(a) integer arithmetic without checking for saturation.
coder.inline('always');
coder.internal.prefer_const(a,b);
c = eml_plus(a,b,class(a),'spill');

function c = iminus(a,b)
% c = a - b in class(a) integer arithmetic without checking for saturation.
coder.inline('always');
coder.internal.prefer_const(a,b);
c = eml_minus(a,b,class(a),'spill');

function c = itimes(a,b)
% c = a.*b in class(a) integer arithmetic without checking for saturation.
coder.inline('always');
coder.internal.prefer_const(a,b);
c = eml_times(a,b,class(a),'spill');

function c = irdivide(a,b)
% c = a./b in class(a) integer arithmetic.  Rounding is to zero and there
% is no saturation checking.
coder.inline('always');
coder.internal.prefer_const(a,b);
c = eml_rdivide(a,b,class(a),'to zero','spill');

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

function y = eml_half_colon(a,d,b)
% For half precision we use double precision arithmetic to fill in the
% values in line with the MATLAB implementation.
coder.internal.prefer_const(a,d,b);
if floor(a) == a && floor(b) == b && floor(d) == d
    y = eml_half_colonobj(a,d,b);
else
    y = eml_float_colon(a,d,b);
end

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

function y = eml_half_colonobj(a,d,b)
coder.internal.prefer_const(a,d,b);
a1 = double(a);
d1 = double(d);
b1 = double(b);
n = coder.internal.indexInt(length(eml_colonobj(a1,d1,b1)));
y = coder.nullcopy(zeros(1,n,'half'));
for k = 1:n
    y(k) = a1 + double(k - 1)*d1;
end

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