changeset 9790:a5035bc7fbfb

rewrite dispatch part & slightly improve min,max,cummin,cummax
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 09 Nov 2009 13:09:49 +0100
parents 97f5de91427b
children 6e425f6be618
files liboctave/CNDArray.h liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/dNDArray.h liboctave/dSparse.cc liboctave/dSparse.h liboctave/fCNDArray.h liboctave/fNDArray.h liboctave/intNDArray.h src/ChangeLog src/DLD-FUNCTIONS/max.cc src/ov-base.cc src/ov-base.h
diffstat 14 files changed, 468 insertions(+), 788 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CNDArray.h
+++ b/liboctave/CNDArray.h
@@ -87,17 +87,17 @@
   ComplexNDArray concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx);
   ComplexNDArray concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx);
 
-  ComplexNDArray max (int dim = 0) const;
-  ComplexNDArray max (Array<octave_idx_type>& index, int dim = 0) const;
-  ComplexNDArray min (int dim = 0) const;
-  ComplexNDArray min (Array<octave_idx_type>& index, int dim = 0) const;
+  ComplexNDArray max (int dim = -1) const;
+  ComplexNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
+  ComplexNDArray min (int dim = -1) const;
+  ComplexNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
 
-  ComplexNDArray cummax (int dim = 0) const;
-  ComplexNDArray cummax (Array<octave_idx_type>& index, int dim = 0) const;
-  ComplexNDArray cummin (int dim = 0) const;
-  ComplexNDArray cummin (Array<octave_idx_type>& index, int dim = 0) const;
+  ComplexNDArray cummax (int dim = -1) const;
+  ComplexNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  ComplexNDArray cummin (int dim = -1) const;
+  ComplexNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
 
-  ComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+  ComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
   ComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   ComplexNDArray& insert (const ComplexNDArray& a, octave_idx_type r, octave_idx_type c);
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -243,25 +243,28 @@
 SparseComplexMatrix
 SparseComplexMatrix::max (int dim) const
 {
-  Array2<octave_idx_type> dummy_idx;
+  Array<octave_idx_type> dummy_idx;
   return max (dummy_idx, dim);
 }
 
 SparseComplexMatrix
-SparseComplexMatrix::max (Array2<octave_idx_type>& idx_arg, int dim) const
+SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
 {
   SparseComplexMatrix result;
   dim_vector dv = dims ();
 
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+  if (dv.numel () == 0 || dim >= dv.length ())
     return result;
  
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
   octave_idx_type nr = dv(0);
   octave_idx_type nc = dv(1);
 
   if (dim == 0)
     {
-      idx_arg.resize (1, nc);
+      idx_arg.clear (1, nc);
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
 	{
@@ -321,7 +324,7 @@
     }
   else
     {
-      idx_arg.resize (nr, 1, 0);
+      idx_arg.resize_fill (nr, 1, 0);
 
       for (octave_idx_type i = cidx(0); i < cidx(1); i++)
 	idx_arg.elem(ridx(i)) = -1;
@@ -395,25 +398,28 @@
 SparseComplexMatrix
 SparseComplexMatrix::min (int dim) const
 {
-  Array2<octave_idx_type> dummy_idx;
+  Array<octave_idx_type> dummy_idx;
   return min (dummy_idx, dim);
 }
 
 SparseComplexMatrix
-SparseComplexMatrix::min (Array2<octave_idx_type>& idx_arg, int dim) const
+SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
 {
   SparseComplexMatrix result;
   dim_vector dv = dims ();
 
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+  if (dv.numel () == 0 || dim >= dv.length ())
     return result;
  
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
   octave_idx_type nr = dv(0);
   octave_idx_type nc = dv(1);
 
   if (dim == 0)
     {
-      idx_arg.resize (1, nc);
+      idx_arg.clear (1, nc);
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
 	{
@@ -473,7 +479,7 @@
     }
   else
     {
-      idx_arg.resize (nr, 1, 0);
+      idx_arg.resize_fill (nr, 1, 0);
 
       for (octave_idx_type i = cidx(0); i < cidx(1); i++)
 	idx_arg.elem(ridx(i)) = -1;
--- a/liboctave/CSparse.h
+++ b/liboctave/CSparse.h
@@ -109,10 +109,10 @@
 
   bool is_hermitian (void) const;
 
-  SparseComplexMatrix max (int dim = 0) const;
-  SparseComplexMatrix max (Array2<octave_idx_type>& index, int dim = 0) const;
-  SparseComplexMatrix min (int dim = 0) const;
-  SparseComplexMatrix min (Array2<octave_idx_type>& index, int dim = 0) const;
+  SparseComplexMatrix max (int dim = -1) const;
+  SparseComplexMatrix max (Array<octave_idx_type>& index, int dim = -1) const;
+  SparseComplexMatrix min (int dim = -1) const;
+  SparseComplexMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
 
   SparseComplexMatrix& insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c);
   SparseComplexMatrix& insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c);
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -4,6 +4,18 @@
 	* dSparse.cc: Update.
 	* CSparse.h (Sparse::max): Use Array<octave_idx_type>.
 	* CSparse.cc: Update.
+	* dNDArray.h (NDArray::max, NDArray::min, NDArray::cummax,
+	NDArray::cummin, NDArray::diff): Use dim = -1 as default.
+	* fNDArray.h (FloatNDArray::max, FloatNDArray::min, FloatNDArray::cummax,
+	FloatNDArray::cummin, FloatNDArray::diff): Use dim = -1 as default.
+	* CNDArray.h (ComplexNDArray::max, ComplexNDArray::min,
+	ComplexNDArray::cummax, ComplexNDArray::cummin, ComplexNDArray::diff):
+	Use dim = -1 as default.
+	* fCNDArray.h (FloatComplexNDArray::max, FloatComplexNDArray::min,
+	FloatComplexNDArray::cummax, FloatComplexNDArray::cummin,
+	FloatComplexNDArray::diff): Use dim = -1 as default.
+	* intNDArray.h (intNDArray::max, intNDArray::min, intNDArray::cummax,
+	intNDArray::cummin, intNDArray::diff): Use dim = -1 as default.
 
 2009-11-03  Jaroslav Hajek  <highegg@gmail.com>
 
--- a/liboctave/dNDArray.h
+++ b/liboctave/dNDArray.h
@@ -99,17 +99,17 @@
   ComplexNDArray concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx);
   charNDArray concat (const charNDArray& rb, const Array<octave_idx_type>& ra_idx);
 
-  NDArray max (int dim = 0) const;
-  NDArray max (Array<octave_idx_type>& index, int dim = 0) const;
-  NDArray min (int dim = 0) const;
-  NDArray min (Array<octave_idx_type>& index, int dim = 0) const;
+  NDArray max (int dim = -1) const;
+  NDArray max (Array<octave_idx_type>& index, int dim = -1) const;
+  NDArray min (int dim = -1) const;
+  NDArray min (Array<octave_idx_type>& index, int dim = -1) const;
   
-  NDArray cummax (int dim = 0) const;
-  NDArray cummax (Array<octave_idx_type>& index, int dim = 0) const;
-  NDArray cummin (int dim = 0) const;
-  NDArray cummin (Array<octave_idx_type>& index, int dim = 0) const;
+  NDArray cummax (int dim = -1) const;
+  NDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  NDArray cummin (int dim = -1) const;
+  NDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
 
-  NDArray diff (octave_idx_type order = 1, int dim = 0) const;
+  NDArray diff (octave_idx_type order = 1, 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);
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -274,25 +274,28 @@
 SparseMatrix
 SparseMatrix::max (int dim) const
 {
-  Array2<octave_idx_type> dummy_idx;
+  Array<octave_idx_type> dummy_idx;
   return max (dummy_idx, dim);
 }
 
 SparseMatrix
-SparseMatrix::max (Array2<octave_idx_type>& idx_arg, int dim) const
+SparseMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
 {
   SparseMatrix result;
   dim_vector dv = dims ();
 
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+  if (dv.numel () == 0 || dim >= dv.length ())
     return result;
  
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
   octave_idx_type nr = dv(0);
   octave_idx_type nc = dv(1);
 
   if (dim == 0)
     {
-      idx_arg.resize (1, nc);
+      idx_arg.clear (1, nc);
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
 	{
@@ -346,7 +349,7 @@
     }
   else
     {
-      idx_arg.resize (nr, 1, 0);
+      idx_arg.resize_fill (nr, 1, 0);
 
       for (octave_idx_type i = cidx(0); i < cidx(1); i++)
 	idx_arg.elem(ridx(i)) = -1;
@@ -420,25 +423,28 @@
 SparseMatrix
 SparseMatrix::min (int dim) const
 {
-  Array2<octave_idx_type> dummy_idx;
+  Array<octave_idx_type> dummy_idx;
   return min (dummy_idx, dim);
 }
 
 SparseMatrix
-SparseMatrix::min (Array2<octave_idx_type>& idx_arg, int dim) const
+SparseMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
 {
   SparseMatrix result;
   dim_vector dv = dims ();
 
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+  if (dv.numel () == 0 || dim >= dv.length ())
     return result;
  
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
   octave_idx_type nr = dv(0);
   octave_idx_type nc = dv(1);
 
   if (dim == 0)
     {
-      idx_arg.resize (1, nc);
+      idx_arg.clear (1, nc);
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
 	{
@@ -492,7 +498,7 @@
     }
   else
     {
-      idx_arg.resize (nr, 1, 0);
+      idx_arg.resize_fill (nr, 1, 0);
 
       for (octave_idx_type i = cidx(0); i < cidx(1); i++)
 	idx_arg.elem(ridx(i)) = -1;
--- a/liboctave/dSparse.h
+++ b/liboctave/dSparse.h
@@ -99,10 +99,10 @@
 
   bool is_symmetric (void) const;
 
-  SparseMatrix max (int dim = 0) const;
-  SparseMatrix max (Array2<octave_idx_type>& index, int dim = 0) const;
-  SparseMatrix min (int dim = 0) const;
-  SparseMatrix min (Array2<octave_idx_type>& index, int dim = 0) const;
+  SparseMatrix max (int dim = -1) const;
+  SparseMatrix max (Array<octave_idx_type>& index, int dim = -1) const;
+  SparseMatrix min (int dim = -1) const;
+  SparseMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
   
   // destructive insert/delete/reorder operations
 
--- a/liboctave/fCNDArray.h
+++ b/liboctave/fCNDArray.h
@@ -87,17 +87,17 @@
   FloatComplexNDArray concat (const FloatComplexNDArray& rb, const Array<octave_idx_type>& ra_idx);
   FloatComplexNDArray concat (const FloatNDArray& rb, const Array<octave_idx_type>& ra_idx);
 
-  FloatComplexNDArray max (int dim = 0) const;
-  FloatComplexNDArray max (Array<octave_idx_type>& index, int dim = 0) const;
-  FloatComplexNDArray min (int dim = 0) const;
-  FloatComplexNDArray min (Array<octave_idx_type>& index, int dim = 0) const;
+  FloatComplexNDArray max (int dim = -1) const;
+  FloatComplexNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatComplexNDArray min (int dim = -1) const;
+  FloatComplexNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
 
-  FloatComplexNDArray cummax (int dim = 0) const;
-  FloatComplexNDArray cummax (Array<octave_idx_type>& index, int dim = 0) const;
-  FloatComplexNDArray cummin (int dim = 0) const;
-  FloatComplexNDArray cummin (Array<octave_idx_type>& index, int dim = 0) const;
+  FloatComplexNDArray cummax (int dim = -1) const;
+  FloatComplexNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatComplexNDArray cummin (int dim = -1) const;
+  FloatComplexNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
 
-  FloatComplexNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+  FloatComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
   FloatComplexNDArray& insert (const NDArray& a, octave_idx_type r, octave_idx_type c);
   FloatComplexNDArray& insert (const FloatComplexNDArray& a, octave_idx_type r, octave_idx_type c);
--- a/liboctave/fNDArray.h
+++ b/liboctave/fNDArray.h
@@ -96,17 +96,17 @@
   FloatComplexNDArray concat (const FloatComplexNDArray& rb, const Array<octave_idx_type>& ra_idx);
   charNDArray concat (const charNDArray& rb, const Array<octave_idx_type>& ra_idx);
 
-  FloatNDArray max (int dim = 0) const;
-  FloatNDArray max (Array<octave_idx_type>& index, int dim = 0) const;
-  FloatNDArray min (int dim = 0) const;
-  FloatNDArray min (Array<octave_idx_type>& index, int dim = 0) const;
+  FloatNDArray max (int dim = -1) const;
+  FloatNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatNDArray min (int dim = -1) const;
+  FloatNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
   
-  FloatNDArray cummax (int dim = 0) const;
-  FloatNDArray cummax (Array<octave_idx_type>& index, int dim = 0) const;
-  FloatNDArray cummin (int dim = 0) const;
-  FloatNDArray cummin (Array<octave_idx_type>& index, int dim = 0) const;
+  FloatNDArray cummax (int dim = -1) const;
+  FloatNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatNDArray cummin (int dim = -1) const;
+  FloatNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
 
-  FloatNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+  FloatNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
   FloatNDArray& insert (const FloatNDArray& a, octave_idx_type r, octave_idx_type c);
   FloatNDArray& insert (const FloatNDArray& a, const Array<octave_idx_type>& ra_idx);
--- a/liboctave/intNDArray.h
+++ b/liboctave/intNDArray.h
@@ -77,21 +77,21 @@
   boolNDArray all (int dim = -1) const;
   boolNDArray any (int dim = -1) const;
 
-  intNDArray max (int dim = 0) const;
-  intNDArray max (Array<octave_idx_type>& index, int dim = 0) const;
-  intNDArray min (int dim = 0) const;
-  intNDArray min (Array<octave_idx_type>& index, int dim = 0) const;
+  intNDArray max (int dim = -1) const;
+  intNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
+  intNDArray min (int dim = -1) const;
+  intNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
   
-  intNDArray cummax (int dim = 0) const;
-  intNDArray cummax (Array<octave_idx_type>& index, int dim = 0) const;
-  intNDArray cummin (int dim = 0) const;
-  intNDArray cummin (Array<octave_idx_type>& index, int dim = 0) const;
+  intNDArray cummax (int dim = -1) const;
+  intNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  intNDArray cummin (int dim = -1) const;
+  intNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
   
   intNDArray sum (int dim) const;
   NDArray dsum (int dim) const;
   intNDArray cumsum (int dim) const;
 
-  intNDArray diff (octave_idx_type order = 1, int dim = 0) const;
+  intNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
   intNDArray abs (void) const;
   intNDArray signum (void) const;
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,16 @@
+2009-11-09  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-base.cc (btyp_mixed_numeric): New function.
+	* ov-base.h: Declare it.
+	(btyp_isnumeric): New inline function.
+	* DLD_FUNCTIONS/max.cc (do_minmax_red_op, do_minmax_bin_op,
+	do_minmax_body, do_cumminmax_red_op): New functions.
+	(MINMAX_DOUBLE_SBODY, MINMAX_DOUBLE_BODY,
+	MINMAX_SINGLE_SBODY, MINMAX_SINGLE_BODY,
+	MINMAX_SPARSE_BODY, MINMAX_INT_SBODY, MINMAX_INT_BODY,
+	MINMAX_BODY, CUMMINMAX_BODY): Remove.
+	(Fmin, Fmax, Fcummin, Fcummax): Update.
+
 2009-11-09  Jaroslav Hajek  <highegg@gmail.com>
 
 	* ov.h: Add sparse matrix extractors.
--- a/src/DLD-FUNCTIONS/max.cc
+++ b/src/DLD-FUNCTIONS/max.cc
@@ -42,646 +42,236 @@
 #include "ov-re-sparse.h"
 #include "ov-cx-sparse.h"
 
-#define MINMAX_DOUBLE_SBODY(FCN) \
-{ \
-  if (nargout == 1 || nargout == 0) \
-    { \
-      if (arg1.is_real_type ()) \
-	{ \
-	  NDArray m = arg1.array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      NDArray n = m. FCN (dim); \
-	      retval(0) = n; \
-	    } \
-	} \
-      else if (arg1.is_complex_type ()) \
-	{ \
-	  ComplexNDArray m = arg1.complex_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      ComplexNDArray n = m. FCN (dim); \
-	      retval(0) = n; \
-	    } \
-	} \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
-    } \
-  else if (nargout == 2) \
-    { \
-      Array<octave_idx_type> index; \
- \
-      if (arg1.is_real_type ()) \
-	{ \
-	  NDArray m = arg1.array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      NDArray n = m. FCN (index, dim);	\
-	      retval(0) = n; \
-	    } \
-	} \
-      else if (arg1.is_complex_type ()) \
-	{ \
-	  ComplexNDArray m = arg1.complex_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      ComplexNDArray n = m. FCN (index, dim);	\
-	      retval(0) = n; \
-	    } \
-	} \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
- \
-      octave_idx_type len = index.numel (); \
- \
-      if (len > 0) \
-	retval(1) = octave_value (index, true, true);	\
-      else \
-	retval(1) = NDArray (); \
-    } \
+template <class ArrayType>
+static octave_value_list
+do_minmax_red_op (const octave_value& arg,
+                  int nargout, int dim, bool ismin)
+{
+  octave_value_list retval;
+  ArrayType array = octave_value_extract<ArrayType> (arg);
+
+  if (error_state)
+    return retval;
+
+  if (nargout == 2)
+    {
+      retval.resize (2);
+      Array<octave_idx_type> idx;
+      if (ismin)
+        retval(0) = array.min (idx, dim);
+      else
+        retval(0) = array.max (idx, dim);
+
+      retval(1) = octave_value (idx, true, true);
+    }
+  else
+    {
+      if (ismin)
+        retval(0) = array.min (dim);
+      else
+        retval(0) = array.max (dim);
+    }
+
+  return retval;
 }
 
-#define MINMAX_DOUBLE_BODY(FCN) \
-{ \
-  bool single_arg = (nargin == 1) || (arg2.is_empty() && nargin == 3);	\
-  if (single_arg) \
-    MINMAX_DOUBLE_SBODY (FCN) \
-  else \
-    { \
-      int arg1_is_scalar = arg1.is_scalar_type (); \
-      int arg2_is_scalar = arg2.is_scalar_type (); \
- \
-      int arg1_is_complex = arg1.is_complex_type (); \
-      int arg2_is_complex = arg2.is_complex_type (); \
- \
-      if (arg1_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      Complex c1 = arg1.complex_value (); \
-	      ComplexNDArray m2 = arg2.complex_array_value (); \
-	      if (! error_state) \
-		{ \
-		  ComplexNDArray result = FCN (c1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      double d1 = arg1.double_value (); \
-	      NDArray m2 = arg2.array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  NDArray result = FCN (d1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else if (arg2_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      ComplexNDArray m1 = arg1.complex_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  Complex c2 = arg2.complex_value (); \
-		  ComplexNDArray result = FCN (m1, c2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      NDArray m1 = arg1.array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  double d2 = arg2.double_value (); \
-		  NDArray result = FCN (m1, d2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      ComplexNDArray m1 = arg1.complex_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  ComplexNDArray m2 = arg2.complex_array_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      ComplexNDArray result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	  else \
-	    { \
-	      NDArray m1 = arg1.array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  NDArray m2 = arg2.array_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      NDArray result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	} \
-    } \
-}
+template <class ArrayType>
+static octave_value
+do_minmax_bin_op (const octave_value& argx, const octave_value& argy,
+                  bool ismin)
+{
+  typedef typename ArrayType::element_type ScalarType;
+
+  octave_value retval;
+
+  if (argx.is_scalar_type () == 1)
+    {
+      ScalarType x = octave_value_extract<ScalarType> (argx);
+      ArrayType y = octave_value_extract<ArrayType> (argy);
 
-#define MINMAX_SINGLE_SBODY(FCN) \
-{ \
-  if (nargout == 1 || nargout == 0) \
-    { \
-      if (arg1.is_real_type ()) \
-	{ \
-	  FloatNDArray m = arg1.float_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      FloatNDArray n = m. FCN (dim); \
-	      retval(0) = n; \
-	    } \
-	} \
-      else if (arg1.is_complex_type ()) \
-	{ \
-	  FloatComplexNDArray m = arg1.float_complex_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      FloatComplexNDArray n = m. FCN (dim); \
-	      retval(0) = n; \
-	    } \
-	} \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
-    } \
-  else if (nargout == 2) \
-    { \
-      Array<octave_idx_type> index; \
- \
-      if (arg1.is_real_type ()) \
-	{ \
-	  FloatNDArray m = arg1.float_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      FloatNDArray n = m. FCN (index, dim);	\
-	      retval(0) = n; \
-	    } \
-	} \
-      else if (arg1.is_complex_type ()) \
-	{ \
-	  FloatComplexNDArray m = arg1.float_complex_array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      FloatComplexNDArray n = m. FCN (index, dim);	\
-	      retval(0) = n; \
-	    } \
-	} \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
- \
-      octave_idx_type len = index.numel (); \
- \
-      if (len > 0) \
-	retval(1) = octave_value (index, true, true);	\
-      else \
-	retval(1) = NDArray (); \
-    } \
+      if (error_state)
+        ;
+      else if (ismin)
+        retval = min (x, y);
+      else
+        retval = max (x, y);
+    }
+  else if (argy.is_scalar_type () == 1)
+    {
+      ArrayType x = octave_value_extract<ArrayType> (argx);
+      ScalarType y = octave_value_extract<ScalarType> (argy);
+
+      if (error_state)
+        ;
+      else if (ismin)
+        retval = min (x, y);
+      else
+        retval = max (x, y);
+    }
+  else
+    {
+      ArrayType x = octave_value_extract<ArrayType> (argx);
+      ArrayType y = octave_value_extract<ArrayType> (argy);
+
+      if (error_state)
+        ;
+      else if (ismin)
+        retval = min (x, y);
+      else
+        retval = max (x, y);
+    }
+
+  return retval;
 }
 
-#define MINMAX_SINGLE_BODY(FCN) \
-{ \
-  bool single_arg = (nargin == 1) || (arg2.is_empty() && nargin == 3);	\
-  if (single_arg) \
-    MINMAX_SINGLE_SBODY(FCN) \
-  else \
-    { \
-      int arg1_is_scalar = arg1.is_scalar_type (); \
-      int arg2_is_scalar = arg2.is_scalar_type (); \
- \
-      int arg1_is_complex = arg1.is_complex_type (); \
-      int arg2_is_complex = arg2.is_complex_type (); \
- \
-      if (arg1_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      FloatComplex c1 = arg1.float_complex_value (); \
-	      FloatComplexNDArray m2 = arg2.float_complex_array_value (); \
-	      if (! error_state) \
-		{ \
-		  FloatComplexNDArray result = FCN (c1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      float d1 = arg1.float_value (); \
-	      FloatNDArray m2 = arg2.float_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  FloatNDArray result = FCN (d1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else if (arg2_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      FloatComplexNDArray m1 = arg1.float_complex_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  FloatComplex c2 = arg2.float_complex_value (); \
-		  FloatComplexNDArray result = FCN (m1, c2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      FloatNDArray m1 = arg1.float_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  float d2 = arg2.float_value (); \
-		  FloatNDArray result = FCN (m1, d2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      FloatComplexNDArray m1 = arg1.float_complex_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  FloatComplexNDArray m2 = arg2.float_complex_array_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      FloatComplexNDArray result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	  else \
-	    { \
-	      FloatNDArray m1 = arg1.float_array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  FloatNDArray m2 = arg2.float_array_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      FloatNDArray result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	} \
-    } \
-}
+static octave_value_list
+do_minmax_body (const octave_value_list& args,
+                int nargout, bool ismin)
+{
+  octave_value_list retval;
+
+  const char *func = ismin ? "min" : "max";
+
+  int nargin = args.length ();
 
-#define MINMAX_INT_SBODY(FCN, TYP) \
-{ \
-  if (nargout == 1 || nargout == 0) \
-    { \
-      TYP ## NDArray m = arg1. TYP ## _array_value (); \
- \
-      if (! error_state) \
-	{ \
-	  TYP ## NDArray n = m. FCN (dim); \
-	  retval(0) = n; \
-	} \
-    } \
-  else if (nargout == 2) \
-    { \
-      Array<octave_idx_type> index; \
- \
-      TYP ## NDArray m = arg1. TYP ## _array_value (); \
- \
-      if (! error_state) \
-        { \
-	  TYP ## NDArray n = m. FCN (index, dim);	\
-	  retval(0) = n; \
-	} \
- \
-      octave_idx_type len = index.numel (); \
- \
-      if (len > 0) \
-	retval(1) = octave_value (index, true, true);	\
-      else \
-	retval(1) = NDArray (); \
-    } \
-}
+  if (nargin == 3 || nargin == 1)
+    {
+      octave_value arg = args(0);
+      int dim = -1;
+      if (nargin == 3)
+        {
+          dim = args(2).int_value (true) - 1;
+          if (error_state || dim < 0)
+            {
+              error ("%s: invalid dimension", func);
+              return retval;
+            }
 
-#define MINMAX_INT_BODY(FCN, TYP) \
- { \
-  bool single_arg = (nargin == 1) || (arg2.is_empty() && nargin == 3);	\
-  if (single_arg) \
-    MINMAX_INT_SBODY (FCN, TYP) \
-  else \
-    { \
-      int arg1_is_scalar = arg1.is_scalar_type (); \
-      int arg2_is_scalar = arg2.is_scalar_type (); \
- \
-      if (arg1_is_scalar) \
-	{ \
-	  octave_ ## TYP d1 = arg1. TYP ## _scalar_value (); \
-	  TYP ## NDArray m2 = arg2. TYP ## _array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      TYP ## NDArray result = FCN (d1, m2); \
-	      if (! error_state) \
-		retval(0) = result; \
-	    } \
-	} \
-      else if (arg2_is_scalar) \
-	{ \
-	  TYP ## NDArray m1 = arg1. TYP ## _array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      octave_ ## TYP d2 = arg2. TYP ## _scalar_value (); \
-	      TYP ## NDArray result = FCN (m1, d2); \
-	      if (! error_state) \
-		retval(0) = result; \
-	    } \
-	} \
-      else \
-	{ \
-	  TYP ## NDArray m1 = arg1. TYP ## _array_value (); \
- \
-	  if (! error_state) \
-	    { \
-	      TYP ## NDArray m2 = arg2. TYP ## _array_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  TYP ## NDArray result = FCN (m1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-    } \
-}
+          if (! args(1).is_empty ())
+            warning ("%s: second argument is ignored");
+        }
 
-#define MINMAX_SPARSE_BODY(FCN) \
-{ \
-  bool single_arg = (nargin == 1) || arg2.is_empty();	\
- \
-  if (single_arg && (nargout == 1 || nargout == 0)) \
-    { \
-      if (arg1.is_real_type ()) \
-	retval(0) = arg1.sparse_matrix_value () .FCN (dim); \
-      else if (arg1.is_complex_type ()) \
-	retval(0) = arg1.sparse_complex_matrix_value () .FCN (dim); \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
-    } \
-  else if (single_arg && nargout == 2) \
-    { \
-      Array2<octave_idx_type> index; \
- \
-      if (arg1.is_real_type ()) \
-	retval(0) = arg1.sparse_matrix_value () .FCN (index, dim); \
-      else if (arg1.is_complex_type ()) \
-	retval(0) = arg1.sparse_complex_matrix_value () .FCN (index, dim); \
-      else \
-	gripe_wrong_type_arg (#FCN, arg1); \
- \
-      octave_idx_type len = index.numel (); \
- \
-      if (len > 0) \
-	retval(1) = octave_value (index, true, true);	\
-      else \
-	retval(1) = NDArray (); \
-    } \
-  else \
-    { \
-      int arg1_is_scalar = arg1.is_scalar_type (); \
-      int arg2_is_scalar = arg2.is_scalar_type (); \
- \
-      int arg1_is_complex = arg1.is_complex_type (); \
-      int arg2_is_complex = arg2.is_complex_type (); \
- \
-      if (arg1_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      Complex c1 = arg1.complex_value (); \
-	      \
-	      SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \
-	      \
-	      if (! error_state) \
-		{ \
-		  SparseComplexMatrix result = FCN (c1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      double d1 = arg1.double_value (); \
-	      SparseMatrix m2 = arg2.sparse_matrix_value (); \
-	      \
-	      if (! error_state) \
-		{ \
-		  SparseMatrix result = FCN (d1, m2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else if (arg2_is_scalar) \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  Complex c2 = arg2.complex_value (); \
-		  SparseComplexMatrix result = FCN (m1, c2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	  else \
-	    { \
-	      SparseMatrix m1 = arg1.sparse_matrix_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  double d2 = arg2.double_value (); \
-		  SparseMatrix result = FCN (m1, d2); \
-		  if (! error_state) \
-		    retval(0) = result; \
-		} \
-	    } \
-	} \
-      else \
-	{ \
-	  if (arg1_is_complex || arg2_is_complex) \
-	    { \
-	      SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      SparseComplexMatrix result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	  else \
-	    { \
-	      SparseMatrix m1 = arg1.sparse_matrix_value (); \
- \
-	      if (! error_state) \
-		{ \
-		  SparseMatrix m2 = arg2.sparse_matrix_value (); \
- \
-		  if (! error_state) \
-		    { \
-		      SparseMatrix result = FCN (m1, m2); \
-		      if (! error_state) \
-			retval(0) = result; \
-		    } \
-		} \
-	    } \
-	} \
-    } \
+      switch (arg.builtin_type ())
+        {
+        case btyp_double:
+          {
+            if (arg.is_range () && (dim == -1 || dim == 1))
+              {
+                Range range = arg.range_value ();
+                if (range.nelem () == 0)
+                  {
+                    retval(0) = arg;
+                    if (nargout > 1)
+                      retval(1) = arg;
+                  }
+                else if (ismin)
+                  {
+                    retval(0) = range.min ();
+                    if (nargout > 1)
+                      retval(1) = static_cast<double> (range.inc () < 0 ? range.nelem () : 1);
+                  }
+                else
+                  {
+                    retval(0) = range.max ();
+                    if (nargout > 1)
+                      retval(1) = static_cast<double> (range.inc () >= 0 ? range.nelem () : 1);
+                  }
+              }
+            else if (arg.is_sparse_type ())
+              retval = do_minmax_red_op<SparseMatrix> (arg, nargout, dim, ismin);
+            else
+              retval = do_minmax_red_op<NDArray> (arg, nargout, dim, ismin);
+            break;
+          }
+        case btyp_complex:
+          {
+            if (arg.is_sparse_type ())
+              retval = do_minmax_red_op<SparseComplexMatrix> (arg, nargout, dim, ismin);
+            else
+              retval = do_minmax_red_op<ComplexNDArray> (arg, nargout, dim, ismin);
+            break;
+          }
+        case btyp_float:
+          retval = do_minmax_red_op<FloatNDArray> (arg, nargout, dim, ismin);
+          break;
+        case btyp_float_complex:
+          retval = do_minmax_red_op<FloatComplexNDArray> (arg, nargout, dim, ismin);
+          break;
+#define MAKE_INT_BRANCH(X) \
+        case btyp_ ## X: \
+          retval = do_minmax_red_op<X ## NDArray> (arg, nargout, dim, ismin); \
+          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:
+          gripe_wrong_type_arg (func, arg);
+      }
+    }
+  else if (nargin == 2)
+    {
+      octave_value argx = args(0), argy = args(1);
+      builtin_type_t xtyp = argx.builtin_type (), ytyp = argy.builtin_type ();
+      builtin_type_t rtyp = btyp_mixed_numeric (xtyp, ytyp);
+
+      switch (rtyp)
+        {
+        case btyp_double:
+          {
+            if ((argx.is_sparse_type () 
+                 && (argy.is_sparse_type () || argy.is_scalar_type ()))
+                || (argy.is_sparse_type () && argx.is_scalar_type ()))
+              retval = do_minmax_bin_op<SparseMatrix> (argx, argy, ismin);
+            else
+              retval = do_minmax_bin_op<NDArray> (argx, argy, ismin);
+            break;
+          }
+        case btyp_complex:
+          {
+            if ((argx.is_sparse_type () 
+                 && (argy.is_sparse_type () || argy.is_scalar_type ()))
+                || (argy.is_sparse_type () && argx.is_scalar_type ()))
+              retval = do_minmax_bin_op<SparseComplexMatrix> (argx, argy, ismin);
+            else
+              retval = do_minmax_bin_op<ComplexNDArray> (argx, argy, ismin);
+            break;
+          }
+        case btyp_float:
+          retval = do_minmax_bin_op<FloatNDArray> (argx, argy, ismin);
+          break;
+        case btyp_float_complex:
+          retval = do_minmax_bin_op<FloatComplexNDArray> (argx, argy, ismin);
+          break;
+#define MAKE_INT_BRANCH(X) \
+        case btyp_ ## X: \
+          retval = do_minmax_bin_op<X ## NDArray> (argx, argy, ismin); \
+          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:
+          error ("%s: cannot compute %s (%s, %s)", func, func,
+                 argx.type_name ().c_str (), argy.type_name ().c_str ());
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
 }
 
-
-#define MINMAX_BODY(FCN) \
- \
-  octave_value_list retval;  \
- \
-  int nargin = args.length (); \
- \
-  if (nargin < 1 || nargin > 3 || nargout > 2) \
-    { \
-      print_usage (); \
-      return retval; \
-    } \
- \
-  octave_value arg1; \
-  octave_value arg2; \
-  octave_value arg3; \
- \
-  switch (nargin) \
-    { \
-    case 3: \
-      arg3 = args(2); \
- \
-    case 2: \
-      arg2 = args(1); \
- \
-    case 1: \
-      arg1 = args(0); \
-      break; \
- \
-    default: \
-      panic_impossible (); \
-      break; \
-    } \
- \
-  int dim; \
-  dim_vector dv = arg1.dims (); \
-  if (error_state) \
-    { \
-      gripe_wrong_type_arg (#FCN, arg1);  \
-      return retval; \
-    } \
- \
-  if (nargin == 3) \
-    { \
-      dim = arg3.nint_value () - 1;  \
-      if (dim < 0 || dim >= dv.length ()) \
-        { \
-	  error ("%s: invalid dimension", #FCN); \
-	  return retval; \
-	} \
-    } \
-  else \
-    { \
-      dim = 0; \
-      while ((dim < dv.length ()) && (dv (dim) <= 1)) \
-	dim++; \
-      if (dim == dv.length ()) \
-	dim = 0; \
-    } \
- \
-  if (arg1.is_integer_type ()) \
-    { \
-      if (arg1.is_uint8_type ()) \
-        MINMAX_INT_BODY (FCN, uint8) \
-      else if (arg1.is_uint16_type ()) \
-        MINMAX_INT_BODY (FCN, uint16) \
-      else if (arg1.is_uint32_type ()) \
-        MINMAX_INT_BODY (FCN, uint32) \
-      else if (arg1.is_uint64_type ()) \
-        MINMAX_INT_BODY (FCN, uint64) \
-      else if (arg1.is_int8_type ()) \
-        MINMAX_INT_BODY (FCN, int8) \
-      else if (arg1.is_int16_type ()) \
-        MINMAX_INT_BODY (FCN, int16) \
-      else if (arg1.is_int32_type ()) \
-        MINMAX_INT_BODY (FCN, int32) \
-      else if (arg1.is_int64_type ()) \
-        MINMAX_INT_BODY (FCN, int64) \
-    } \
-  else if (arg1.is_sparse_type ()) \
-    MINMAX_SPARSE_BODY (FCN) \
-  else if (arg1.is_single_type ()) \
-    MINMAX_SINGLE_BODY (FCN) \
-  else \
-    MINMAX_DOUBLE_BODY (FCN) \
- \
- return retval;
-
 DEFUN_DLD (min, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn  {Loadable Function} {} min (@var{x})\n\
@@ -728,7 +318,7 @@
 @seealso{max, cummin, cummax}\n\
 @end deftypefn")
 {
-  MINMAX_BODY (min);
+  return do_minmax_body (args, nargout, true);
 }
 
 /*
@@ -804,7 +394,7 @@
 @seealso{min, cummax, cummin}\n\
 @end deftypefn")
 {
-  MINMAX_BODY (max);
+  return do_minmax_body (args, nargout, false);
 }
 
 /* 
@@ -835,86 +425,99 @@
 
 */
 
-#define CUMMINMAX_BODY(FCN) \
- \
-  octave_value_list retval;  \
- \
-  int nargin = args.length (); \
- \
-  if (nargin < 1 || nargin > 2 || nargout > 2) \
-    { \
-      print_usage (); \
-      return retval; \
-    } \
- \
-  octave_value arg1; \
-  octave_value arg2; \
- \
-  switch (nargin) \
-    { \
-    case 2: \
-      arg2 = args(1); \
- \
-    case 1: \
-      arg1 = args(0); \
-      break; \
- \
-    default: \
-      panic_impossible (); \
-      break; \
-    } \
- \
-  int dim; \
-  dim_vector dv = arg1.dims (); \
-  if (error_state) \
-    { \
-      gripe_wrong_type_arg (#FCN, arg1);  \
-      return retval; \
-    } \
- \
-  if (nargin == 2) \
-    { \
-      dim = arg2.nint_value () - 1;  \
-      if (dim < 0 || dim >= dv.length ()) \
-        { \
-	  error ("%s: invalid dimension", #FCN); \
-	  return retval; \
-	} \
-    } \
-  else \
-    { \
-      dim = 0; \
-      while ((dim < dv.length ()) && (dv (dim) <= 1)) \
-	dim++; \
-      if (dim == dv.length ()) \
-	dim = 0; \
-    } \
- \
-  if (arg1.is_integer_type ()) \
-    { \
-      if (arg1.is_uint8_type ()) \
-        MINMAX_INT_SBODY (FCN, uint8) \
-      else if (arg1.is_uint16_type ()) \
-        MINMAX_INT_SBODY (FCN, uint16) \
-      else if (arg1.is_uint32_type ()) \
-        MINMAX_INT_SBODY (FCN, uint32) \
-      else if (arg1.is_uint64_type ()) \
-        MINMAX_INT_SBODY (FCN, uint64) \
-      else if (arg1.is_int8_type ()) \
-        MINMAX_INT_SBODY (FCN, int8) \
-      else if (arg1.is_int16_type ()) \
-        MINMAX_INT_SBODY (FCN, int16) \
-      else if (arg1.is_int32_type ()) \
-        MINMAX_INT_SBODY (FCN, int32) \
-      else if (arg1.is_int64_type ()) \
-        MINMAX_INT_SBODY (FCN, int64) \
-    } \
-  else if (arg1.is_single_type ()) \
-    MINMAX_SINGLE_SBODY (FCN) \
-  else \
-    MINMAX_DOUBLE_SBODY (FCN) \
- \
- return retval;
+template <class ArrayType>
+static octave_value_list
+do_cumminmax_red_op (const octave_value& arg,
+                     int nargout, int dim, bool ismin)
+{
+  octave_value_list retval;
+  ArrayType array = octave_value_extract<ArrayType> (arg);
+
+  if (error_state)
+    return retval;
+
+  if (nargout == 2)
+    {
+      retval.resize (2);
+      Array<octave_idx_type> idx;
+      if (ismin)
+        retval(0) = array.cummin (idx, dim);
+      else
+        retval(0) = array.cummax (idx, dim);
+
+      retval(1) = octave_value (idx, true, true);
+    }
+  else
+    {
+      if (ismin)
+        retval(0) = array.cummin (dim);
+      else
+        retval(0) = array.cummax (dim);
+    }
+
+  return retval;
+}
+
+static octave_value_list
+do_cumminmax_body (const octave_value_list& args,
+                   int nargout, bool ismin)
+{
+  octave_value_list retval;
+
+  const char *func = ismin ? "cummin" : "cummax";
+
+  int nargin = args.length ();
+
+  if (nargin == 1 || nargin == 2)
+    {
+      octave_value arg = args(0);
+      int dim = -1;
+      if (nargin == 2)
+        {
+          dim = args(1).int_value (true) - 1;
+          if (error_state || dim < 0)
+            {
+              error ("%s: invalid dimension", func);
+              return retval;
+            }
+        }
+
+      switch (arg.builtin_type ())
+        {
+        case btyp_double:
+          retval = do_cumminmax_red_op<NDArray> (arg, nargout, dim, ismin);
+          break;
+        case btyp_complex:
+          retval = do_cumminmax_red_op<ComplexNDArray> (arg, nargout, dim, ismin);
+          break;
+        case btyp_float:
+          retval = do_cumminmax_red_op<FloatNDArray> (arg, nargout, dim, ismin);
+          break;
+        case btyp_float_complex:
+          retval = do_cumminmax_red_op<FloatComplexNDArray> (arg, nargout, dim, ismin);
+          break;
+#define MAKE_INT_BRANCH(X) \
+        case btyp_ ## X: \
+          retval = do_cumminmax_red_op<X ## NDArray> (arg, nargout, dim, ismin); \
+          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:
+          gripe_wrong_type_arg (func, arg);
+      }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
 
 DEFUN_DLD (cummin, args, nargout,
   "-*- texinfo -*-\n\
@@ -955,7 +558,7 @@
 @seealso{cummax, min, max}\n\
 @end deftypefn")
 {
-  CUMMINMAX_BODY (cummin);
+  return do_cumminmax_body (args, nargout, true);
 }
 
 DEFUN_DLD (cummax, args, nargout,
@@ -996,7 +599,7 @@
 @seealso{cummin, max, min}\n\
 @end deftypefn")
 {
-  CUMMINMAX_BODY (cummax);
+  return do_cumminmax_body (args, nargout, false);
 }
 
 /*
--- a/src/ov-base.cc
+++ b/src/ov-base.cc
@@ -54,6 +54,30 @@
 #include "utils.h"
 #include "variables.h"
 
+builtin_type_t btyp_mixed_numeric (builtin_type_t x, builtin_type_t y)
+{
+  builtin_type_t retval = btyp_unknown;
+
+  if (x == btyp_bool)
+    x = btyp_double;
+  if (y == btyp_bool)
+    y = btyp_double;
+
+  if (x <= btyp_float_complex && y <= btyp_float_complex)
+    retval = static_cast<builtin_type_t> (x | y);
+  else if (x <= btyp_uint64 && y <= btyp_float)
+    retval = x;
+  else if (x <= btyp_float && y <= btyp_uint64)
+    retval = y;
+  else if ((x >= btyp_int8 && x <= btyp_int64 
+            && y >= btyp_int8 && y <= btyp_int64)
+           || (x >= btyp_uint8 && x <= btyp_uint64
+               && y >= btyp_uint8 && y <= btyp_uint64))
+    retval = (x > y) ? x : y;
+
+  return retval;
+}
+
 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_base_value,
 				     "<unknown type>", "unknown");
 
--- a/src/ov-base.h
+++ b/src/ov-base.h
@@ -69,12 +69,28 @@
   btyp_uint16,
   btyp_uint32,
   btyp_uint64,
+  btyp_bool,
   btyp_char,
-  btyp_bool,
   btyp_unknown,
   btyp_num_types = btyp_unknown
 };
 
+inline bool btyp_isnumeric (builtin_type_t btyp)
+{ return btyp <= btyp_uint64; }
+
+// Compute a numeric type for a possibly mixed-type operation, using these rules:
+// bool -> double
+// single + double -> single
+// real + complex -> complex
+// integer + real -> integer
+// uint + uint -> uint (the bigger one)
+// sint + sint -> sint (the bigger one)
+//
+// failing otherwise.
+
+extern OCTINTERP_API
+builtin_type_t btyp_mixed_numeric (builtin_type_t x, builtin_type_t y);
+
 template <class T>
 struct class_to_btyp
 {