function c = bsxfun(fun,a,b)
%MATLAB Code Generation Library Function

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

coder.internal.allowHalfInputs;
eml_invariant(nargin == 3,'MATLAB:minrhs');
eml_invariant(isa(fun,'function_handle'), ...
    'Coder:toolbox:bsxfun_2');
ONE = coder.internal.indexInt(1);
nda = coder.internal.indexInt(eml_ndims(a));
ndb = coder.internal.indexInt(eml_ndims(b));
ND = coder.const(eml_max(nda,ndb));
% Compute the output size vector and check for consistency.
csz = zeros(1,ND,coder.internal.indexIntClass);
iscompatible = true;
for k = coder.unroll(ONE:ND)
    sak = coder.internal.indexInt(size(a,k));
    sbk = coder.internal.indexInt(size(b,k));
    if coder.internal.isConst(sak)
        if coder.const(sak == ONE)
            sck = sbk;
        else
            sck = sak;
            iscompatible = iscompatible && (sbk == ONE || sak == sbk);
        end
    elseif coder.internal.isConst(sbk)
        if coder.const(sbk == ONE)
            sck = sak;
        else
            sck = sbk;
            iscompatible = iscompatible && (sak == ONE || sak == sbk);
        end
    elseif sbk == ONE
        sck = sak;
    elseif sak == ONE
        sck = sbk;
    elseif sak == sbk
        sck = sak;
    else
        % This case is a caller error, but we use sck = min(sak,sbk) for
        % safety when error checks are removed.
        iscompatible = false;
        if sbk < sak
            sck = sbk;
        else
            sck = sak;
        end
    end
    csz(k) = sck;
end
eml_invariant(iscompatible,'MATLAB:bsxfun:arrayDimensionsMustMatch');
% Determine the output type.
a0 = coder.internal.scalarEg(a);
b0 = coder.internal.scalarEg(b);
% Pass empty arrays to FUN in case FUN(A0,B0) would generate an error.
c0 = coder.internal.scalarEg(fun(a0([]),b0([])));
% Allocate the output array.
c = coder.nullcopy(eml_expand(c0,csz));
if isempty(c)
    % Nothing to compute.
elseif coder.isRowMajor
    c = applyBinaryScalarFunction(fun,a,b,c,ONE);
else
    c = applyBinaryScalarFunction(fun,a,b,c,ND);
end

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

function c = applyBinaryScalarFunction(fun,a,b,c,dim,varargin)
% Recursively generate the nested loops required to evaluate c = fun(a,b)
% in contiguous fashion.
coder.inline('always');
coder.internal.prefer_const(dim);
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
ND = coder.const(coder.internal.indexInt(eml_ndims(c)));
% Carve varargin up into 3 sets of subscripts.
NSUBSCRIPTS = coder.const(coder.internal.indexInt(length(varargin)/3));
IC1 = ONE;
IC2 = NSUBSCRIPTS;
IA1 = coder.const(NSUBSCRIPTS + ONE);
IA2 = coder.const(IA1 + NSUBSCRIPTS - ONE);
IB1 = coder.const(IA2 + ONE);
IB2 = coder.const(IB1 + NSUBSCRIPTS - ONE);
if coder.const(dim >= ONE && dim <= ND)
    % Generate a loop for this dimension and recurse.
    sza = coder.internal.indexInt(size(a,dim));
    szb = coder.internal.indexInt(size(b,dim));
    szc = coder.internal.indexInt(size(c,dim));
    % acoef == 0 --> expand this dimension of A.
    % acoef == 1 --> index along this dimension of A.
    acoef = coder.internal.indexInt(sza ~= ONE);
    % bcoef == 0 --> expand this dimension of B.
    % bcoef == 1 --> index along this dimension of B.
    bcoef = coder.internal.indexInt(szb ~= ONE);
    for k = ZERO:szc - ONE
        ia = acoef*k + ONE;
        ib = bcoef*k + ONE;
        ic = k + ONE;
        if coder.isRowMajor
            c = applyBinaryScalarFunction(fun,a,b,c,dim + ONE, ...
                varargin{IC1:IC2},ic, ...
                varargin{IA1:IA2},ia, ...
                varargin{IB1:IB2},ib);
        else
            c = applyBinaryScalarFunction(fun,a,b,c,dim - ONE, ...
                ic,varargin{IC1:IC2}, ...
                ia,varargin{IA1:IA2}, ...
                ib,varargin{IB1:IB2});
        end
    end
else
    % This is the body of the innermost loop.
    c(varargin{IC1:IC2}) = fun(a(varargin{IA1:IA2}),b(varargin{IB1:IB2}));
end

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