diff -r e909207e9ff1 libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/libinterp/corefcn/data.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -7926,3 +7926,519 @@
 %!error <input was not valid base64> base64_decode ("AQ=")
 %!error <incorrect input size> base64_decode ("AQ==")
 */
+
+DEFUN (movsum, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} movsum (@var{x}, @var{k})
+@deftypefnx {} {} movsum (@var{x}, @var{[kb,kf]})
+@deftypefnx {} {} movsum (@var{x}, @dots{}, @var{dim})
+Return sums of local k elements along dimension @var{dim}.
+
+If @var{dim} is omitted, it defaults to the first non-singleton dimension.
+For example:
+
+@example
+@group
+movsum ([1, 2; 3, 4; 5, 6], [1,1])
+   @result{}  4   6
+       9  12
+       8  10
+@end group
+@end example
+@seealso{movmean, movmax, movmin, cumsum}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  bool isnative = false;
+  bool isdouble = false;
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+
+  int dim = -1;
+  if (nargin == 3)
+    {
+      dim = args(2).int_value () - 1;
+      if (dim < 0)
+        error ("movsum: invalid dimension argument = %d", dim + 1);
+    }
+
+  octave_value retval;
+  octave_value arg = args(0);
+  
+  octave_idx_type kb,kf;
+  Array<octave_idx_type> n = args(1).octave_idx_type_vector_value ();
+  if (n.numel () == 1 && n(0) > 0)
+    {
+      kb = n(0) / 2;
+      kf = (n(0) - 1) / 2;
+    }
+  else if (n.numel () == 2 && n(0) >= 0 && n(0) >= 1)
+    {
+      kb = n(0);
+      kf = n(1);
+    }
+  else
+    {
+      error ("movsum: Window length must be a finite positive scalar or 2-element vector of finite nonnegative scalars");
+    }
+
+  switch (arg.builtin_type ())
+    {
+    case btyp_double:
+      if (arg.issparse ())
+        retval = arg.sparse_matrix_value ().movsum (kb, kf, dim);
+      else
+        retval = arg.array_value ().movsum (kb, kf, dim);
+      break;
+    case btyp_complex:
+      if (arg.issparse ())
+        retval = arg.sparse_complex_matrix_value ().movsum (kb, kf, dim);
+      else
+        retval = arg.complex_array_value ().movsum (kb, kf, dim);
+      break;
+    case btyp_float:
+      if (isdouble)
+        retval = arg.array_value ().movsum (kb, kf, dim);
+      else
+        retval = arg.float_array_value ().movsum (kb, kf, dim);
+      break;
+    case btyp_float_complex:
+      if (isdouble)
+        retval = arg.complex_array_value ().movsum (kb, kf, dim);
+      else
+        retval = arg.float_complex_array_value ().movsum (kb, kf, dim);
+      break;
+
+#define MAKE_INT_BRANCH(X)                                              \
+      case btyp_ ## X:                                                  \
+        if (isnative)                                                   \
+          retval = arg.X ## _array_value ().movsum (kb, kf, dim);       \
+        else                                                            \
+          retval = arg.array_value ().movsum (kb, kf, dim);             \
+        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
+
+    case btyp_bool:
+      if (arg.issparse ())
+        {
+          SparseMatrix cs = arg.sparse_matrix_value ().movsum (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      else
+        {
+          NDArray cs = arg.array_value ().movsum (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      break;
+
+    default:
+      err_wrong_type_arg ("movsum", arg);
+    }
+
+  return retval;
+}
+
+DEFUN (movprod, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} movprod (@var{x}, @var{k})
+@deftypefnx {} {} movprod (@var{x}, @var{[kb,kf]})
+@deftypefnx {} {} movprod (@var{x}, @dots{}, @var{dim})
+Return products of local k points along dimension @var{dim}.
+
+If @var{dim} is omitted, it defaults to the first non-singleton dimension.
+For example:
+
+@example
+@group
+movprod ([1, 2; 3, 4; 5, 6], [1,1])
+   @result{}  3   8
+      15  48
+      15  24
+@end group
+@end example
+@seealso{movmean, movmax, movmin, cumprod}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  bool isnative = false;
+  bool isdouble = false;
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+
+  int dim = -1;
+  if (nargin == 3)
+    {
+      dim = args(2).int_value () - 1;
+      if (dim < 0)
+        error ("movprod: invalid dimension argument = %d", dim + 1);
+    }
+
+  octave_value retval;
+  octave_value arg = args(0);
+
+  octave_idx_type kb,kf;
+  Array<octave_idx_type> n = args(1).octave_idx_type_vector_value ();
+  if (n.numel () == 1 && n(0) > 0)
+    {
+      kb = n(0) / 2;
+      kf = (n(0) - 1) / 2;
+    }
+  else if (n.numel () == 2 && n(0) >= 0 && n(0) >= 1)
+    {
+      kb = n(0);
+      kf = n(1);
+    }
+  else
+    {
+      error ("movprod: Window length must be a finite positive scalar or 2-element vector of finite nonnegative scalars");
+    }
+
+  switch (arg.builtin_type ())
+    {
+    case btyp_double:
+      if (arg.issparse ())
+        retval = arg.sparse_matrix_value ().movprod (kb, kf, dim);
+      else
+        retval = arg.array_value ().movprod (kb, kf, dim);
+      break;
+    case btyp_complex:
+      if (arg.issparse ())
+        retval = arg.sparse_complex_matrix_value ().movprod (kb, kf, dim);
+      else
+        retval = arg.complex_array_value ().movprod (kb, kf, dim);
+      break;
+    case btyp_float:
+      if (isdouble)
+        retval = arg.array_value ().movprod (kb, kf, dim);
+      else
+        retval = arg.float_array_value ().movprod (kb, kf, dim);
+      break;
+    case btyp_float_complex:
+      if (isdouble)
+        retval = arg.complex_array_value ().movprod (kb, kf, dim);
+      else
+        retval = arg.float_complex_array_value ().movprod (kb, kf, dim);
+      break;
+
+#define MAKE_INT_BRANCH(X)                                               \
+      case btyp_ ## X:                                                   \
+        if (isnative)                                                    \
+          retval = arg.X ## _array_value ().movprod (kb, kf, dim);       \
+        else                                                             \
+          retval = arg.array_value ().movprod (kb, kf, dim);             \
+        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
+
+    case btyp_bool:
+      if (arg.issparse ())
+        {
+          SparseMatrix cs = arg.sparse_matrix_value ().movprod (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      else
+        {
+          NDArray cs = arg.array_value ().movprod (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      break;
+
+    default:
+      err_wrong_type_arg ("movprod", arg);
+    }
+
+  return retval;
+}
+
+DEFUN (movmax, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} movmax (@var{x}, @var{k})
+@deftypefnx {} {} movmax (@var{x}, @var{[kb,kf]})
+@deftypefnx {} {} movmax (@var{x}, @dots{}, @var{dim})
+Return products of local k points along dimension @var{dim}.
+
+If @var{dim} is omitted, it defaults to the first non-singleton dimension.
+For example:
+
+@example
+@group
+movmax ([1, 2; 3, 4; 5, 6], [1,1])
+   @result{}  3   4
+       5   4
+       5   6
+@end group
+@end example
+@seealso{movmean, movprod, movmin, cumprod}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  bool isnative = false;
+  bool isdouble = false;
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+
+  int dim = -1;
+  if (nargin == 3)
+    {
+      dim = args(2).int_value () - 1;
+      if (dim < 0)
+        error ("movmax: invalid dimension argument = %d", dim + 1);
+    }
+
+  octave_value retval;
+  octave_value arg = args(0);
+
+  octave_idx_type kb,kf;
+  Array<octave_idx_type> n = args(1).octave_idx_type_vector_value ();
+  if (n.numel () == 1 && n(0) > 0)
+    {
+      kb = n(0) / 2;
+      kf = (n(0) - 1) / 2;
+    }
+  else if (n.numel () == 2 && n(0) >= 0 && n(0) >= 1)
+    {
+      kb = n(0);
+      kf = n(1);
+    }
+  else
+    {
+      error ("movmax: Window length must be a finite positive scalar or 2-element vector of finite nonnegative scalars");
+    }
+
+  switch (arg.builtin_type ())
+    {
+    case btyp_double:
+      if (arg.issparse ())
+        retval = arg.sparse_matrix_value ().movmax (kb, kf, dim);
+      else
+        retval = arg.array_value ().movmax (kb, kf, dim);
+      break;
+    case btyp_complex:
+      if (arg.issparse ())
+        retval = arg.sparse_complex_matrix_value ().movmax (kb, kf, dim);
+      else
+        retval = arg.complex_array_value ().movmax (kb, kf, dim);
+      break;
+    case btyp_float:
+      if (isdouble)
+        retval = arg.array_value ().movmax (kb, kf, dim);
+      else
+        retval = arg.float_array_value ().movmax (kb, kf, dim);
+      break;
+    case btyp_float_complex:
+      if (isdouble)
+        retval = arg.complex_array_value ().movmax (kb, kf, dim);
+      else
+        retval = arg.float_complex_array_value ().movmax (kb, kf, dim);
+      break;
+
+#define MAKE_INT_BRANCH(X)                                               \
+      case btyp_ ## X:                                                   \
+        if (isnative)                                                    \
+          retval = arg.X ## _array_value ().movmax (kb, kf, dim);       \
+        else                                                             \
+          retval = arg.array_value ().movmax (kb, kf, dim);             \
+        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
+
+    case btyp_bool:
+      if (arg.issparse ())
+        {
+          SparseMatrix cs = arg.sparse_matrix_value ().movmax (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      else
+        {
+          NDArray cs = arg.array_value ().movmax (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      break;
+
+    default:
+      err_wrong_type_arg ("movmax", arg);
+    }
+
+  return retval;
+}
+
+DEFUN (movmin, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} movmin (@var{x}, @var{k})
+@deftypefnx {} {} movmin (@var{x}, @var{[kb,kf]})
+@deftypefnx {} {} movmin (@var{x}, @dots{}, @var{dim})
+Return products of local k points along dimension @var{dim}.
+
+If @var{dim} is omitted, it defaults to the first non-singleton dimension.
+For example:
+
+@example
+@group
+movmin ([1, 2; 3, 4; 5, 6], [1,1])
+   @result{}  1   2
+       1   2
+       3   4
+@end group
+@end example
+@seealso{movmean, movprod, movsum, cumprod}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  bool isnative = false;
+  bool isdouble = false;
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+
+  int dim = -1;
+  if (nargin == 3)
+    {
+      dim = args(2).int_value () - 1;
+      if (dim < 0)
+        error ("movmin: invalid dimension argument = %d", dim + 1);
+    }
+
+  octave_value retval;
+  octave_value arg = args(0);
+
+  octave_idx_type kb,kf;
+  Array<octave_idx_type> n = args(1).octave_idx_type_vector_value ();
+  if (n.numel () == 1 && n(0) > 0)
+    {
+      kb = n(0) / 2;
+      kf = (n(0) - 1) / 2;
+    }
+  else if (n.numel () == 2 && n(0) >= 0 && n(0) >= 1)
+    {
+      kb = n(0);
+      kf = n(1);
+    }
+  else
+    {
+      error ("movmin: Window length must be a finite positive scalar or 2-element vector of finite nonnegative scalars");
+    }
+
+  switch (arg.builtin_type ())
+    {
+    case btyp_double:
+      if (arg.issparse ())
+        retval = arg.sparse_matrix_value ().movmin (kb, kf, dim);
+      else
+        retval = arg.array_value ().movmin (kb, kf, dim);
+      break;
+    case btyp_complex:
+      if (arg.issparse ())
+        retval = arg.sparse_complex_matrix_value ().movmin (kb, kf, dim);
+      else
+        retval = arg.complex_array_value ().movmin (kb, kf, dim);
+      break;
+    case btyp_float:
+      if (isdouble)
+        retval = arg.array_value ().movmin (kb, kf, dim);
+      else
+        retval = arg.float_array_value ().movmin (kb, kf, dim);
+      break;
+    case btyp_float_complex:
+      if (isdouble)
+        retval = arg.complex_array_value ().movmin (kb, kf, dim);
+      else
+        retval = arg.float_complex_array_value ().movmin (kb, kf, dim);
+      break;
+
+#define MAKE_INT_BRANCH(X)                                               \
+      case btyp_ ## X:                                                   \
+        if (isnative)                                                    \
+          retval = arg.X ## _array_value ().movmin (kb, kf, dim);       \
+        else                                                             \
+          retval = arg.array_value ().movmin (kb, kf, dim);             \
+        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
+
+    case btyp_bool:
+      if (arg.issparse ())
+        {
+          SparseMatrix cs = arg.sparse_matrix_value ().movmin (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      else
+        {
+          NDArray cs = arg.array_value ().movmin (kb, kf, dim);
+          if (isnative)
+            retval = cs != 0.0;
+          else
+            retval = cs;
+        }
+      break;
+
+    default:
+      err_wrong_type_arg ("movmin", arg);
+    }
+
+  return retval;
+}
diff -r e909207e9ff1 liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CMatrix.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -2996,6 +2996,30 @@
 }
 
 ComplexMatrix
+ComplexMatrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return ComplexNDArray::movsum (kb, kf, dim);
+}
+
+ComplexMatrix
+ComplexMatrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return ComplexNDArray::movprod (kb, kf, dim);
+}
+
+ComplexMatrix
+ComplexMatrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return ComplexNDArray::movmax (kb, kf, dim);
+}
+
+ComplexMatrix
+ComplexMatrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return ComplexNDArray::movmin (kb, kf, dim);
+}
+
+ComplexMatrix
 ComplexMatrix::diag (octave_idx_type k) const
 {
   return ComplexNDArray::diag (k);
diff -r e909207e9ff1 liboctave/array/CMatrix.h
--- a/liboctave/array/CMatrix.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CMatrix.h	Mon Apr 02 03:59:40 2018 +0800
@@ -388,6 +388,11 @@
   ComplexMatrix sumsq (int dim = -1) const;
   Matrix abs (void) const;
 
+  ComplexMatrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexMatrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexMatrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexMatrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   ComplexMatrix diag (octave_idx_type k = 0) const;
 
   ComplexDiagMatrix diag (octave_idx_type m, octave_idx_type n) const;
diff -r e909207e9ff1 liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CNDArray.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -632,6 +632,30 @@
 }
 
 ComplexNDArray
+ComplexNDArray::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<Complex, Complex> (*this, dim, kb, kf, mx_inline_movsum);
+}
+
+ComplexNDArray
+ComplexNDArray::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<Complex, Complex> (*this, dim, kb, kf, mx_inline_movprod);
+}
+
+ComplexNDArray
+ComplexNDArray::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<Complex, Complex> (*this, dim, kb, kf, mx_inline_movmax);
+}
+
+ComplexNDArray
+ComplexNDArray::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<Complex, Complex> (*this, dim, kb, kf, mx_inline_movmin);
+}
+
+ComplexNDArray
 ComplexNDArray::concat (const ComplexNDArray& rb,
                         const Array<octave_idx_type>& ra_idx)
 {
diff -r e909207e9ff1 liboctave/array/CNDArray.h
--- a/liboctave/array/CNDArray.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CNDArray.h	Mon Apr 02 03:59:40 2018 +0800
@@ -97,6 +97,11 @@
 
   ComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
+  ComplexNDArray movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexNDArray movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexNDArray movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  ComplexNDArray movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   ComplexNDArray& insert (const NDArray& a,
                           octave_idx_type r, octave_idx_type c);
   ComplexNDArray& insert (const ComplexNDArray& a,
diff -r e909207e9ff1 liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CSparse.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -7465,6 +7465,30 @@
   return MSparse<Complex>::diag (k);
 }
 
+SparseComplexMatrix
+SparseComplexMatrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVSUM (SparseComplexMatrix, Complex, movsum);
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVPROD (SparseComplexMatrix, Complex, movprod);
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVMAX (SparseComplexMatrix, Complex, movmax);
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVMIN (SparseComplexMatrix, Complex, movmin);
+}
+
 std::ostream&
 operator << (std::ostream& os, const SparseComplexMatrix& a)
 {
diff -r e909207e9ff1 liboctave/array/CSparse.h
--- a/liboctave/array/CSparse.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/CSparse.h	Mon Apr 02 03:59:40 2018 +0800
@@ -445,6 +445,11 @@
 
   SparseComplexMatrix diag (octave_idx_type k = 0) const;
 
+  SparseComplexMatrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseComplexMatrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseComplexMatrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseComplexMatrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   // i/o
   friend OCTAVE_API std::ostream& operator << (std::ostream& os,
                                                const SparseComplexMatrix& a);
diff -r e909207e9ff1 liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dMatrix.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -2580,6 +2580,30 @@
 }
 
 Matrix
+Matrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return NDArray::movsum (kb, kf, dim);
+}
+
+Matrix
+Matrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return NDArray::movprod (kb, kf, dim);
+}
+
+Matrix
+Matrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return NDArray::movmax (kb, kf, dim);
+}
+
+Matrix
+Matrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return NDArray::movmin (kb, kf, dim);
+}
+
+Matrix
 Matrix::diag (octave_idx_type k) const
 {
   return NDArray::diag (k);
diff -r e909207e9ff1 liboctave/array/dMatrix.h
--- a/liboctave/array/dMatrix.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dMatrix.h	Mon Apr 02 03:59:40 2018 +0800
@@ -332,6 +332,11 @@
   Matrix sumsq (int dim = -1) const;
   Matrix abs (void) const;
 
+  Matrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  Matrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  Matrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  Matrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   Matrix diag (octave_idx_type k = 0) const;
 
   DiagMatrix diag (octave_idx_type m, octave_idx_type n) const;
diff -r e909207e9ff1 liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dNDArray.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -729,6 +729,30 @@
 }
 
 NDArray
+NDArray::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<double, double> (*this, dim, kb, kf, mx_inline_movsum);
+}
+
+NDArray
+NDArray::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<double, double> (*this, dim, kb, kf, mx_inline_movprod);
+}
+
+NDArray
+NDArray::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<double, double> (*this, dim, kb, kf, mx_inline_movmax);
+}
+
+NDArray
+NDArray::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<double, double> (*this, dim, kb, kf, mx_inline_movmin);
+}
+
+NDArray
 NDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx)
 {
   if (rb.numel () > 0)
diff -r e909207e9ff1 liboctave/array/dNDArray.h
--- a/liboctave/array/dNDArray.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dNDArray.h	Mon Apr 02 03:59:40 2018 +0800
@@ -116,6 +116,11 @@
 
   NDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
+  NDArray movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  NDArray movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  NDArray movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  NDArray movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   NDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   NDArray& insert (const NDArray& a, const Array<octave_idx_type>& ra_idx);
 
diff -r e909207e9ff1 liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dSparse.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -7469,6 +7469,30 @@
 }
 
 SparseMatrix
+SparseMatrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVSUM (SparseMatrix, double, movsum);
+}
+
+SparseMatrix
+SparseMatrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVPROD (SparseMatrix, double, movprod);
+}
+
+SparseMatrix
+SparseMatrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVMAX (SparseMatrix, double, movmax);
+}
+
+SparseMatrix
+SparseMatrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  SPARSE_MOVMIN (SparseMatrix, double, movmin);
+}
+
+SparseMatrix
 SparseMatrix::diag (octave_idx_type k) const
 {
   return MSparse<double>::diag (k);
diff -r e909207e9ff1 liboctave/array/dSparse.h
--- a/liboctave/array/dSparse.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/dSparse.h	Mon Apr 02 03:59:40 2018 +0800
@@ -418,6 +418,11 @@
   SparseMatrix sumsq (int dim = -1) const;
   SparseMatrix abs (void) const;
 
+  SparseMatrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseMatrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseMatrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  SparseMatrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   SparseMatrix diag (octave_idx_type k = 0) const;
 
   Matrix matrix_value (void) const;
diff -r e909207e9ff1 liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fCMatrix.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -3019,6 +3019,30 @@
 }
 
 FloatComplexMatrix
+FloatComplexMatrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatComplexNDArray::movsum (kb, kf, dim);
+}
+
+FloatComplexMatrix
+FloatComplexMatrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatComplexNDArray::movprod (kb, kf, dim);
+}
+
+FloatComplexMatrix
+FloatComplexMatrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatComplexNDArray::movmax (kb, kf, dim);
+}
+
+FloatComplexMatrix
+FloatComplexMatrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatComplexNDArray::movmin (kb, kf, dim);
+}
+
+FloatComplexMatrix
 FloatComplexMatrix::diag (octave_idx_type k) const
 {
   return FloatComplexNDArray::diag (k);
diff -r e909207e9ff1 liboctave/array/fCMatrix.h
--- a/liboctave/array/fCMatrix.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fCMatrix.h	Mon Apr 02 03:59:40 2018 +0800
@@ -410,6 +410,11 @@
   FloatComplexMatrix sumsq (int dim = -1) const;
   FloatMatrix abs (void) const;
 
+  FloatComplexMatrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexMatrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexMatrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexMatrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   FloatComplexMatrix diag (octave_idx_type k = 0) const;
 
   FloatComplexDiagMatrix diag (octave_idx_type m, octave_idx_type n) const;
diff -r e909207e9ff1 liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fCNDArray.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -634,6 +634,30 @@
 }
 
 FloatComplexNDArray
+FloatComplexNDArray::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<FloatComplex, FloatComplex> (*this, dim, kb, kf, mx_inline_movsum);
+}
+
+FloatComplexNDArray
+FloatComplexNDArray::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<FloatComplex, FloatComplex> (*this, dim, kb, kf, mx_inline_movprod);
+}
+
+FloatComplexNDArray
+FloatComplexNDArray::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<FloatComplex, FloatComplex> (*this, dim, kb, kf, mx_inline_movmax);
+}
+
+FloatComplexNDArray
+FloatComplexNDArray::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<FloatComplex, FloatComplex> (*this, dim, kb, kf, mx_inline_movmin);
+}
+
+FloatComplexNDArray
 FloatComplexNDArray::concat (const FloatComplexNDArray& rb,
                              const Array<octave_idx_type>& ra_idx)
 {
diff -r e909207e9ff1 liboctave/array/fCNDArray.h
--- a/liboctave/array/fCNDArray.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fCNDArray.h	Mon Apr 02 03:59:40 2018 +0800
@@ -101,6 +101,12 @@
 
   FloatComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
+  FloatComplexNDArray movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexNDArray movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexNDArray movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatComplexNDArray movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
+
   FloatComplexNDArray& insert (const NDArray& a,
                                octave_idx_type r, octave_idx_type c);
   FloatComplexNDArray& insert (const FloatComplexNDArray& a,
diff -r e909207e9ff1 liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fMatrix.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -2589,6 +2589,30 @@
 }
 
 FloatMatrix
+FloatMatrix::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatNDArray::movsum (kb, kf, dim);
+}
+
+FloatMatrix
+FloatMatrix::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatNDArray::movprod (kb, kf, dim);
+}
+
+FloatMatrix
+FloatMatrix::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatNDArray::movmax (kb, kf, dim);
+}
+
+FloatMatrix
+FloatMatrix::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return FloatNDArray::movmin (kb, kf, dim);
+}
+
+FloatMatrix
 FloatMatrix::diag (octave_idx_type k) const
 {
   return FloatNDArray::diag (k);
diff -r e909207e9ff1 liboctave/array/fMatrix.h
--- a/liboctave/array/fMatrix.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fMatrix.h	Mon Apr 02 03:59:40 2018 +0800
@@ -346,6 +346,11 @@
   FloatMatrix sumsq (int dim = -1) const;
   FloatMatrix abs (void) const;
 
+  FloatMatrix movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatMatrix movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatMatrix movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatMatrix movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   FloatMatrix diag (octave_idx_type k = 0) const;
 
   FloatDiagMatrix diag (octave_idx_type m, octave_idx_type n) const;
diff -r e909207e9ff1 liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fNDArray.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -687,6 +687,30 @@
 }
 
 FloatNDArray
+FloatNDArray::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<float, float> (*this, dim, kb, kf, mx_inline_movsum);
+}
+
+FloatNDArray
+FloatNDArray::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<float, float> (*this, dim, kb, kf, mx_inline_movprod);
+}
+
+FloatNDArray
+FloatNDArray::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<float, float> (*this, dim, kb, kf, mx_inline_movmax);
+}
+
+FloatNDArray
+FloatNDArray::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<float, float> (*this, dim, kb, kf, mx_inline_movmin);
+}
+
+FloatNDArray
 FloatNDArray::concat (const FloatNDArray& rb,
                       const Array<octave_idx_type>& ra_idx)
 {
diff -r e909207e9ff1 liboctave/array/fNDArray.h
--- a/liboctave/array/fNDArray.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/fNDArray.h	Mon Apr 02 03:59:40 2018 +0800
@@ -110,6 +110,11 @@
 
   FloatNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
+  FloatNDArray movsum (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatNDArray movprod (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatNDArray movmax (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+  FloatNDArray movmin (octave_idx_type kb, octave_idx_type kf, int dim = -1) const;
+
   FloatNDArray& insert (const FloatNDArray& a,
                         octave_idx_type r, octave_idx_type c);
   FloatNDArray& insert (const FloatNDArray& a,
diff -r e909207e9ff1 liboctave/array/intNDArray.cc
--- a/liboctave/array/intNDArray.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/intNDArray.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -298,3 +298,31 @@
 {
   return do_mx_diff_op<T> (*this, dim, order, mx_inline_diff);
 }
+
+template <typename T>
+intNDArray<T>
+intNDArray<T>::movsum (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<T, T> (*this, dim, kb, kf, mx_inline_movsum);
+}
+
+template <typename T>
+intNDArray<T>
+intNDArray<T>::movprod (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<T, T> (*this, dim, kb, kf, mx_inline_movprod);
+}
+
+template <typename T>
+intNDArray<T>
+intNDArray<T>::movmax (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<T, T> (*this, dim, kb, kf, mx_inline_movmax);
+}
+
+template <typename T>
+intNDArray<T>
+intNDArray<T>::movmin (octave_idx_type kb, octave_idx_type kf, int dim) const
+{
+  return do_mx_mov_op<T, T> (*this, dim, kb, kf, mx_inline_movmin);
+}
diff -r e909207e9ff1 liboctave/array/intNDArray.h
--- a/liboctave/array/intNDArray.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/array/intNDArray.h	Mon Apr 02 03:59:40 2018 +0800
@@ -99,6 +99,11 @@
 
   intNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
+  intNDArray movsum (octave_idx_type kb, octave_idx_type kf, int dim) const;
+  intNDArray movprod (octave_idx_type kb, octave_idx_type kf, int dim) const;
+  intNDArray movmax (octave_idx_type kb, octave_idx_type kf, int dim) const;
+  intNDArray movmin (octave_idx_type kb, octave_idx_type kf, int dim) const;
+
   intNDArray abs (void) const;
   intNDArray signum (void) const;
 
diff -r e909207e9ff1 liboctave/operators/Sparse-op-defs.h
--- a/liboctave/operators/Sparse-op-defs.h	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/operators/Sparse-op-defs.h	Mon Apr 02 03:59:40 2018 +0800
@@ -1524,6 +1524,322 @@
                                                                         \
   return retval
 
+#define SPARSE_MOVSUM(RET_TYPE, ELT_TYPE, FCN)                                \
+                                                                              \
+  octave_idx_type nr = rows ();                                               \
+  octave_idx_type nc = cols ();                                               \
+                                                                              \
+  RET_TYPE retval;                                                            \
+                                                                              \
+  if (nr > 0 && nc > 0)                                                       \
+    {                                                                         \
+      if ((nr == 1 && dim == -1) || dim == 1)                                 \
+        /* Ugly!! Is there a better way? */                                   \
+        retval = transpose (). FCN (kb, kf, 0) .transpose ();                 \
+      else                                                                    \
+        {                                                                     \
+          octave_idx_type nel = 0;                                            \
+          for (octave_idx_type i = 0; i < nc; i++)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; j++)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kf;                          \
+                  octave_idx_type clamp_kf = j + kb + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          ++nel;                                              \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                }                                                             \
+            }                                                                 \
+                                                                              \
+          retval = RET_TYPE(nr, nc, nel);                                     \
+          retval.cidx (0) = 0;                                                \
+          octave_idx_type ii = 0;                                             \
+          for (octave_idx_type i = 0; i < nc; ++i)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; ++j)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kf;                          \
+                  octave_idx_type clamp_kf = j + kb + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  bool find_flag = false;                                     \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if (!find_flag)                                     \
+                            {                                                 \
+                              retval.data(ii) = ELT_TYPE();                   \
+                              retval.ridx(ii) = j;                            \
+                            }                                                 \
+                          retval.data(ii) += data(k);                         \
+                          find_flag = true;                                   \
+                        }                                                     \
+                    }                                                         \
+                  if (find_flag) ++ii;                                        \
+                }                                                             \
+              retval.cidx (i + 1) = ii;                                       \
+            }                                                                 \
+          }                                                                   \
+    }                                                                         \
+  else                                                                        \
+    retval = RET_TYPE (nr,nc);                                                \
+                                                                              \
+  return retval
+
+#define SPARSE_MOVPROD(RET_TYPE, ELT_TYPE, FCN)                               \
+                                                                              \
+  octave_idx_type nr = rows ();                                               \
+  octave_idx_type nc = cols ();                                               \
+                                                                              \
+  RET_TYPE retval;                                                            \
+                                                                              \
+  if (nr > 0 && nc > 0)                                                       \
+    {                                                                         \
+      if ((nr == 1 && dim == -1) || dim == 1)                                 \
+        /* Ugly!! Is there a better way? */                                   \
+        retval = transpose (). FCN (kb, kf, 0) .transpose ();                 \
+      else                                                                    \
+        {                                                                     \
+          octave_idx_type nel = 0;                                            \
+          for (octave_idx_type i = 0; i < nc; i++)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; j++)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kf;                          \
+                  octave_idx_type clamp_kf = j + kb + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  bool zero_flag = false;                                     \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) < clamp_kb || ridx(k) >= clamp_kf)          \
+                        {                                                     \
+                          zero_flag = true;                                   \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                  if (!zero_flag) ++nel;                                      \
+                }                                                             \
+            }                                                                 \
+                                                                              \
+          retval = RET_TYPE(nr, nc, nel);                                     \
+          retval.cidx (0) = 0;                                                \
+          octave_idx_type ii = 0;                                             \
+          for (octave_idx_type i = 0; i < nc; ++i)                            \
+            {                                                                 \
+              for (octave_idx_type j = cidx(i); j < cidx(i + 1); ++j)         \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kf;                          \
+                  octave_idx_type clamp_kf = j + kb + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  bool zero_flag = false;                                     \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if (!zero_flag)                                     \
+                            {                                                 \
+                              retval.data(ii) = ELT_TYPE(1.);                 \
+                              retval.ridx(ii) = j;                            \
+                            }                                                 \
+                          retval.data(ii) *= data(k);                         \
+                        }                                                     \
+                      else                                                    \
+                        {                                                     \
+                          zero_flag = true;                                   \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                  if (!zero_flag) ++ii;                                       \
+                }                                                             \
+              retval.cidx (i + 1) = ii;                                       \
+            }                                                                 \
+          }                                                                   \
+    }                                                                         \
+  else                                                                        \
+    retval = RET_TYPE (nr,nc);                                                \
+                                                                              \
+  return retval
+
+#define SPARSE_MOVMAX(RET_TYPE, ELT_TYPE, FCN)                                \
+                                                                              \
+  octave_idx_type nr = rows ();                                               \
+  octave_idx_type nc = cols ();                                               \
+                                                                              \
+  RET_TYPE retval;                                                            \
+                                                                              \
+  if (nr > 0 && nc > 0)                                                       \
+    {                                                                         \
+      if ((nr == 1 && dim == -1) || dim == 1)                                 \
+        /* Ugly!! Is there a better way? */                                   \
+        retval = transpose (). FCN (kb, kf, 0) .transpose ();                 \
+      else                                                                    \
+        {                                                                     \
+          octave_idx_type nel = 0;                                            \
+          for (octave_idx_type i = 0; i < nc; i++)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; j++)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kb;                          \
+                  octave_idx_type clamp_kf = j + kf + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  bool greater_zero_flag = false;                             \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if(data(k) > ELT_TYPE()) greater_zero_flag = true;  \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                  if (greater_zero_flag ||                                    \
+                    cidx(i + 1) - cidx(i) == nr) ++nel;                       \
+                }                                                             \
+            }                                                                 \
+                                                                              \
+          retval = RET_TYPE(nr, nc, nel);                                     \
+          retval.cidx (0) = 0;                                                \
+          octave_idx_type ii = 0;                                             \
+          for (octave_idx_type i = 0; i < nc; ++i)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; ++j)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kb;                          \
+                  octave_idx_type clamp_kf = j + kf + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  ELT_TYPE max = ELT_TYPE();                                  \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          max = data(k);                                      \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                                                                              \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if (data(k) > max) max = data(k);                   \
+                        }                                                     \
+                    }                                                         \
+                  if (max != ELT_TYPE())                                      \
+                    {                                                         \
+                      retval.data(ii) = max;                                  \
+                      retval.ridx(ii) = j;                                    \
+                      ++ii;                                                   \
+                    }                                                         \
+                }                                                             \
+              retval.cidx (i + 1) = ii;                                       \
+            }                                                                 \
+          }                                                                   \
+    }                                                                         \
+  else                                                                        \
+    retval = RET_TYPE (nr,nc);                                                \
+                                                                              \
+  return retval
+
+#define SPARSE_MOVMIN(RET_TYPE, ELT_TYPE, FCN)                                \
+                                                                              \
+  octave_idx_type nr = rows ();                                               \
+  octave_idx_type nc = cols ();                                               \
+                                                                              \
+  RET_TYPE retval;                                                            \
+                                                                              \
+  if (nr > 0 && nc > 0)                                                       \
+    {                                                                         \
+      if ((nr == 1 && dim == -1) || dim == 1)                                 \
+        /* Ugly!! Is there a better way? */                                   \
+        retval = transpose (). FCN (kb, kf, 0) .transpose ();                 \
+      else                                                                    \
+        {                                                                     \
+          octave_idx_type nel = 0;                                            \
+          for (octave_idx_type i = 0; i < nc; i++)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; j++)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kb;                          \
+                  octave_idx_type clamp_kf = j + kf + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  bool lesser_zero_flag = false;                              \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if(data(k) < ELT_TYPE()) lesser_zero_flag = true;   \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                  if (lesser_zero_flag ||                                     \
+                    cidx(i + 1) - cidx(i) == nr) ++nel;                       \
+                }                                                             \
+            }                                                                 \
+                                                                              \
+          retval = RET_TYPE(nr, nc, nel);                                     \
+          retval.cidx (0) = 0;                                                \
+          octave_idx_type ii = 0;                                             \
+          for (octave_idx_type i = 0; i < nc; ++i)                            \
+            {                                                                 \
+              for (octave_idx_type j = 0; j < nr; ++j)                        \
+                {                                                             \
+                  octave_idx_type clamp_kb = j - kb;                          \
+                  octave_idx_type clamp_kf = j + kf + 1;                      \
+                  if (clamp_kb < 0) clamp_kb = 0;                             \
+                  if (clamp_kf > nr) clamp_kf = nr;                           \
+                                                                              \
+                  ELT_TYPE min = ELT_TYPE();                                  \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          min = data(k);                                      \
+                          break;                                              \
+                        }                                                     \
+                    }                                                         \
+                                                                              \
+                  for (octave_idx_type k = cidx(i); k < cidx(i + 1) ;++k)     \
+                    {                                                         \
+                      if (ridx(k) >= clamp_kb && ridx(k) < clamp_kf)          \
+                        {                                                     \
+                          if (data(k) < min) min = data(k);                   \
+                        }                                                     \
+                    }                                                         \
+                  if (min != ELT_TYPE())                                      \
+                    {                                                         \
+                      retval.data(ii) = min;                                  \
+                      retval.ridx(ii) = j;                                    \
+                      ++ii;                                                   \
+                    }                                                         \
+                }                                                             \
+              retval.cidx (i + 1) = ii;                                       \
+            }                                                                 \
+          }                                                                   \
+    }                                                                         \
+  else                                                                        \
+    retval = RET_TYPE (nr,nc);                                                \
+                                                                              \
+  return retval
+
 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \
                                  INIT_VAL, MT_RESULT)                   \
                                                                         \
diff -r e909207e9ff1 liboctave/operators/mx-inlines.cc
--- a/liboctave/operators/mx-inlines.cc	Fri Mar 30 17:47:12 2018 -0400
+++ b/liboctave/operators/mx-inlines.cc	Mon Apr 02 03:59:40 2018 +0800
@@ -1500,6 +1500,126 @@
     }
 }
 
+// movsum
+template <typename TSRC, typename TRES>
+  inline void
+  mx_inline_movsum (const TSRC *v, TRES *r,
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,
+     octave_idx_type kb, octave_idx_type kf)
+  {
+    for (octave_idx_type i = 0; i < u; i++)
+      {
+        const TSRC* dim_begin = v;
+        const TSRC* dim_end = v + n * l;
+        for(octave_idx_type j = 0; j < n; j++)
+          {
+            for(octave_idx_type k = 0; k < l; k++)
+              {
+                const TSRC* clamp_kb = v - kb * l;
+                const TSRC* clamp_kf = v + (kf + 1) * l;
+                if (clamp_kb < dim_begin + k) clamp_kb = dim_begin + k;
+                if (clamp_kf > dim_end + k) clamp_kf = dim_end + k;
+
+                *r = *clamp_kb;
+                for (const TSRC* h = clamp_kb + l; h < clamp_kf; h += l)
+                    *r += *h;
+                ++r;
+                ++v;
+              }
+          }
+      }
+  }
+
+// movprod
+template <typename TSRC, typename TRES>
+  inline void
+  mx_inline_movprod (const TSRC *v, TRES *r,
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,
+     octave_idx_type kb, octave_idx_type kf)
+  {
+    for (octave_idx_type i = 0; i < u; i++)
+      {
+        const TSRC* dim_begin = v;
+        const TSRC* dim_end = v + n * l;
+        for(octave_idx_type j = 0; j < n; j++)
+          {
+            for(octave_idx_type k = 0; k < l; k++)
+              {
+                const TSRC* clamp_kb = v - kb * l;
+                const TSRC* clamp_kf = v + (kf + 1) * l;
+                if (clamp_kb < dim_begin + k) clamp_kb = dim_begin + k;
+                if (clamp_kf > dim_end + k) clamp_kf = dim_end + k;
+
+                *r = *clamp_kb;
+                for (const TSRC* h = clamp_kb + l; h < clamp_kf; h += l)
+                    *r *= *h;
+                ++r;
+                ++v;
+              }
+          }
+      }
+  }
+
+// movmax
+template <typename TSRC, typename TRES>
+  inline void
+  mx_inline_movmax (const TSRC *v, TRES *r,
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,
+     octave_idx_type kb, octave_idx_type kf)
+  {
+    for (octave_idx_type i = 0; i < u; i++)
+      {
+        const TSRC* dim_begin = v;
+        const TSRC* dim_end = v + n * l;
+        for(octave_idx_type j = 0; j < n; j++)
+          {
+            for(octave_idx_type k = 0; k < l; k++)
+              {
+                const TSRC* clamp_kb = v - kb * l;
+                const TSRC* clamp_kf = v + (kf + 1) * l;
+                if (clamp_kb < dim_begin + k) clamp_kb = dim_begin + k;
+                if (clamp_kf > dim_end + k) clamp_kf = dim_end + k;
+
+                *r = *clamp_kb;
+                for (const TSRC* h = clamp_kb + l; h < clamp_kf; h += l)
+                    if (*h > *r) *r = *h;
+                ++r;
+                ++v;
+              }
+          }
+      }
+  }
+
+  // movmin
+template <typename TSRC, typename TRES>
+  inline void
+  mx_inline_movmin (const TSRC *v, TRES *r,
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,
+     octave_idx_type kb, octave_idx_type kf)
+  {
+    for (octave_idx_type i = 0; i < u; i++)
+      {
+        const TSRC* dim_begin = v;
+        const TSRC* dim_end = v + n * l;
+        for(octave_idx_type j = 0; j < n; j++)
+          {
+            for(octave_idx_type k = 0; k < l; k++)
+              {
+                const TSRC* clamp_kb = v - kb * l;
+                const TSRC* clamp_kf = v + (kf + 1) * l;
+                if (clamp_kb < dim_begin + k) clamp_kb = dim_begin + k;
+                if (clamp_kf > dim_end + k) clamp_kf = dim_end + k;
+
+                *r = *clamp_kb;
+                for (const TSRC* h = clamp_kb + l; h < clamp_kf; h += l)
+                    if (*h < *r) *r = *h;
+                ++r;
+                ++v;
+              }
+          }
+      }
+  }
+
 // Assistant function
 
 inline void
@@ -1684,6 +1804,25 @@
   return ret;
 }
 
+template <typename R, typename T>
+inline Array<R>
+do_mx_mov_op (const Array<T>& src, int dim,
+              octave_idx_type kb, octave_idx_type kf,
+              void (*mx_mov_op) (const T *, R *,
+                                 octave_idx_type, octave_idx_type, octave_idx_type,
+                                 octave_idx_type, octave_idx_type))
+{
+  octave_idx_type l, n, u;
+  dim_vector dims = src.dims ();
+  get_extent_triplet (dims, dim, l, n, u);
+
+  // Moving operation doesn't reduce the array size.
+  Array<R> ret (dims);
+  mx_mov_op (src.data (), ret.fortran_vec (), l, n, u, kb, kf);
+
+  return ret;
+}
+
 // Fast extra-precise summation.  According to
 // T. Ogita, S. M. Rump, S. Oishi:
 // Accurate Sum And Dot Product,
