function y = diff(x,n,dim)
% Embedded MATLAB Library function.
% 
% Limitations:
% 1) N and DIM must be constants.
% 2) N-D arrays with N>2 are not supported.
%
% $INCLUDE(DOC) toolbox/eml/lib/matlab/diff.m $
% Copyright 2002-2004 The MathWorks, Inc.

eml_assert(nargin > 0, 'Not enough input arguments.');

if (nargin < 2)
    n = 1;
end

eml_assert(eml_is_const(n), 'Difference order N must be a constant.');
eml_assert(isPosIntScalar(n),...
    'Difference order N must be a positive integer scalar in the range 1 to intmax.');

xSize = size(x);
fixDim = nargin == 3;

if fixDim
    eml_assert(eml_is_const(dim), 'Dimension argument DIM must be a constant.');
    eml_assert(isPosIntScalar(dim),...
        'Dimension argument DIM must be a positive integer scalar in the range 1 to intmax');
    eml_assert(dim <= numel(xSize),'Unsupported return type (empty N-D matrix, N>2).');
%     if (dim > ndims(x))  %Matlab N-D behavior.
%         y = zeros([xSize,ones(1,dim-ndims(x)-1),0]);
%         return
%     end
end

%Empty matrix logic.
if isempty(x)
    if (nargin <= 2)
        dim = findFirstIndex(xSize);
        if (dim == 0) 
            y = x;
        else
            dimSize = xSize(dim);
            k = findFirstZero(xSize);
            if (k == 0) || (k > dim)
                if (dimSize > n)
                    xSize = changeDimension(xSize,dim,dimSize-n);
                    y = zeros(xSize,class(x));
                else
                    xSize = changeDimension(xSize,dim,1);
                    %This will be a shallow recursion if the number of 
                    %dimensions is small.
                    n1 = n - dimSize + 1;
                    eml_assert(n1 < n);
                    y = diff(zeros(xSize,class(x)),n1);
                end
            else
                y = x;
            end
        end
    else %decrement specified dimension
        xSize = changeDimension(xSize,dim,max(0,xSize(dim)-n));
        y = zeros(xSize,class(x));
    end
    return
end

if fixDim
    dimSize = xSize(dim);
    if (dimSize <= n)
        xSize = changeDimension(xSize,dim,0);
        y = zeros(xSize,class(x));
        return;
    end
else
    dim = findFirstIndex(xSize);
    if dim == 0
        y = [];
        return
    end
    dimSize = xSize(dim);
end

order = min(dimSize-1,n);
work = zeros(order,1,class(x));
if (dim >= 2)
    stride = prod(xSize(1:dim-1));
else
    stride = 1;
end
newDimSize = dimSize - order;
ySize = changeDimension(xSize,dim,newDimSize);
y1 = zeros(ySize,class(x));
xNext = stride*dimSize;
yNext = stride*newDimSize;
if (dim >= size(xSize,2))
    nHigh = 1;
else
    nHigh = prod(xSize(dim+1:end));
end
ixStart = 1;
iyStart = 1;
for r = 1:nHigh
    ix = ixStart;
    iy = iyStart;
    for s = 1:stride
        ixLead = ix + stride;
        iyLead = iy;
        work(1) = x(ix);
        if (order >= 2)
            for m = 1:order-1
                tmp1 = x(ixLead);
                for k = 1:m
                    tmp2 = work(k);
                    work(k) = tmp1;
                    tmp1 = tmp1 - tmp2;
                end
                work(m+1) = tmp1;
                ixLead = ixLead + stride;
            end
        end
        for m = order+1:dimSize
            tmp1 = x(ixLead);
            for k = 1:order
                tmp2 = work(k);
                work(k) = tmp1;
                tmp1 = tmp1 - tmp2;
            end
            ixLead = ixLead + stride;
            y1(iyLead) = tmp1;
            iyLead = iyLead + stride;
        end
        ix = ix + 1;
        iy = iy + 1;
    end
    ixStart = ixStart + xNext;
    iyStart = iyStart + yNext;
end

if (order < n) %Diff along another dimension.
    n2 = n - order;
    eml_assert(n2 < n);
    y = diff(y1,n2);
else
    y = y1;
end

return


function k = findFirstIndex(x)
% Find the first index and value in the vector x that is >= 2.
% This can be replaced with a call to find() when it becomes
% available in eML.
len = size(x,2);
k = 0;
for m = 1:len
    if x(m) >= 2
        k = m;
        break;
    end
end
return


function k = findFirstZero(x)
% Find the first index k such that x(k) == 0.
% This can be replaced with a call to find() when it becomes
% available in eML.
len = size(x,2);
k = 0;
for m = 1:len
    if x(m) == 0
        k = m;
        break;
    end
end
return


function y = changeDimension(xSize,dim,newDimSize)
% Circumvents a compiler limitation.
y = xSize;
y(dim) = newDimSize;
return


function p = isPosIntScalar(x)
p = false;
if (numel(x) == 1) && (x >= 1) && (x <= intmax)
    ix = int32(x);
    %Like x == floor(x), but floor doesn't work in eml_assert yet.
    p = (x == cast(ix,class(x)));
end
return