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.
%   
%   2. The argument k cannot be a 64-bit integer.  If x is a scalar then it
%   is not allowed to be a 64-bit integer.

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

eml_invariant(nargin == 2, ...
    'MATLAB:minrhs');

eml_invariant(eml_is_const(isscalar(k)), ...
    'Coder:toolbox:InputMustBeFixedSizeScalar','k');
kIsInt64 = coder.internal.is64BitInteger(k);
eml_invariant(~kIsInt64, ...
    'Coder:toolbox:Int64NotSupported', ...
    'IfNotConst', 'Fail');
eml_invariant(coder.internal.isBuiltInNumeric(k) && isscalar(k) && isreal(k) && ...
    k >= 0 && k == floor(k), ...
    'MATLAB:nchoosek:InvalidArg2');
eml_invariant(isvector(x), ...
    'MATLAB:nchoosek:InvalidArg1');
n = numel(x);
if isscalar(x) && isnumeric(x)
    xIsInt64 = coder.internal.is64BitInteger(x);
    eml_invariant(~xIsInt64, ...
                  'Coder:toolbox:Int64NotSupported', ...
                  'IfNotConst', 'Fail');

    hasInt64 = kIsInt64 || xIsInt64;
    if k > x(1)
        coder.internal.error('MATLAB:nchoosek:KOutOfRange');
        y = coder.internal.scalarEg(x,k);
        y(1) = coder.internal.nan;
    else
        eml_invariant(hasInt64 || x(1) <= flintmax, ...
                              'MATLAB:nchoosek:NOutOfRange');
        y = nCk(x(1),k);
    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(nCk(n,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)
% 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 isinteger(n) || isinteger(k)
    y(1) = nCk(double(n),double(k));
elseif ~isa(k,class(n))
    y(1) = nCk(single(n),single(k));
else
    if isfinite(n) && isfinite(k)
        if k > eml_rdivide(n,2)
            k = n - k;
        end
        if k > 1000
            y(1) = coder.internal.inf;
        else
            y = ones(class(y));
            nmk = n - k;
            for j = 1:k
                y = y * eml_rdivide(j+nmk,j);
            end
            y = round(y);
        end
    else
        y(1) = coder.internal.nan;
    end
    % Note that MATLAB has an algorithm to produce correct results between
    % flintmax and intmax for 64-bit integers. Consider adding that when we
    % support int64 and uint64.
    if ~(y <= flintmax(class(y)))
        if ~coder.setwarningflag
            coder.internal.warning('Coder:MATLAB:nchoosek_LargeCoefficient');
        end
    end
end

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