--- a/liboctave/ChangeLog	Tue Mar 23 09:30:46 2010 +0100
+++ b/liboctave/ChangeLog	Tue Mar 23 13:01:34 2010 +0100
@@ -1,3 +1,11 @@
+2010-03-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* oct-inttypes.h (octave_int_arith_base::rem, octave_int_base::mod):
+	New methods.
+	(rem, mod): New template functions.
+	* lo-mappers.cc (rem, mod): New overloaded functions.
+	* lo-mappers.h: Declare them.
+
 2010-03-23  Jaroslav Hajek  <highegg@gmail.com>
 
 	* oct-binmap.h: New source.

--- a/liboctave/lo-mappers.cc	Tue Mar 23 09:30:46 2010 +0100
+++ b/liboctave/lo-mappers.cc	Tue Mar 23 13:01:34 2010 +0100
@@ -131,6 +131,16 @@
 }
 
 double
+mod (double x, double y)
+{
+  if (y == 0)
+    return x;
+
+  double r = fmod (x, y);
+  return ((r < 0) != (y < 0)) ? y+r : r;
+}
+
+double
 xlog2 (double x)
 {
 #if defined (HAVE_LOG2)
@@ -435,6 +445,16 @@
 }
 
 float
+mod (float x, float y)
+{
+  if (y == 0)
+    return x;
+
+  float r = fmodf (x, y);
+  return ((r < 0) != (y < 0)) ? y+r : r;
+}
+
+float
 xlog2 (float x)
 {
 #if defined (HAVE_LOG2)

--- a/liboctave/lo-mappers.h	Tue Mar 23 09:30:46 2010 +0100
+++ b/liboctave/lo-mappers.h	Tue Mar 23 13:01:34 2010 +0100
@@ -37,6 +37,7 @@
 extern OCTAVE_API double xroundb (double x);
 extern OCTAVE_API double signum (double x);
 extern OCTAVE_API double xtrunc (double x);
+extern OCTAVE_API double mod (double x, double y);
 extern OCTAVE_API double xlog2 (double x); 
 extern OCTAVE_API Complex xlog2 (const Complex& x); 
 extern OCTAVE_API double xlog2 (double x, int& exp);
@@ -125,6 +126,7 @@
 extern OCTAVE_API float xroundb (float x);
 extern OCTAVE_API float signum (float x);
 extern OCTAVE_API float xtrunc (float x);
+extern OCTAVE_API float mod (float x, float y);
 extern OCTAVE_API float xlog2 (float x); 
 extern OCTAVE_API FloatComplex xlog2 (const FloatComplex& x); 
 extern OCTAVE_API float xlog2 (float x, int& exp);

--- a/liboctave/oct-inttypes.h	Tue Mar 23 09:30:46 2010 +0100
+++ b/liboctave/oct-inttypes.h	Tue Mar 23 13:01:34 2010 +0100
@@ -386,6 +386,20 @@
           return x ? octave_int_base<T>::max_val () : 0;
         }
     }
+
+  // Remainder.
+  static T
+  rem (T x, T y)
+    {
+      return y != 0 ? x % y : 0;
+    }
+
+  // Modulus. Note the weird y = 0 case for Matlab compatibility.
+  static T
+  mod (T x, T y)
+    {
+      return y != 0 ? x % y : x;
+    }
 };
 
 #ifdef OCTAVE_INT_USE_LONG_DOUBLE
@@ -654,6 +668,25 @@
       return z;
     }
 
+  // Remainder.
+  static T
+  rem (T x, T y)
+    {
+      return y != 0 ? x % y : 0;
+    }
+
+  // Modulus. Note the weird y = 0 case for Matlab compatibility.
+  static T
+  mod (T x, T y)
+    {
+      if (y != 0)
+        {
+          T r = x % y;
+          return ((r < 0) != (y < 0)) ? r + y : r;
+        }
+      else
+        return x;
+    }
 };
 
 #ifdef OCTAVE_INT_USE_LONG_DOUBLE
@@ -781,6 +814,7 @@
   OCTAVE_INT_BIN_OP(-, sub, octave_int<T>)
   OCTAVE_INT_BIN_OP(*, mul, octave_int<T>)
   OCTAVE_INT_BIN_OP(/, div, octave_int<T>)
+  OCTAVE_INT_BIN_OP(%, rem, octave_int<T>)
   OCTAVE_INT_BIN_OP(<<, lshift, int)
   OCTAVE_INT_BIN_OP(>>, rshift, int)
 
@@ -807,6 +841,16 @@
   T ival;
 };
 
+template <class T>
+inline octave_int<T>
+rem (const octave_int<T>& x, const octave_int<T>& y)
+{ return octave_int_arith<T>::rem (x.value (), y.value ()); }
+
+template <class T>
+inline octave_int<T>
+mod (const octave_int<T>& x, const octave_int<T>& y)
+{ return octave_int_arith<T>::mod (x.value (), y.value ()); }
+
 // No mixed integer binary operations!
 
 template <class T>

--- a/scripts/general/mod.m	Tue Mar 23 09:30:46 2010 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,130 +0,0 @@
-## Copyright (C) 1999, 2000, 2002, 2004, 2005, 2006, 2007, 2008,
-##               2009 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 3 of the License, 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, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Mapping Function} {} mod (@var{x}, @var{y})
-## Compute the modulo of @var{x} and @var{y}.  Conceptually this is given by
-##
-## @example
-## x - y .* floor (x ./ y)
-## @end example
-##
-## and is written such that the correct modulus is returned for
-## integer types.  This function handles negative values correctly.  That
-## is, @code{mod (-1, 3)} is 2, not -1, as @code{rem (-1, 3)} returns.
-## @code{mod (@var{x}, 0)} returns @var{x}.
-##
-## An error results if the dimensions of the arguments do not agree, or if
-## either of the arguments is complex.
-## @seealso{rem, fmod}
-## @end deftypefn
-
-## 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)
-    print_usage ();
-  endif
-
-  if (! size_equal (x, y) && ! (isscalar (x) || isscalar (y)))
-    error ("mod: argument sizes must agree");
-  endif
-
-  if (isreal (x) && isreal (y))
-    nz = y != 0.0;
-    if (all (nz(:)))
-      ## No elements of y are zero.
-      if (isinteger(x) || isinteger(y))
-	if (isinteger (x))
-	  typ = class (x);
-	else
-	  typ = class (y);
-	endif
-	r = x - y .* cast (floor (double(x) ./ double(y)), typ);
-      else
-	r = x - y .* floor (x ./ y);
-      endif
-    elseif (isscalar (y))
-      ## y must be zero.
-      r = x;
-    else
-      ## Some elements of y are zero.
-      if (isscalar (x))
-	r = x * ones (size(y), class(y));
-      else
-	r = x;
-	x = x(nz);
-      endif
-      y = y(nz);
-      if (isinteger(x) || isinteger(y))
-	if (isinteger (x))
-	  typ = class (x);
-	else
-	  typ = class (y);
-	endif
-	r(nz) = x - y .* floor (double(x) ./ double(y));
-      else
-	r(nz) = x - y .* floor (x ./ y);
-      endif
-    endif
-  else
-    error ("mod: complex arguments are not allowed");
-  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]);
-
-## integer types
-%!assert (mod(uint8(5),uint8(4)),uint8(1))
-%!assert (mod(uint8([1:5]),uint8(4)),uint8([1,2,3,0,1]))
-%!assert (mod(uint8([1:5]),uint8(0)),uint8([1:5]))
-%!error (mod(uint8(5),int8(4)))
-
-## mixed integer/real types
-%!assert (mod(uint8(5),4),uint8(1))
-%!assert (mod(5,uint8(4)),uint8(1))
-%!assert (mod(uint8([1:5]),4),uint8([1,2,3,0,1]))

--- a/scripts/general/module.mk	Tue Mar 23 09:30:46 2010 +0100
+++ b/scripts/general/module.mk	Tue Mar 23 13:01:34 2010 +0100
@@ -50,7 +50,6 @@
   general/isvector.m \
   general/loadobj.m \
   general/logspace.m \
-  general/mod.m \
   general/nargchk.m \
   general/nargoutchk.m \
   general/nextpow2.m \
@@ -64,7 +63,6 @@
   general/quadl.m \
   general/quadv.m \
   general/rat.m \
-  general/rem.m \
   general/repmat.m \
   general/rot90.m \
   general/rotdim.m \

--- a/scripts/general/rem.m	Tue Mar 23 09:30:46 2010 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-## Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2002, 2004,
-##               2005, 2006, 2007, 2008, 2009 John W. Eaton
-##
-## 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 3 of the License, 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, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Mapping Function} {} rem (@var{x}, @var{y})
-## Return the remainder of the division @code{@var{x} / @var{y}}, computed 
-## using the expression
-##
-## @example
-## x - y .* fix (x ./ y)
-## @end example
-##
-## An error message is printed if the dimensions of the arguments do not
-## agree, or if either of the arguments is complex.
-## @seealso{mod, fmod}
-## @end deftypefn
-
-## Author: jwe
-
-function r = rem (x, y)
-
-  if (nargin != 2)
-    print_usage ();
-  endif
-
-  if (! size_equal (x, y) && ! (isscalar (x) || isscalar (y)))
-    error ("rem: argument sizes must agree");
-  endif
-
-  if (isreal (x) && isreal (y))
-      if (isinteger(x) || isinteger(y))
-	if (isinteger (x))
-	  typ = class (x);
-	else
-	  typ = class (y);
-	endif
-	r = x - y .* cast (fix (double (x) ./ double (y)), typ);
-      else
-	r = x - y .* fix (x ./ y);
-      endif
-  else
-    error ("rem: complex arguments are not allowed");
-  endif
-
-endfunction
-
-%!assert(rem ([1, 2, 3; -1, -2, -3], 2), [1, 0, 1; -1, 0, -1]);
-
-%!assert(rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3)),[1, 0, 1; -1, 0, -1]);
-
-%!error rem ();
-
-%!error rem (1, 2, 3);
-
-%!error rem ([1, 2], [3, 4, 5]);
-
-%!error rem (i, 1);
-
-%!assert(rem (uint8([1, 2, 3; -1, -2, -3]), uint8 (2)), uint8([1, 0, 1; -1, 0, -1]));
-
-%!assert(uint8(rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3))),uint8([1, 0, 1; -1, 0, -1]));
-
-%!error rem (uint(8),int8(5));
-
-%!error rem (uint8([1, 2]), uint8([3, 4, 5]));

--- a/src/ChangeLog	Tue Mar 23 09:30:46 2010 +0100
+++ b/src/ChangeLog	Tue Mar 23 13:01:34 2010 +0100
@@ -1,3 +1,8 @@
+2010-03-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* data.cc (Frem, Fmod): New DEFUNs.
+	(Ffmod): Remove. Make an alias of Frem.
+
 2010-03-23  Jaroslav Hajek  <highegg@gmail.com>
 
 	* data.cc (map_d_m, map_m_d, map_m_m, map_f_fm, map_fm_f, map_fm_fm,

--- a/src/data.cc	Tue Mar 23 09:30:46 2010 +0100
+++ b/src/data.cc	Tue Mar 23 13:01:34 2010 +0100
@@ -473,13 +473,20 @@
 %! assert (e(1:2,:), [0,1; 2,3]);
 */
 
-DEFUN (fmod, args, ,
+DEFUN (rem, args, ,
   "-*- texinfo -*-\n\
-@deftypefn {Mapping Function} {} fmod (@var{x}, @var{y})\n\
-Compute the floating point remainder of dividing @var{x} by @var{y}\n\
-using the C library function @code{fmod}.  The result has the same\n\
-sign as @var{x}.  If @var{y} is zero, the result is implementation-dependent.\n\
-@seealso{mod, rem}\n\
+@deftypefn {Mapping Function} {} rem (@var{x}, @var{y})\n\
+@deftypefnx {Mapping Function} {} fmod (@var{x}, @var{y})\n\
+Return the remainder of the division @code{@var{x} / @var{y}}, computed \n\
+using the expression\n\
+\n\
+@example\n\
+x - y .* fix (x ./ y)\n\
+@end example\n\
+\n\
+An error message is printed if the dimensions of the arguments do not\n\
+agree, or if either of the arguments is complex.\n\
+@seealso{mod}\n\
 @end deftypefn")
 {
   octave_value retval;
@@ -489,11 +496,49 @@
   if (nargin == 2)
     {
       if (! args(0).is_numeric_type ())
-        gripe_wrong_type_arg ("fmod", args(0));
+        gripe_wrong_type_arg ("rem", args(0));
       else if (! args(1).is_numeric_type ())
-        gripe_wrong_type_arg ("fmod", args(1));
+        gripe_wrong_type_arg ("rem", args(1));
       else if (args(0).is_complex_type () || args(1).is_complex_type ())
-        error ("fmod: not defined for complex numbers");
+        error ("rem: not defined for complex numbers");
+      else if (args(0).is_integer_type () || args(1).is_integer_type ())
+        {
+          builtin_type_t btyp0 = args(0).builtin_type ();
+          builtin_type_t btyp1 = args(1).builtin_type ();
+          if (btyp0 == btyp_double || btyp0 == btyp_float)
+            btyp0 = btyp1;
+          if (btyp1 == btyp_double || btyp1 == btyp_float)
+            btyp1 = btyp0;
+
+          if (btyp0 == btyp1)
+            {
+              switch (btyp0)
+                {
+#define MAKE_INT_BRANCH(X) \
+                case btyp_ ## X: \
+                    { \
+                    X##NDArray a0 = args(0).X##_array_value (); \
+                    X##NDArray a1 = args(1).X##_array_value (); \
+                    retval = binmap<octave_##X> (a0, a1, rem, "rem"); \
+                    } \
+                  break
+                MAKE_INT_BRANCH (int8);
+                MAKE_INT_BRANCH (int16);
+                MAKE_INT_BRANCH (int32);
+                MAKE_INT_BRANCH (int64);
+                MAKE_INT_BRANCH (uint8);
+                MAKE_INT_BRANCH (uint16);
+                MAKE_INT_BRANCH (uint32);
+                MAKE_INT_BRANCH (uint64);
+#undef MAKE_INT_BRANCH
+                default:
+                  panic_impossible ();
+                }
+            }
+          else
+            error ("rem: cannot combine %s and %d", 
+                   args(0).class_name ().c_str (), args(1).class_name ().c_str ());
+        }
       else if (args(0).is_single_type () || args(1).is_single_type ())
         {
           if (args(0).is_scalar_type () && args(1).is_scalar_type ())
@@ -502,7 +547,7 @@
             {
               FloatNDArray a0 = args(0).float_array_value ();
               FloatNDArray a1 = args(1).float_array_value ();
-              retval = binmap<float> (a0, a1, ::fmodf, "fmod");
+              retval = binmap<float> (a0, a1, fmodf, "rem");
             }
         }
       else
@@ -516,13 +561,13 @@
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
-              retval = binmap<double> (m0, m1, ::fmod, "fmod");
+              retval = binmap<double> (m0, m1, fmod, "rem");
             }
           else
             {
               NDArray a0 = args(0).array_value ();
               NDArray a1 = args(1).array_value ();
-              retval = binmap<double> (a0, a1, ::fmod, "fmod");
+              retval = binmap<double> (a0, a1, fmod, "rem");
             }
         }
     }
@@ -533,6 +578,21 @@
 }
 
 /*
+
+%!assert(rem ([1, 2, 3; -1, -2, -3], 2), [1, 0, 1; -1, 0, -1]);
+%!assert(rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3)),[1, 0, 1; -1, 0, -1]);
+%!error rem ();
+%!error rem (1, 2, 3);
+%!error rem ([1, 2], [3, 4, 5]);
+%!error rem (i, 1);
+%!assert(rem (uint8([1, 2, 3; -1, -2, -3]), uint8 (2)), uint8([1, 0, 1; -1, 0, -1]));
+%!assert(uint8(rem ([1, 2, 3; -1, -2, -3], 2 * ones (2, 3))),uint8([1, 0, 1; -1, 0, -1]));
+%!error rem (uint(8),int8(5));
+%!error rem (uint8([1, 2]), uint8([3, 4, 5]));
+
+*/
+
+/*
 %!assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
 %!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
 %!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
@@ -540,6 +600,154 @@
 %!assert (size (fmod (1, 2)), [1, 1])
 */
 
+DEFALIAS (fmod, rem)
+
+DEFUN (mod, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} mod (@var{x}, @var{y})\n\
+Compute the modulo of @var{x} and @var{y}.  Conceptually this is given by\n\
+\n\
+@example\n\
+x - y .* floor (x ./ y)\n\
+@end example\n\
+\n\
+and is written such that the correct modulus is returned for\n\
+integer types.  This function handles negative values correctly.  That\n\
+is, @code{mod (-1, 3)} is 2, not -1, as @code{rem (-1, 3)} returns.\n\
+@code{mod (@var{x}, 0)} returns @var{x}.\n\
+\n\
+An error results if the dimensions of the arguments do not agree, or if\n\
+either of the arguments is complex.\n\
+@seealso{rem}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 2)
+    {
+      if (! args(0).is_numeric_type ())
+        gripe_wrong_type_arg ("mod", args(0));
+      else if (! args(1).is_numeric_type ())
+        gripe_wrong_type_arg ("mod", args(1));
+      else if (args(0).is_complex_type () || args(1).is_complex_type ())
+        error ("mod: not defined for complex numbers");
+      else if (args(0).is_integer_type () || args(1).is_integer_type ())
+        {
+          builtin_type_t btyp0 = args(0).builtin_type ();
+          builtin_type_t btyp1 = args(1).builtin_type ();
+          if (btyp0 == btyp_double || btyp0 == btyp_float)
+            btyp0 = btyp1;
+          if (btyp1 == btyp_double || btyp1 == btyp_float)
+            btyp1 = btyp0;
+
+          if (btyp0 == btyp1)
+            {
+              switch (btyp0)
+                {
+#define MAKE_INT_BRANCH(X) \
+                case btyp_ ## X: \
+                    { \
+                    X##NDArray a0 = args(0).X##_array_value (); \
+                    X##NDArray a1 = args(1).X##_array_value (); \
+                    retval = binmap<octave_##X> (a0, a1, mod, "mod"); \
+                    } \
+                  break
+                MAKE_INT_BRANCH (int8);
+                MAKE_INT_BRANCH (int16);
+                MAKE_INT_BRANCH (int32);
+                MAKE_INT_BRANCH (int64);
+                MAKE_INT_BRANCH (uint8);
+                MAKE_INT_BRANCH (uint16);
+                MAKE_INT_BRANCH (uint32);
+                MAKE_INT_BRANCH (uint64);
+#undef MAKE_INT_BRANCH
+                default:
+                  panic_impossible ();
+                }
+            }
+          else
+            error ("mod: cannot combine %s and %d", 
+                   args(0).class_name ().c_str (), args(1).class_name ().c_str ());
+        }
+      else if (args(0).is_single_type () || args(1).is_single_type ())
+        {
+          if (args(0).is_scalar_type () && args(1).is_scalar_type ())
+            retval = mod (args(0).float_value (), args(1).float_value ());
+          else
+            {
+              FloatNDArray a0 = args(0).float_array_value ();
+              FloatNDArray a1 = args(1).float_array_value ();
+              retval = binmap<float> (a0, a1, mod, "mod");
+            }
+        }
+      else
+        {
+          bool a0_scalar = args(0).is_scalar_type ();
+          bool a1_scalar = args(1).is_scalar_type ();
+          if (a0_scalar && a1_scalar)
+            retval = mod (args(0).scalar_value (), args(1).scalar_value ());
+          else if ((a0_scalar || args(0).is_sparse_type ()) 
+                   && (a1_scalar || args(1).is_sparse_type ()))
+            {
+              SparseMatrix m0 = args(0).sparse_matrix_value ();
+              SparseMatrix m1 = args(1).sparse_matrix_value ();
+              retval = binmap<double> (m0, m1, mod, "mod");
+            }
+          else
+            {
+              NDArray a0 = args(0).array_value ();
+              NDArray a1 = args(1).array_value ();
+              retval = binmap<double> (a0, a1, mod, "mod");
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+## 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]);
+
+## integer types
+%!assert (mod(uint8(5),uint8(4)),uint8(1))
+%!assert (mod(uint8([1:5]),uint8(4)),uint8([1,2,3,0,1]))
+%!assert (mod(uint8([1:5]),uint8(0)),uint8([1:5]))
+%!error (mod(uint8(5),int8(4)))
+
+## mixed integer/real types
+%!assert (mod(uint8(5),4),uint8(1))
+%!assert (mod(5,uint8(4)),uint8(1))
+%!assert (mod(uint8([1:5]),4),uint8([1,2,3,0,1]))
+*/
+
 // FIXME Need to convert the reduction functions of this file for single precision
 
 #define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \