# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1234770720 -3600
# Node ID 9f7ce4bf76508cbe02dab056b8b1cb0bdf65ef00
# Parent  8af4ba6b4216491222a4ed282b9edd779b58e7ed
optimize min/max functions

diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/CNDArray.cc
--- a/liboctave/CNDArray.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/CNDArray.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -692,191 +692,25 @@
 ComplexNDArray
 ComplexNDArray::max (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return do_mx_minmax_op<ComplexNDArray> (*this, dim, mx_inline_max);
 }
 
 ComplexNDArray
 ComplexNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return ComplexNDArray ();
-  
-  dr(dim) = 1;
-
-  ComplexNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      Complex tmp_max;
-
-      double abs_max = octave_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_max = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_max))
-	    {
-	      abs_max = std::abs(tmp_max);
-	      break;
-	    }
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  Complex tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-
-	  double abs_tmp = std::abs (tmp);
-
-	  if (abs_tmp > abs_max)
-	    {
-	      idx_j = j;
-	      tmp_max = tmp;
-	      abs_max = abs_tmp;
-	    }
-	}
-
-      if (xisnan (tmp_max))
-	{
-	  result.elem (i) = Complex_NaN_result;
-	  idx_arg.elem (i) = 0;
-	}
-      else
-	{
-	  result.elem (i) = tmp_max;
-	  idx_arg.elem (i) = idx_j;
-	}
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<ComplexNDArray> (*this, idx_arg, dim, mx_inline_max);
 }
 
 ComplexNDArray
 ComplexNDArray::min (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return do_mx_minmax_op<ComplexNDArray> (*this, dim, mx_inline_min);
 }
 
 ComplexNDArray
 ComplexNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return ComplexNDArray ();
-  
-  dr(dim) = 1;
-
-  ComplexNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      Complex tmp_min;
-
-      double abs_min = octave_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_min = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_min))
-	    {
-	      abs_min = std::abs(tmp_min);
-	      break;
-	    }
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  Complex tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-
-	  double abs_tmp = std::abs (tmp);
-
-	  if (abs_tmp < abs_min)
-	    {
-	      idx_j = j;
-	      tmp_min = tmp;
-	      abs_min = abs_tmp;
-	    }
-	}
-
-      if (xisnan (tmp_min))
-	{
-	  result.elem (i) = Complex_NaN_result;
-	  idx_arg.elem (i) = 0;
-	}
-      else
-	{
-	  result.elem (i) = tmp_min;
-	  idx_arg.elem (i) = idx_j;
-	}
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<ComplexNDArray> (*this, idx_arg, dim, mx_inline_min);
 }
 
 NDArray
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/ChangeLog	Mon Feb 16 08:52:00 2009 +0100
@@ -1,3 +1,16 @@
+2009-02-16  Jaroslav Hajek  <highegg@gmail.com>
+
+	* oct-cmplx.h (operator <, operator >): New operators.
+	* mx-inlines.cc (OP_MINMAX_FCN, OP_MINMAX_FCN2, OP_MINMAX_FCNN):
+	New macros.
+	(mx_inline_min, mx_inline_max, do_mx_minmax_op): New overloaded
+	template functions.
+	* dNDArray (NDArray::min, NDArray::max): Use do_mx_minmax_op.
+	* fNDArray (FloatNDArray::min, FloatNDArray::max): Ditto.
+	* CNDArray (ComplexNDArray::min, ComplexNDArray::max): Ditto.
+	* fCNDArray (FloatComplexNDArray::min, FloatComplexNDArray::max):
+	Ditto.
+
 2009-02-16  Jaroslav Hajek  <highegg@gmail.com>
 
 	* chMatrix.cc (charMatrix::all, charMatrix::any): Use do_mx_red_op.
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/dNDArray.cc
--- a/liboctave/dNDArray.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/dNDArray.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -731,157 +731,25 @@
 NDArray
 NDArray::max (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_max);
 }
 
 NDArray
 NDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return NDArray ();
-  
-  dr(dim) = 1;
-
-  NDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      double tmp_max = octave_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_max = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_max))
-	    break;
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  double tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-	  else if (tmp > tmp_max)
-	    {
-	      idx_j = j;
-	      tmp_max = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_max;
-      idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_max);
 }
 
 NDArray
 NDArray::min (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_min);
 }
 
 NDArray
 NDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return NDArray ();
-  
-  dr(dim) = 1;
-
-  NDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      double tmp_min = octave_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_min = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_min))
-	    break;
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  double tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-	  else if (tmp < tmp_min)
-	    {
-	      idx_j = j;
-	      tmp_min = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_min;
-      idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_min);
 }
 
 NDArray
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/fCNDArray.cc
--- a/liboctave/fCNDArray.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/fCNDArray.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -687,191 +687,25 @@
 FloatComplexNDArray
 FloatComplexNDArray::max (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return do_mx_minmax_op<FloatComplexNDArray> (*this, dim, mx_inline_max);
 }
 
 FloatComplexNDArray
 FloatComplexNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return FloatComplexNDArray ();
-  
-  dr(dim) = 1;
-
-  FloatComplexNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      FloatComplex tmp_max;
-
-      float abs_max = octave_Float_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_max = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_max))
-	    {
-	      abs_max = std::abs(tmp_max);
-	      break;
-	    }
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  FloatComplex tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-
-	  float abs_tmp = std::abs (tmp);
-
-	  if (abs_tmp > abs_max)
-	    {
-	      idx_j = j;
-	      tmp_max = tmp;
-	      abs_max = abs_tmp;
-	    }
-	}
-
-      if (xisnan (tmp_max))
-	{
-	  result.elem (i) = FloatComplex_NaN_result;
-	  idx_arg.elem (i) = 0;
-	}
-      else
-	{
-	  result.elem (i) = tmp_max;
-	  idx_arg.elem (i) = idx_j;
-	}
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<FloatComplexNDArray> (*this, idx_arg, dim, mx_inline_max);
 }
 
 FloatComplexNDArray
 FloatComplexNDArray::min (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return do_mx_minmax_op<FloatComplexNDArray> (*this, dim, mx_inline_min);
 }
 
 FloatComplexNDArray
 FloatComplexNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return FloatComplexNDArray ();
-  
-  dr(dim) = 1;
-
-  FloatComplexNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      FloatComplex tmp_min;
-
-      float abs_min = octave_Float_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_min = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_min))
-	    {
-	      abs_min = std::abs(tmp_min);
-	      break;
-	    }
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  FloatComplex tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-
-	  float abs_tmp = std::abs (tmp);
-
-	  if (abs_tmp < abs_min)
-	    {
-	      idx_j = j;
-	      tmp_min = tmp;
-	      abs_min = abs_tmp;
-	    }
-	}
-
-      if (xisnan (tmp_min))
-	{
-	  result.elem (i) = FloatComplex_NaN_result;
-	  idx_arg.elem (i) = 0;
-	}
-      else
-	{
-	  result.elem (i) = tmp_min;
-	  idx_arg.elem (i) = idx_j;
-	}
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<FloatComplexNDArray> (*this, idx_arg, dim, mx_inline_min);
 }
 
 FloatNDArray
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/fNDArray.cc
--- a/liboctave/fNDArray.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/fNDArray.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -686,157 +686,25 @@
 FloatNDArray
 FloatNDArray::max (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return do_mx_minmax_op<FloatNDArray> (*this, dim, mx_inline_max);
 }
 
 FloatNDArray
 FloatNDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return FloatNDArray ();
-  
-  dr(dim) = 1;
-
-  FloatNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      float tmp_max = octave_Float_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_max = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_max))
-	    break;
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  float tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-	  else if (tmp > tmp_max)
-	    {
-	      idx_j = j;
-	      tmp_max = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_max;
-      idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<FloatNDArray> (*this, idx_arg, dim, mx_inline_max);
 }
 
 FloatNDArray
 FloatNDArray::min (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return do_mx_minmax_op<FloatNDArray> (*this, dim, mx_inline_min);
 }
 
 FloatNDArray
 FloatNDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = dims ();
-  dim_vector dr = dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return FloatNDArray ();
-  
-  dr(dim) = 1;
-
-  FloatNDArray result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j;
-
-      float tmp_min = octave_Float_NaN;
-
-      for (idx_j = 0; idx_j < x_len; idx_j++)
-	{
-	  tmp_min = elem (idx_j * x_stride + x_offset);
-	  
-	  if (! xisnan (tmp_min))
-	    break;
-	}
-
-      for (octave_idx_type j = idx_j+1; j < x_len; j++)
-	{
-	  float tmp = elem (j * x_stride + x_offset);
-
-	  if (xisnan (tmp))
-	    continue;
-	  else if (tmp < tmp_min)
-	    {
-	      idx_j = j;
-	      tmp_min = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_min;
-      idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<FloatNDArray> (*this, idx_arg, dim, mx_inline_min);
 }
 
 FloatNDArray
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/intNDArray.cc
--- a/liboctave/intNDArray.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/intNDArray.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -213,140 +213,28 @@
 intNDArray<T>
 intNDArray<T>::max (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return do_mx_minmax_op<intNDArray<T> > (*this, dim, mx_inline_max);
 }
 
 template <class T>
 intNDArray<T>
 intNDArray<T>::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = this->dims ();
-  dim_vector dr = this->dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return intNDArray<T> ();
-  
-  dr(dim) = 1;
-
-  intNDArray<T> result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j = 0;
-
-      T tmp_max = this->elem (x_offset);
-
-      for (octave_idx_type j = 1; j < x_len; j++)
-	{
-	  T tmp = this->elem (j * x_stride + x_offset);
-
-	  if (tmp > tmp_max)
-	    {
-	      idx_j = j;
-	      tmp_max = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_max;
-      idx_arg.elem (i) = idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_max);
 }
 
 template <class T>
 intNDArray<T>
 intNDArray<T>::min (int dim) const
 {
-  ArrayN<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return do_mx_minmax_op<intNDArray<T> > (*this, dim, mx_inline_min);
 }
 
 template <class T>
 intNDArray<T>
 intNDArray<T>::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
 {
-  dim_vector dv = this->dims ();
-  dim_vector dr = this->dims ();
-
-  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
-    return intNDArray<T> ();
-  
-  dr(dim) = 1;
-
-  intNDArray<T> result (dr);
-  idx_arg.resize (dr);
-
-  octave_idx_type x_stride = 1;
-  octave_idx_type x_len = dv(dim);
-  for (int i = 0; i < dim; i++)
-    x_stride *= dv(i);
-
-  for (octave_idx_type i = 0; i < dr.numel (); i++)
-    {
-      octave_idx_type x_offset;
-      if (x_stride == 1)
-	x_offset = i * x_len;
-      else
-	{
-	  octave_idx_type x_offset2 = 0;
-	  x_offset = i;
-	  while (x_offset >= x_stride)
-	    {
-	      x_offset -= x_stride;
-	      x_offset2++;
-	    }
-	  x_offset += x_offset2 * x_stride * x_len;
-	}
-
-      octave_idx_type idx_j = 0;
-
-      T tmp_min = this->elem (x_offset);
-
-      for (octave_idx_type j = 1; j < x_len; j++)
-	{
-	  T tmp = this->elem (j * x_stride + x_offset);
-
-	  if (tmp < tmp_min)
-	    {
-	      idx_j = j;
-	      tmp_min = tmp;
-	    }
-	}
-
-      result.elem (i) = tmp_min;
-      idx_arg.elem (i) = idx_j;
-    }
-
-  result.chop_trailing_singletons ();
-  idx_arg.chop_trailing_singletons ();
-
-  return result;
+  return do_mx_minmax_op<intNDArray<T> > (*this, idx_arg, dim, mx_inline_min);
 }
 
 /*
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/lo-mappers.cc
--- a/liboctave/lo-mappers.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/lo-mappers.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -222,44 +222,18 @@
 
 // (double, double) -> double mappers.
 
-// FIXME -- need to handle NA too?
+// According to Matlab, is both args are NaN, the first one is returned.
 
 double
 xmin (double x, double y)
 {
-  if (x < y)
-    return x;
-
-  if (y <= x)
-    return y;
-
-  if (xisnan (x) && ! xisnan (y))
-    return y;
-  else if (xisnan (y) && ! xisnan (x))
-    return x;
-  else if (octave_is_NA (x) || octave_is_NA (y))
-    return octave_NA;
-  else
-    return octave_NaN;
+  return  xisnan (y) ? x : (x <= y ? x : y);
 }
 
 double
 xmax (double x, double y)
 {
-  if (x > y)
-    return x;
-
-  if (y >= x)
-    return y;
-
-  if (xisnan (x) && ! xisnan (y))
-    return y;
-  else if (xisnan (y) && ! xisnan (x))
-    return x;
-  else if (octave_is_NA (x) || octave_is_NA (y))
-    return octave_NA;
-  else
-    return octave_NaN;
+  return  xisnan (y) ? x : (x >= y ? x : y);
 }
 
 // complex -> complex mappers.
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/mx-inlines.cc
--- a/liboctave/mx-inlines.cc	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/mx-inlines.cc	Mon Feb 16 08:52:00 2009 +0100
@@ -467,6 +467,163 @@
 OP_CUM_FCNN (mx_inline_cumsum)
 OP_CUM_FCNN (mx_inline_cumprod)
 
+#define OP_MINMAX_FCN(F, OP) \
+template <class T> \
+void F (const T *v, T *r, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  T tmp = v[0]; \
+  octave_idx_type i = 1; \
+  while (xisnan (tmp) && i < n) tmp = v[i++]; \
+  for (i = 1; i < n; i++) \
+    if (v[i] OP tmp) tmp = v[i]; \
+  *r = tmp; \
+} \
+template <class T> \
+void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  T tmp = v[0]; \
+  octave_idx_type tmpi = 0; \
+  octave_idx_type i = 1; \
+  while (xisnan (tmp) && i < n) tmp = v[i++]; \
+  for (i = 1; i < n; i++) \
+    if (v[i] OP tmp) { tmp = v[i]; tmpi = i; }\
+  *r = tmp; \
+  *ri = tmpi; \
+}
+
+OP_MINMAX_FCN (mx_inline_min, <)
+OP_MINMAX_FCN (mx_inline_max, >)
+
+// Row reductions will be slightly complicated.  We will proceed with checks
+// for NaNs until we detect that no row will yield a NaN, in which case we
+// proceed to a faster code.
+
+#define OP_MINMAX_FCN2(F, OP) \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  bool nan = false; \
+  octave_idx_type j = 0; \
+  for (octave_idx_type i = 0; i < m; i++) \
+    {  \
+      r[i] = v[i]; \
+      if (xisnan (v[i])) nan = true;  \
+    } \
+  j++; v += m; \
+  while (nan && j < n) \
+    { \
+      nan = false; \
+      for (octave_idx_type i = 0; i < m; i++) \
+        {  \
+          if (xisnan (v[i])) \
+            nan = true;  \
+          else if (xisnan (r[i]) || v[i] OP r[i]) \
+            r[i] = v[i]; \
+        } \
+      j++; v += m; \
+    } \
+  while (j < n) \
+    { \
+      for (octave_idx_type i = 0; i < m; i++) \
+        if (v[i] OP r[i]) r[i] = v[i]; \
+      j++; v += m; \
+    } \
+} \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type *ri, \
+   octave_idx_type m, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  bool nan = false; \
+  octave_idx_type j = 0; \
+  for (octave_idx_type i = 0; i < m; i++) \
+    {  \
+      r[i] = v[i]; ri[i] = j; \
+      if (xisnan (v[i])) nan = true;  \
+    } \
+  j++; v += m; \
+  while (nan && j < n) \
+    { \
+      nan = false; \
+      for (octave_idx_type i = 0; i < m; i++) \
+        {  \
+          if (xisnan (v[i])) \
+            nan = true;  \
+          else if (xisnan (r[i]) || v[i] OP r[i]) \
+            { r[i] = v[i]; ri[i] = j; } \
+        } \
+      j++; v += m; \
+    } \
+  while (j < n) \
+    { \
+      for (octave_idx_type i = 0; i < m; i++) \
+        if (v[i] OP r[i]) \
+          { r[i] = v[i]; ri[i] = j; } \
+      j++; v += m; \
+    } \
+}
+
+OP_MINMAX_FCN2 (mx_inline_min, <)
+OP_MINMAX_FCN2 (mx_inline_max, >)
+
+#define OP_MINMAX_FCNN(F) \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type l, \
+   octave_idx_type n, octave_idx_type u) \
+{ \
+  if (! n) return; \
+  if (l == 1) \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, n); \
+          v += n; r++; \
+        } \
+    } \
+  else \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, l, n); \
+          v += l*n; \
+          r += l; \
+        } \
+    } \
+} \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type *ri, \
+   octave_idx_type l, octave_idx_type n, octave_idx_type u) \
+{ \
+  if (! n) return; \
+  if (l == 1) \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, ri, n); \
+          v += n; r++; ri++; \
+        } \
+    } \
+  else \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, ri, l, n); \
+          v += l*n; \
+          r += l; ri += l; \
+        } \
+    } \
+}
+
+OP_MINMAX_FCNN (mx_inline_min)
+OP_MINMAX_FCNN (mx_inline_max)
+
 // Assistant function
 
 inline void
@@ -542,6 +699,52 @@
   return ret;
 }
 
+template <class ArrayType>
+inline ArrayType
+do_mx_minmax_op (const ArrayType& src, int dim,
+                 void (*mx_minmax_op) (const typename ArrayType::element_type *, 
+                                       typename ArrayType::element_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);
+
+  // If the dimension is zero, we don't do anything.
+  if (dim < dims.length () && dims(dim) != 0) dims(dim) = 1;
+  dims.chop_trailing_singletons ();
+
+  ArrayType ret (dims);
+  mx_minmax_op (src.data (), ret.fortran_vec (), l, n, u);
+
+  return ret;
+}
+
+template <class ArrayType>
+inline ArrayType
+do_mx_minmax_op (const ArrayType& src, Array<octave_idx_type>& idx, int dim,
+                 void (*mx_minmax_op) (const typename ArrayType::element_type *, 
+                                       typename ArrayType::element_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);
+
+  // If the dimension is zero, we don't do anything.
+  if (dim < dims.length () && dims(dim) != 0) dims(dim) = 1;
+  dims.chop_trailing_singletons ();
+
+  ArrayType ret (dims);
+  if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
+
+  mx_minmax_op (src.data (), ret.fortran_vec (), idx.fortran_vec (),
+                l, n, u);
+
+  return ret;
+}
+
 // Avoid some code duplication.  Maybe we should use templates.
 
 #define MX_CUMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \
diff -r 8af4ba6b4216 -r 9f7ce4bf7650 liboctave/oct-cmplx.h
--- a/liboctave/oct-cmplx.h	Mon Feb 16 07:57:03 2009 +0100
+++ b/liboctave/oct-cmplx.h	Mon Feb 16 08:52:00 2009 +0100
@@ -29,6 +29,24 @@
 typedef std::complex<double> Complex;
 typedef std::complex<float> FloatComplex;
 
+// The default comparison of complex number is to compare by abs, then by arg.
+// FIXME: this could be speeded up significantly.
+template <class T>
+inline bool operator < (const std::complex<T>& a,
+                        const std::complex<T>& b)
+{
+  T ax = std::abs (a), bx = std::abs (b);
+  return ax < bx || (ax == bx && std::arg (a) < std::arg (b));
+}
+
+template <class T>
+inline bool operator > (const std::complex<T>& a,
+                        const std::complex<T>& b)
+{
+  T ax = std::abs (a), bx = std::abs (b);
+  return ax > bx || (ax == bx && std::arg (a) > std::arg (b));
+}
+
 #endif
 
 /*