--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,11 @@
+2002-07-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* general/mod.m: Additional error checks, Texinfoize doc string.
+
+2002-04-29  Paul Kienzle <pkienzle@users.sf.net>
+
+	* general/mod.m: New function.
+
 2002-07-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* control/system/sysmult.m: Use varargin instead of old style varargs.

new file mode 100644
--- /dev/null
+++ b/scripts/general/mod.m
@@ -0,0 +1,103 @@
+## Copyright (C) 1999, 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Mapping Function} {} mod (@var{x}, @var{y})
+## Compute modulo function, using
+##
+## @example
+## x - y .* floor (x ./ y)
+## @end example
+##
+## Note that this handles negative numbers correctly:
+## @code{mod (-1, 3)} is 2, not -1 as @code{rem (-1, 3)} returns.
+## Also, @code{mod (@var{x}, 0)} returns @var{x}.
+##
+## An error message is printed if the dimensions of the arguments do not
+## agree, or if either of the arguments is complex.
+## @end deftypefn
+## @seealso{rem, round}
+
+## Author: Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
+## Modified by: Teemu Ikonen <tpikonen@pcu.helsinki.fi>
+## Adapted by: jwe
+
+function r = mod (x, y)
+
+  if (nargin != 2)
+    usage ("r = mod (x, y)");
+  endif
+
+  if (any (size (x) != size (y)) && ! (is_scalar (x) || is_scalar (y)))
+    error ("mod: argument sizes must agree");
+  endif
+
+  ## Matlab allows complex arguments, but as far as I can tell, that's a
+  ## bunch of hooey.
+
+  if (any (any (imag (x))) || any (any (imag (y))))
+    error ("mod: complex arguments are not allowed");
+  endif
+
+  nz = y != 0.0;
+  if (all (nz(:)))
+    ## No elements of y are zero.
+    r = x - y .* floor (x ./ y);
+  elseif (is_scalar (y))
+    ## y must be zero.
+    r = x;
+  else
+    ## Some elements of y are zero.
+    if (is_scalar (x))
+      r = x * ones (size (y));
+    else
+      r = x;
+      x = x(nz);
+    endif
+    r(nz) = x - y .* floor (x ./ y);
+  endif
+
+endfunction
+  
+## empty input test
+%!assert (isempty(mod([], [])));
+
+## x mod y, y != 0 tests
+%!assert (mod(5, 3), 2);
+%!assert (mod(-5, 3), 1);
+%!assert (mod(0, 3), 0);
+%!assert (mod([-5, 5, 0], [3, 3, 3]), [1, 2, 0]);
+%!assert (mod([-5; 5; 0], [3; 3; 3]), [1; 2; 0]);
+%!assert (mod([-5, 5; 0, 3], [3, 3 ; 3, 1]), [1, 2 ; 0, 0]);
+
+## x mod 0 tests
+%!assert (mod(5, 0), 5);
+%!assert (mod(-5, 0), -5);
+%!assert (mod([-5, 5, 0], [3, 0, 3]), [1, 5, 0]);
+%!assert (mod([-5; 5; 0], [3; 0; 3]), [1; 5; 0]);
+%!assert (mod([-5, 5; 0, 3], [3, 0 ; 3, 1]), [1, 5 ; 0, 0]);
+%!assert (mod([-5, 5; 0, 3], [0, 0 ; 0, 0]), [-5, 5; 0, 3]);
+
+## mixed scalar/matrix tests
+%!assert (mod([-5, 5; 0, 3], 0), [-5, 5; 0, 3]); 
+%!assert (mod([-5, 5; 0, 3], 3), [1, 2; 0, 0]);
+%!assert (mod(-5,[0,0; 0,0]), [-5, -5; -5, -5]);
+%!assert (mod(-5,[3,0; 3,1]), [1, -5; 1, 0]);
+%!assert (mod(-5,[3,2; 3,1]), [1, 1; 1, 0]);
+

--- a/scripts/general/rem.m
+++ b/scripts/general/rem.m
@@ -33,12 +33,10 @@
 
 ## Author: jwe
 
-function retval = rem (x, y)
-
-  usage_msg = "rem (x, y)";
+function r = rem (x, y)
 
   if (nargin != 2)
-    usage (usage_msg);
+    usage ("rem (x, y)");
   endif
 
   if (any (size (x) != size (y)) && ! (is_scalar (x) || is_scalar (y)))
@@ -52,10 +50,6 @@
     error ("rem: complex arguments are not allowed");
   endif
 
-  if (nargin == 2)
-    retval = x - y .* fix (x ./ y);
-  else
-    error (usage_msg);
-  endif
+  r = x - y .* fix (x ./ y);
 
 endfunction