function z = mod(x,y)
% Does the same as the MATLAB mod function, but this one is correct (hopefully)

% Check for nonreal arguments.

if ~isreal(x) | ~isreal(y)
    error('Arguments must be real.')
end

% Make arguments the same size.

if length(y) == 1
   y = y(ones(size(x)));
elseif length(x) == 1
   x = x(ones(size(y)));
end
z = zeros(size(x));

if ~isequal(size(x),size(y)), 
  error('X and Y must be the same size.'); 
end
if isempty(x) | isempty(y), return, end

% Integer denominator.
% Use the conventional formula.

m = (y == fix(y)) & (y ~= 0);
z(m) = x(m) - floor(x(m)./y(m)).*y(m);

% Noninteger denominator.
% Adjust any quotient within roundoff error of an integer.

m = (y ~= fix(y));
q = x(m)./y(m);
r = abs(q - round(q)) <= eps*abs(q);
q(r) = round(q(r));
z(m) = (q - floor(q)).*y(m);

% Zero denominator.

m = (y == 0);
z(m) = x(m);