function y = nchoosek(x,k)
%MATLAB Code Generation Library Function

%   Limitations:
%   1. Known difference from MATLAB when x is a scalar and not a
%   non-negative integer. MATLAB interprets x as a set in that case whereas
%   MATLAB Coder interprets x is an invalid scalar input for a binomial
%   coefficient.

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

coder.internal.assert(nargin == 2,'MATLAB:minrhs');

coder.internal.assert(coder.internal.isConstTrue(isscalar(k)), ...
    'Coder:toolbox:InputMustBeFixedSizeScalar','k');
coder.internal.assert(coder.internal.isBuiltInNumeric(k) && isscalar(k) && isreal(k) && ...
    k >= 0 && k == floor(k), ...
    'MATLAB:nchoosek:InvalidArg2');
coder.internal.assert(isvector(x), ...
    'MATLAB:nchoosek:InvalidArg1');
n = numel(x);
if isscalar(x) && isnumeric(x)
   
   if ~coder.target('MATLAB') && coder.internal.isConst(x) && ...
        coder.internal.isConst(k) && coder.internal.isCompiled
       y = coder.const(@feval,'nchoosek',x,k);
       return;
   end
   if isinteger(x)
       coder.internal.assert(isa(x,class(k)) || isa(k,'double'),...
           'MATLAB:nchoosek:mixedIntegerTypes');
   elseif isinteger(k)
       coder.internal.assert(isa(x,'double'),...
           'MATLAB:nchoosek:mixedIntegerTypes');
   end
   eg = coder.internal.scalarEg(x,k);
   coder.internal.assert(k <= x,'MATLAB:nchoosek:KOutOfRange');
   isInt64 = isa(eg,'int64') || isa(eg,'uint64');
   m = cast(x(1),'like',eg);
   r = cast(k(1),'like',eg);
   if r > m - r
       r = m - r;
   end
   coder.internal.errorIf(~isInt64 && m > flintmax,...
       'MATLAB:nchoosek:NOutOfRange');
   if r == 0
       y = cast(1,'like',eg);
   elseif r == 1
       y = m;
   elseif isInt64
       y = nCkInt64(m,r);
   elseif isinteger(eg)
       y = cast(nCk(double(m),double(r),true),'like',eg);
   else
       y = nCk(m,r,false);
   end
   
elseif k > numel(x)
    y = eml_expand(coder.internal.scalarEg(x),[0,k]);
else
    ONE = ones(coder.internal.indexIntClass);
    nint = eml_cast(n,coder.internal.indexIntClass,'to zero','wrap');
    kint = eml_cast(k,coder.internal.indexIntClass,'to zero','wrap');
    nrows = eml_cast(nchoosek(n(1),double(k)),coder.internal.indexIntClass,'to zero','wrap');
    y = coder.nullcopy(eml_expand(coder.internal.scalarEg(x),[nrows,kint]));
    
    comb = ONE:kint;
    icomb = kint; % index into comb
    nmkpi = nint; % n-k+i
    for row = ONE:nrows
        for col = ONE:kint
            y(row,col) = x(comb(col));
        end
        % Compute the next combination of n things taken k at a time.
        % The "first" combination should be 1:k.  The "last" combination
        % will be (n-k+1):n.
        if icomb > 0
            combj = coder.internal.indexPlus(comb(icomb),ONE);
            comb(icomb) = combj;
            if combj < nmkpi
                for j = coder.internal.indexPlus(icomb,ONE):kint
                    combj = coder.internal.indexPlus(combj,ONE);
                    comb(j) = combj;
                end
                icomb = kint;
                nmkpi = nint;
            else
                icomb = coder.internal.indexMinus(icomb,ONE);
                nmkpi = coder.internal.indexMinus(nmkpi,ONE);
            end
        end
    end
end

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

function y = nCk(n,k,inttype)
% n!/(k!*(n-k)!), 0 <= k <= n.
% Warn about inexact result if the computed value exceeds warnabove.
y = coder.nullcopy(coder.internal.scalarEg(n,k));
if ~isfinite(n) || ~isfinite(k)
    y(1) = coder.internal.nan;
elseif k > 1000
    y(1) = coder.internal.inf;
else
    if isa(y,'single')
        y = n;
        for j = 2:k
            y = y * ((n-j+1)/j);
        end
        y = round(y);
        if ~(y <= flintmax(class(y)))
            if ~coder.setwarningflag
                coder.internal.warning('Coder:MATLAB:nchoosek_LargeCoefficient');
            end
        end
    else
        maxRelErr = 0;
        y = n;
        for j = 2:k
            y = y*((n-j+1)/j);
            if ~(y < 1/(4*eps))
                maxRelErr = maxRelErr + 2*eps;
            end
            y = round(y);
        end
        if maxRelErr ~= 0 && y <= 2*double(intmax('uint64')) && ~coder.internal.isAmbiguousTypes
            yint = nCkInt64(uint64(n),uint64(k));
            if yint <= intmax('uint64')
                y = double(yint);
                maxRelErr = 0;
            end
        end
        if ~inttype
            if y > flintmax
                maxRelErr = max(maxRelErr,eps);
            end

            if maxRelErr ~= 0 && isfinite(y) && ~coder.setwarningflag
                coder.internal.warning('MATLAB:nchoosek:LargeCoefficient',...
                    coder.internal.flt2str(maxRelErr), coder.internal.flt2str(ceil(maxRelErr*y)))
            end
        end
    end
end

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

function y = nCkInt64(n,k)
imax = intmax(class(n));
y = ones(class(n));
for i = 1:k
    q = eml_rdivide(y,i,class(n),'spill','to zero');
    if q >= eml_rdivide(imax,n,class(n),'spill','to zero')
        y = imax;
        return
    end
    r = mod(y,i);
    y = eml_plus(eml_times(q,n,class(n),'spill','to zero'),...
        eml_rdivide(eml_times(r,n,class(n),'spill','to zero'),i,....
        'spill','to zero'),'spill',class(n), 'to zero');
    n = eml_minus(n,1,class(n),'spill','to zero');
end
