--- a/liboctave/Array-d.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-d.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,352 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+#if defined (HAVE_IEEE754_DATA_FORMAT)
+
+static inline uint64_t
+FloatFlip (uint64_t f)
+{
+  uint64_t mask
+    = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL;
+
+  return f ^ mask;
+}
+
+static inline uint64_t
+IFloatFlip (uint64_t f)
+{
+  uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL;
+
+  return f ^ mask;
+}
+
+template <>
+bool
+ascending_compare (double a, double b)
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+ascending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec < b->vec));
+}
+
+template <>
+bool
+descending_compare (double a, double b)
+{
+  return (xisnan (a) || (a > b));
+}
+
+template <>
+bool
+descending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec > b->vec));
+}
+
+INSTANTIATE_ARRAY_SORT (uint64_t);
+
+template <>
+Array<double>
+Array<double>::sort (octave_idx_type dim, sortmode mode) const
+{
+  Array<double> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    return m;
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  double *v = m.fortran_vec ();
+
+  uint64_t *p = reinterpret_cast<uint64_t *> (v);
+
+  octave_sort<uint64_t> lsort;
+
+  if (mode == ASCENDING)
+    lsort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    lsort.set_compare (descending_compare);
+
+  if (stride == 1)
+    {
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  // Flip the data in the vector so that int compares on
+	  // IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i] = FloatFlip (p[i]);
+	      
+	  lsort.sort (p, ns);
+
+	  // Flip the data out of the vector so that int compares
+	  // on IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i] = IFloatFlip (p[i]);
+
+	  // There are two representations of NaN.  One will be
+	  // sorted to the beginning of the vector and the other
+	  // to the end.  If it will be sorted incorrectly, fix
+	  // things up.
+
+	  if (lo_ieee_signbit (octave_NaN))
+	    {
+	      if (mode == UNDEFINED || mode == ASCENDING)
+		{
+		  octave_idx_type i = 0;
+		  double *vtmp = reinterpret_cast<double *> (p);
+		  while (xisnan (vtmp[i++]) && i < ns);
+		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		    vtmp[l] = vtmp[l+i-1];
+		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
+		    vtmp[l] = octave_NaN;
+		}
+	      else
+		{
+		  octave_idx_type i = ns;
+		  double *vtmp = reinterpret_cast<double *> (p);
+		  while (xisnan (vtmp[--i]) && i > 0);
+		  for (octave_idx_type l = i; l >= 0; l--)
+		    vtmp[l-i+ns-1] = vtmp[l];
+		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
+		    vtmp[l] = octave_NaN;
+		}
+	    }
+
+	  p += ns;
+	}
+    }
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns);
+
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  octave_idx_type offset = j;
+	  octave_idx_type offset2 = 0;
+	  while (offset >= stride)
+	    {
+	      offset -= stride;
+	      offset2++;
+	    }
+	  offset += offset2 * stride * ns;
+
+	  // Flip the data in the vector so that int compares on
+	  // IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    vi[i] = FloatFlip (p[i*stride + offset]);
+
+	  lsort.sort (vi, ns);
+
+	  // Flip the data out of the vector so that int compares
+	  // on IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i*stride + offset] = IFloatFlip (vi[i]);
+	      
+	  // There are two representations of NaN. One will be
+	  // sorted to the beginning of the vector and the other
+	  // to the end. If it will be sorted to the beginning,
+	  // fix things up.
+
+	  if (lo_ieee_signbit (octave_NaN))
+	    {
+	      if (mode == UNDEFINED || mode == ASCENDING)
+		{
+		   octave_idx_type i = 0;
+		  while (xisnan (v[i++*stride + offset]) && i < ns);
+		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		    v[l*stride + offset] = v[(l+i-1)*stride + offset];
+		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
+		    v[l*stride + offset] = octave_NaN;
+		}
+	      else
+		{
+		   octave_idx_type i = ns;
+		  while (xisnan (v[--i*stride + offset]) && i > 0);
+		  for (octave_idx_type l = i; l >= 0; l--)
+		    v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
+		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
+		    v[l*stride + offset] = octave_NaN;
+		}
+	    }
+	}
+    }
+
+  return m;
+}
+
+template <>
+Array<double>
+Array<double>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, 
+		     sortmode mode) const
+{
+  Array<double> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    {
+      sidx = Array<octave_idx_type> (dv);
+      return m;
+    }
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  double *v = m.fortran_vec ();
+
+  uint64_t *p = reinterpret_cast<uint64_t *> (v);
+
+  octave_sort<vec_index<uint64_t> *> indexed_sort;
+
+  if (mode == ASCENDING)
+    indexed_sort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    indexed_sort.set_compare (descending_compare);
+
+  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns);
+  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns);
+  
+  for (octave_idx_type i = 0; i < ns; i++)
+    vi[i] = &vix[i];
+
+  sidx = Array<octave_idx_type> (dv);
+      
+  for (octave_idx_type j = 0; j < iter; j++)
+    {
+      octave_idx_type offset = j;
+      octave_idx_type offset2 = 0;
+      while (offset >= stride)
+	{
+	  offset -= stride;
+	  offset2++;
+	}
+      offset += offset2 * stride * ns;
+
+      // Flip the data in the vector so that int compares on
+      // IEEE754 give the correct ordering.
+
+      for (octave_idx_type i = 0; i < ns; i++)
+	{
+	  vi[i]->vec = FloatFlip (p[i*stride + offset]);
+	  vi[i]->indx = i;
+	}
+
+      indexed_sort.sort (vi, ns);
+
+      // Flip the data out of the vector so that int compares on
+      // IEEE754 give the correct ordering
+
+      for (octave_idx_type i = 0; i < ns; i++)
+	{
+	  p[i*stride + offset] = IFloatFlip (vi[i]->vec);
+	  sidx(i*stride + offset) = vi[i]->indx;
+	}
+
+      // There are two representations of NaN.  One will be sorted
+      // to the beginning of the vector and the other to the end.
+      // If it will be sorted to the beginning, fix things up.
+
+      if (lo_ieee_signbit (octave_NaN))
+	{
+	  if (mode == UNDEFINED || mode == ASCENDING)
+	    {
+	      octave_idx_type i = 0;
+	      while (xisnan (v[i++*stride+offset]) && i < ns);
+	      OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
+	      for (octave_idx_type l = 0; l < i -1; l++)
+		itmp[l] = sidx(l*stride + offset);
+	      for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		{
+		  v[l*stride + offset] = v[(l+i-1)*stride + offset];
+		  sidx(l*stride + offset) = sidx((l+i-1)*stride + offset);
+		}
+	      for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++)
+		{
+		  v[l*stride + offset] = octave_NaN;
+		  sidx(l*stride + offset) = 
+		    static_cast<octave_idx_type>(itmp[k]);
+		}
+	    }
+	  else 
+	    {
+	      octave_idx_type i = ns;
+	      while (xisnan (v[--i*stride+offset]) && i > 0);
+	      OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1);
+	      for (octave_idx_type l = 0; l < ns - i -1; l++)
+		itmp[l] = sidx((l+i+1)*stride + offset);
+	      for (octave_idx_type l = i; l >= 0; l--)
+		{
+		  v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
+		  sidx((l-i+ns-1)*stride + offset) = sidx(l*stride + offset);
+		}
+	      for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++)
+		{
+		  v[l*stride + offset] = octave_NaN;
+		  sidx(l*stride + offset) = 
+		    static_cast<octave_idx_type>(itmp[k]);
+		}
+	    }
+	}
+    }
+
+  return m;
+}
+
+#else
+
+template <>
+bool
+Array<double>::ascending_compare (double a, double b) const
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+Array<double>::ascending_compare (vec_index<double> *a, 
+				  vec_index<double> *b) const
+{
+  return (xisnan (b->vec) || (a->vec < b->vec));
+}
+
+template <>
+bool
+Array<double>::descending_compare (double a, double b) const
+{
+  return (xisnan (a) || (a > b));
+}
+
+template <>
+bool
+Array<double>::descending_compare (vec_index<double> *a, 
+				   vec_index<double> *b) const
+{
+  return (xisnan (b->vec) || (a->vec > b->vec));
+}
+
+INSTANTIATE_ARRAY_SORT (double);
+
+#endif
 
 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API);