# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1235053199 -3600
# Node ID de16ebeef93dc628cee477332567752e5414b130
# Parent  70d06ed27c08b49b0b7b2e5bfbe60d77f3d1e110
improve lookup, provide Array<T>::lookup

diff -r 70d06ed27c08 -r de16ebeef93d liboctave/Array.cc
--- a/liboctave/Array.cc	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/Array.cc	Thu Feb 19 15:19:59 2009 +0100
@@ -2383,6 +2383,67 @@
 
 }
 
+// Do a binary lookup in a sorted array.
+template <class T>
+octave_idx_type 
+Array<T>::lookup (const T& value, sortmode mode) const
+{
+  octave_idx_type n = numel ();
+  octave_sort<T> lsort;
+
+  if (mode == UNSORTED)
+    {
+      // auto-detect mode
+      if (n > 1 && lsort.descending_compare (elem (0), elem (n-1)))
+        mode = DESCENDING;
+      else
+        mode = ASCENDING;
+    }
+
+  lsort.set_compare (mode);
+
+  return lsort.lookup (data (), n, value);
+}
+
+// Ditto, but for an array of values, specializing on long runs.
+// Adds optional offset to all indices.
+template <class T>
+Array<octave_idx_type> 
+Array<T>::lookup (const Array<T>& values, sortmode mode, 
+                  bool linf, bool rinf) const
+{
+  octave_idx_type n = numel ();
+  octave_sort<T> lsort;
+  Array<octave_idx_type> idx (values.dims ());
+
+  if (mode == UNSORTED)
+    {
+      // auto-detect mode
+      if (n > 1 && lsort.descending_compare (elem (0), elem (n-1)))
+        mode = DESCENDING;
+      else
+        mode = ASCENDING;
+    }
+
+  lsort.set_compare (mode);
+
+  // set offset and shift size.
+  octave_idx_type offset = 0;
+
+  if (linf && n > 0)
+    {
+      offset++;
+      n--;
+    }
+  if (rinf && n > 0)
+    n--;
+
+  lsort.lookup (data () + offset, n, values.data (), values.numel (),
+                idx.fortran_vec (), offset);
+
+  return idx;
+}
+
 #define INSTANTIATE_ARRAY_SORT(T) template class octave_sort<T>;
 
 #define NO_INSTANTIATE_ARRAY_SORT(T) \
@@ -2409,6 +2470,13 @@
 template <> sortmode  \
 Array<T>::is_sorted_rows (sortmode) const \
 { return UNSORTED; } \
+ \
+template <> octave_idx_type  \
+Array<T>::lookup (const T&, sortmode) const \
+{ return 0; } \
+template <> Array<octave_idx_type>  \
+Array<T>::lookup (const Array<T>&, sortmode, bool, bool) const \
+{ return Array<octave_idx_type> (); } \
 
 
 
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/Array.h
--- a/liboctave/Array.h	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/Array.h	Thu Feb 19 15:19:59 2009 +0100
@@ -567,6 +567,18 @@
   // Ordering is auto-detected or can be specified.
   sortmode is_sorted_rows (sortmode mode = UNSORTED) const;
 
+  // Do a binary lookup in a sorted array.
+  // Mode can be specified or is auto-detected by comparing 1st and last element.
+  octave_idx_type lookup (const T& value, sortmode mode = UNSORTED) const;
+
+  // Ditto, but for an array of values, specializing on long runs.
+  // If linf is true, the leftmost interval is extended to infinity 
+  // (indices will be >= 1).
+  // If rinf is true, the rightmost interval is extended to infinity 
+  // (indices will be <= length ()-1).
+  Array<octave_idx_type> lookup (const Array<T>& values, sortmode mode = UNSORTED, 
+                                 bool linf = false, bool rinf = false) const;
+
   Array<T> diag (octave_idx_type k = 0) const;
 
   template <class U, class F>
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/ChangeLog
--- a/liboctave/ChangeLog	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/ChangeLog	Thu Feb 19 15:19:59 2009 +0100
@@ -1,3 +1,16 @@
+2009-02-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* oct-types.h (sortmode): Move enum here.
+	* oct-sort.h (octave_sort<T>::ms): Declare as pointer.
+	(octave_sort<T>::lookup): New overloaded method.
+	* oct-sort.cc: Reflect change to ms.
+	(octave_sort<T>::lookup): New overloaded method.
+	(out_of_range_pred): New helper class.
+	(out_of_range): New helper function.
+	* oct-lookup.h: Remove file.
+	* Array.cc (Array<T>::lookup): New overloaded method.
+	* Array.h: Declare it.
+
 2009-02-18  John W. Eaton  <jwe@octave.org>
 
 	* dbleQR.cc (QR::init, QR::form): Cast int to octave_idx_type in
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/Makefile.in
--- a/liboctave/Makefile.in	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/Makefile.in	Thu Feb 19 15:19:59 2009 +0100
@@ -88,7 +88,7 @@
 	lo-ieee.h lo-mappers.h lo-math.h lo-specfun.h lo-sysdep.h \
 	lo-traits.h lo-utils.h mach-info.h md5.h oct-alloc.h oct-cmplx.h \
 	oct-env.h oct-fftw.h oct-getopt.h oct-group.h oct-inttypes.h \
-	oct-locbuf.h oct-lookup.h oct-md5.h oct-mutex.h oct-norm.h \
+	oct-locbuf.h oct-md5.h oct-mutex.h oct-norm.h \
         oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h oct-shlib.h \
         oct-sort.h oct-spparms.h oct-syscalls.h oct-sparse.h oct-time.h \
         oct-uname.h pathlen.h pathsearch.h prog-args.h \
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/oct-lookup.h
--- a/liboctave/oct-lookup.h	Thu Feb 19 07:34:15 2009 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,175 +0,0 @@
-/*
-
-Copyright (C) 2008 VZLU Prague, a.s., Czech Republic
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-// Author: Jaroslav Hajek <highegg@gmail.com>
-
-#if !defined (octave_oct_lookup)
-#define octave_oct_lookup 1
-
-#include <algorithm>
-#include <functional>
-
-#include "oct-types.h"
-
-// a simple binary lookup
-template<typename T, typename bpred>
-inline octave_idx_type
-bin_lookup (const T *table, octave_idx_type size, 
-            const T& val,
-            bpred comp)
-{
-  return std::upper_bound (table, table + size, val, comp) - table;
-}
-
-// version using < operator
-template<typename T>
-inline octave_idx_type
-bin_lookup (const T *table, octave_idx_type size,
-            const T& val)
-{
-  return std::upper_bound (table, table + size, val) - table;
-}
-
-// a unary functor that checks whether a value is outside [a,b) range
-template<class T, class bpred>
-class out_range : public std::unary_function<T, bool>
-{
-public:
-  out_range (const T& aa, const T& bb, const bpred& ccomp) 
-    : a(aa), b(bb), comp(ccomp) { }
-
-  bool operator() (const T& x) { return comp (x, a) || ! comp (x, b); }
-
-private:
-  T a;
-  T b;
-
-  bpred comp;
-};
-
-// conveniently constructs the above functor
-// NOTE: with SGI extensions, this can be written as
-// compose2 (logical_and<bool>(), 
-//           bind2nd (less<T>(), a),
-//           not1 (bind2nd (less<T>(), b)))
-template<class T, class bpred>
-inline out_range<T, bpred> 
-chk_out_range (const T& a, const T& b, bpred comp)
-{
-  return out_range<T, bpred> (a, b, comp);
-}
-
-template<typename T, typename bpred>
-inline void 
-seq_lookup (const T *table, octave_idx_type offset, octave_idx_type size,
-            const T *vals, octave_idx_type nvals,
-            octave_idx_type *idx, bpred comp)
-{
-  const T *begin = table + offset;
-
-  if (size == 0)
-    // the trivial case of empty table
-    std::fill_n (idx, nvals, offset);
-  else
-    {
-      const T *vcur = vals;
-      const T *vend = vals + nvals;
-
-      const T *cur = begin;
-      const T *end = begin + size;
-
-      while (vcur < vend)
-        {
-          // determine the enclosing interval for next value, trying
-          // ++cur as a special case;
-          if (cur == end || comp (*vcur, *cur))
-            cur = std::upper_bound (begin, cur, *vcur, comp);
-          else
-            {
-              ++cur;
-              if (cur < end && ! comp (*vcur, *cur))
-                cur = std::upper_bound (cur + 1, end, *vcur, comp);
-            }
-
-          // store index of the current interval.
-          *(idx++) = (cur - table);
-          ++vcur;
-
-          // find first value not in current subrange
-          const T *vnew;
-          if (cur < end)
-            if (cur > begin)
-              // inner interval
-              vnew = std::find_if (vcur, vend,
-                                   chk_out_range (*(cur-1), *cur, comp));
-
-            else
-              // special case: lowermost range (-Inf, min) 
-              vnew = std::find_if (vcur, vend,
-                                   std::not1 (std::bind2nd (comp, *cur)));
-          else
-            // special case: uppermost range [max, Inf)
-            vnew = std::find_if (vcur, vend,
-                                 std::bind2nd (comp, *(cur-1)));
-
-          // store index of the current interval.
-          std::fill_n (idx, vnew - vcur, cur - table);
-	  idx += (vnew - vcur);
-          vcur = vnew;
-
-        }
-    }
-}
-
-// overload using < operator
-template<typename T, typename bpred>
-inline void 
-seq_lookup (const T *table, octave_idx_type offset, octave_idx_type size,
-            const T *vals, octave_idx_type nvals,
-            octave_idx_type *idx)
-{
-  seq_lookup (table, offset, size, vals, nvals, idx, std::less<T>());
-}
-
-// helper functions - determine whether an array is descending
-template<typename T>
-inline bool 
-is_descending (const T *table, octave_idx_type size)
-{
-  return size > 1 && table[size-1] < table[0];
-}
-
-template<typename T, typename bpred>
-inline bool 
-is_descending (const T *table, octave_idx_type size,
-                    bpred comp)
-{
-  return size > 1 && comp (table[size-1], table[0]);
-}
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/oct-sort.cc
--- a/liboctave/oct-sort.cc	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/oct-sort.cc	Thu Feb 19 15:19:59 2009 +0100
@@ -96,6 +96,7 @@
 
 #include <cassert>
 #include <algorithm>
+#include <functional>
 #include <cstring>
 #include <stack>
 
@@ -105,15 +106,22 @@
 #include "oct-locbuf.h"
 
 template <class T>
-octave_sort<T>::octave_sort (void) : compare (ascending_compare)
+octave_sort<T>::octave_sort (void) : 
+  compare (ascending_compare), ms (0)
 { 
-  merge_init ();
 }
 
 template <class T>
-octave_sort<T>::octave_sort (compare_fcn_type comp) : compare (comp) 
+octave_sort<T>::octave_sort (compare_fcn_type comp) 
+  : compare (comp), ms (0)
 { 
-  merge_init (); 
+}
+
+template <class T>
+octave_sort<T>::~octave_sort () 
+{ 
+  merge_freemem ();
+  delete ms;
 }
 
 template <class T>
@@ -476,11 +484,12 @@
 void
 octave_sort<T>::merge_init (void)
 {
-  ms.a = 0;
-  ms.ia = 0;
-  ms.alloced = 0;
-  ms.n = 0;
-  ms.min_gallop = MIN_GALLOP;
+  if (! ms) ms = new MergeState;
+  ms->a = 0;
+  ms->ia = 0;
+  ms->alloced = 0;
+  ms->n = 0;
+  ms->min_gallop = MIN_GALLOP;
 }
 
 /* Free all the temp memory owned by the MergeState.  This must be called
@@ -491,10 +500,13 @@
 void
 octave_sort<T>::merge_freemem (void)
 {
-  delete [] ms.a;
-  delete [] ms.ia;
-  ms.alloced = 0;
-  ms.a = 0;
+  if (ms)
+    {
+      delete [] ms->a;
+      delete [] ms->ia;
+      ms->alloced = 0;
+      ms->a = 0;
+    }
 }
 
 static inline octave_idx_type
@@ -541,7 +553,7 @@
 int
 octave_sort<T>::merge_getmem (octave_idx_type need)
 {
-  if (need <= ms.alloced)
+  if (need <= ms->alloced)
     return 0;
 
   need = roundupsize (need); 
@@ -549,10 +561,10 @@
    * we don't care what's in the block.
    */
   merge_freemem ();
-  ms.a = new T[need];
-  if (ms.a)
+  ms->a = new T[need];
+  if (ms->a)
     {
-      ms.alloced = need;
+      ms->alloced = need;
       return 0;
     }
   merge_freemem ();	/* reset to sane state */
@@ -564,7 +576,7 @@
 int
 octave_sort<T>::merge_getmemi (octave_idx_type need)
 {
-  if (need <= ms.alloced && ms.a && ms.ia)
+  if (need <= ms->alloced && ms->a && ms->ia)
     return 0;
 
   need = roundupsize (need); 
@@ -572,11 +584,11 @@
    * we don't care what's in the block.
    */
   merge_freemem ();
-  ms.a = new T[need];
-  ms.ia = new octave_idx_type[need];
-  if (ms.a && ms.ia)
+  ms->a = new T[need];
+  ms->ia = new octave_idx_type[need];
+  if (ms->a && ms->ia)
     {
-      ms.alloced = need;
+      ms->alloced = need;
       return 0;
     }
   merge_freemem ();	/* reset to sane state */
@@ -600,13 +612,13 @@
   octave_idx_type k;
   T *dest;
   int result = -1;	/* guilty until proved innocent */
-  octave_idx_type min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms->min_gallop;
 
   if (merge_getmem (na) < 0)
     return -1;
-  std::copy (pa, pa + na, ms.a);
+  std::copy (pa, pa + na, ms->a);
   dest = pa;
-  pa = ms.a;
+  pa = ms->a;
 
   *dest++ = *pb++;
   --nb;
@@ -662,7 +674,7 @@
       do
 	{
 	  min_gallop -= min_gallop > 1;
-	  ms.min_gallop = min_gallop;
+	  ms->min_gallop = min_gallop;
 	  k = gallop_right (*pb, pa, na, 0, comp);
 	  acount = k;
 	  if (k)
@@ -706,7 +718,7 @@
       while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP);
 
       ++min_gallop;	/* penalize it for leaving galloping mode */
-      ms.min_gallop = min_gallop;
+      ms->min_gallop = min_gallop;
     }
 
  Succeed:
@@ -736,14 +748,14 @@
   T *dest;
   octave_idx_type *idest;
   int result = -1;	/* guilty until proved innocent */
-  octave_idx_type min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms->min_gallop;
 
   if (merge_getmemi (na) < 0)
     return -1;
-  std::copy (pa, pa + na, ms.a);
-  std::copy (ipa, ipa + na, ms.ia);
+  std::copy (pa, pa + na, ms->a);
+  std::copy (ipa, ipa + na, ms->ia);
   dest = pa; idest = ipa;
-  pa = ms.a; ipa = ms.ia;
+  pa = ms->a; ipa = ms->ia;
 
   *dest++ = *pb++; *idest++ = *ipb++;
   --nb;
@@ -796,7 +808,7 @@
       do
 	{
 	  min_gallop -= min_gallop > 1;
-	  ms.min_gallop = min_gallop;
+	  ms->min_gallop = min_gallop;
 	  k = gallop_right (*pb, pa, na, 0, comp);
 	  acount = k;
 	  if (k)
@@ -842,7 +854,7 @@
       while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP);
 
       ++min_gallop;	/* penalize it for leaving galloping mode */
-      ms.min_gallop = min_gallop;
+      ms->min_gallop = min_gallop;
     }
 
  Succeed:
@@ -883,15 +895,15 @@
   T *dest;
   int result = -1;	/* guilty until proved innocent */
   T *basea, *baseb;
-  octave_idx_type min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms->min_gallop;
 
   if (merge_getmem (nb) < 0)
     return -1;
   dest = pb + nb - 1;
-  std::copy (pb, pb + nb, ms.a);
+  std::copy (pb, pb + nb, ms->a);
   basea = pa;
-  baseb = ms.a;
-  pb = ms.a + nb - 1;
+  baseb = ms->a;
+  pb = ms->a + nb - 1;
   pa += na - 1;
 
   *dest-- = *pa--;
@@ -944,7 +956,7 @@
       do 
 	{
 	  min_gallop -= min_gallop > 1;
-	  ms.min_gallop = min_gallop;
+	  ms->min_gallop = min_gallop;
 	  k = gallop_right (*pb, basea, na, na-1, comp);
 	  if (k < 0)
 	    goto Fail;
@@ -989,7 +1001,7 @@
 	    goto Succeed;
 	} while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP);
       ++min_gallop;	/* penalize it for leaving galloping mode */
-      ms.min_gallop = min_gallop;
+      ms->min_gallop = min_gallop;
     }
 
 Succeed:
@@ -1022,17 +1034,17 @@
   int result = -1;	/* guilty until proved innocent */
   T *basea, *baseb;
   octave_idx_type *ibasea, *ibaseb;
-  octave_idx_type min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms->min_gallop;
 
   if (merge_getmemi (nb) < 0)
     return -1;
   dest = pb + nb - 1;
   idest = ipb + nb - 1;
-  std::copy (pb, pb + nb, ms.a);
-  std::copy (ipb, ipb + nb, ms.ia);
+  std::copy (pb, pb + nb, ms->a);
+  std::copy (ipb, ipb + nb, ms->ia);
   basea = pa; ibasea = ipa;
-  baseb = ms.a; ibaseb = ms.ia;
-  pb = ms.a + nb - 1; ipb = ms.ia + nb - 1;
+  baseb = ms->a; ibaseb = ms->ia;
+  pb = ms->a + nb - 1; ipb = ms->ia + nb - 1;
   pa += na - 1; ipa += na - 1;
 
   *dest-- = *pa--; *idest-- = *ipa--;
@@ -1085,7 +1097,7 @@
       do 
 	{
 	  min_gallop -= min_gallop > 1;
-	  ms.min_gallop = min_gallop;
+	  ms->min_gallop = min_gallop;
 	  k = gallop_right (*pb, basea, na, na-1, comp);
 	  if (k < 0)
 	    goto Fail;
@@ -1132,7 +1144,7 @@
 	    goto Succeed;
 	} while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP);
       ++min_gallop;	/* penalize it for leaving galloping mode */
-      ms.min_gallop = min_gallop;
+      ms->min_gallop = min_gallop;
     }
 
 Succeed:
@@ -1169,19 +1181,19 @@
   octave_idx_type na, nb;
   octave_idx_type k;
 
-  pa = data + ms.pending[i].base;
-  na = ms.pending[i].len;
-  pb = data + ms.pending[i+1].base;
-  nb = ms.pending[i+1].len;
+  pa = data + ms->pending[i].base;
+  na = ms->pending[i].len;
+  pb = data + ms->pending[i+1].base;
+  nb = ms->pending[i+1].len;
 
   /* Record the length of the combined runs; if i is the 3rd-last
    * run now, also slide over the last run (which isn't involved
    * in this merge).  The current run i+1 goes away in any case.
    */
-  ms.pending[i].len = na + nb;
-  if (i == ms.n - 3)
-    ms.pending[i+1] = ms.pending[i+2];
-  --ms.n;
+  ms->pending[i].len = na + nb;
+  if (i == ms->n - 3)
+    ms->pending[i+1] = ms->pending[i+2];
+  ms->n--;
 
   /* Where does b start in a?  Elements in a before that can be
    * ignored (already in place).
@@ -1221,21 +1233,21 @@
   octave_idx_type na, nb;
   octave_idx_type k;
 
-  pa = data + ms.pending[i].base;
-  ipa = idx + ms.pending[i].base;
-  na = ms.pending[i].len;
-  pb = data + ms.pending[i+1].base;
-  ipb = idx + ms.pending[i+1].base;
-  nb = ms.pending[i+1].len;
+  pa = data + ms->pending[i].base;
+  ipa = idx + ms->pending[i].base;
+  na = ms->pending[i].len;
+  pb = data + ms->pending[i+1].base;
+  ipb = idx + ms->pending[i+1].base;
+  nb = ms->pending[i+1].len;
 
   /* Record the length of the combined runs; if i is the 3rd-last
    * run now, also slide over the last run (which isn't involved
    * in this merge).  The current run i+1 goes away in any case.
    */
-  ms.pending[i].len = na + nb;
-  if (i == ms.n - 3)
-    ms.pending[i+1] = ms.pending[i+2];
-  --ms.n;
+  ms->pending[i].len = na + nb;
+  if (i == ms->n - 3)
+    ms->pending[i+1] = ms->pending[i+2];
+  ms->n--;
 
   /* Where does b start in a?  Elements in a before that can be
    * ignored (already in place).
@@ -1279,11 +1291,11 @@
 int
 octave_sort<T>::merge_collapse (T *data, Comp comp)
 {
-  struct s_slice *p = ms.pending;
+  struct s_slice *p = ms->pending;
 
-  while (ms.n > 1) 
+  while (ms->n > 1) 
     {
-      octave_idx_type n = ms.n - 2;
+      octave_idx_type n = ms->n - 2;
       if (n > 0 && p[n-1].len <= p[n].len + p[n+1].len) 
 	{
 	  if (p[n-1].len < p[n+1].len)
@@ -1308,11 +1320,11 @@
 int
 octave_sort<T>::merge_collapse (T *data, octave_idx_type *idx, Comp comp)
 {
-  struct s_slice *p = ms.pending;
+  struct s_slice *p = ms->pending;
 
-  while (ms.n > 1) 
+  while (ms->n > 1) 
     {
-      octave_idx_type n = ms.n - 2;
+      octave_idx_type n = ms->n - 2;
       if (n > 0 && p[n-1].len <= p[n].len + p[n+1].len) 
 	{
 	  if (p[n-1].len < p[n+1].len)
@@ -1342,11 +1354,11 @@
 int
 octave_sort<T>::merge_force_collapse (T *data, Comp comp)
 {
-  struct s_slice *p = ms.pending;
+  struct s_slice *p = ms->pending;
 
-  while (ms.n > 1) 
+  while (ms->n > 1) 
     {
-      octave_idx_type n = ms.n - 2;
+      octave_idx_type n = ms->n - 2;
       if (n > 0 && p[n-1].len < p[n+1].len)
 	--n;
       if (merge_at (n, data, comp) < 0)
@@ -1361,11 +1373,11 @@
 int
 octave_sort<T>::merge_force_collapse (T *data, octave_idx_type *idx, Comp comp)
 {
-  struct s_slice *p = ms.pending;
+  struct s_slice *p = ms->pending;
 
-  while (ms.n > 1) 
+  while (ms->n > 1) 
     {
-      octave_idx_type n = ms.n - 2;
+      octave_idx_type n = ms->n - 2;
       if (n > 0 && p[n-1].len < p[n+1].len)
 	--n;
       if (merge_at (n, data, idx, comp) < 0)
@@ -1406,8 +1418,13 @@
 octave_sort<T>::sort (T *data, octave_idx_type nel, Comp comp)
 {
   /* Re-initialize the Mergestate as this might be the second time called */
-  ms.n = 0;
-  ms.min_gallop = MIN_GALLOP;
+  if (ms)
+    {
+      ms->n = 0;
+      ms->min_gallop = MIN_GALLOP;
+    }
+  else
+    merge_init ();
 
   if (nel > 1)
     {
@@ -1437,10 +1454,10 @@
 	      n = force;
 	    }
 	  /* Push run onto pending-runs stack, and maybe merge. */
-	  assert (ms.n < MAX_MERGE_PENDING);
-	  ms.pending[ms.n].base = lo;
-	  ms.pending[ms.n].len = n;
-	  ++ms.n;
+	  assert (ms->n < MAX_MERGE_PENDING);
+	  ms->pending[ms->n].base = lo;
+	  ms->pending[ms->n].len = n;
+	  ms->n++;
 	  if (merge_collapse (data, comp) < 0)
 	    goto fail;
 	  /* Advance to find next run. */
@@ -1462,10 +1479,6 @@
 octave_sort<T>::sort (T *data, octave_idx_type *idx, octave_idx_type nel, 
                       Comp comp)
 {
-  /* Re-initialize the Mergestate as this might be the second time called */
-  ms.n = 0;
-  ms.min_gallop = MIN_GALLOP;
-
   if (nel > 1)
     {
       octave_idx_type nremaining = nel; 
@@ -1497,10 +1510,10 @@
 	      n = force;
 	    }
 	  /* Push run onto pending-runs stack, and maybe merge. */
-	  assert (ms.n < MAX_MERGE_PENDING);
-	  ms.pending[ms.n].base = lo;
-	  ms.pending[ms.n].len = n;
-	  ++ms.n;
+	  assert (ms->n < MAX_MERGE_PENDING);
+	  ms->pending[ms->n].base = lo;
+	  ms->pending[ms->n].len = n;
+	  ms->n++;
 	  if (merge_collapse (data, idx, comp) < 0)
 	    goto fail;
 	  /* Advance to find next run. */
@@ -1520,7 +1533,17 @@
 void
 octave_sort<T>::sort (T *data, octave_idx_type nel)
 {
+  /* Re-initialize the Mergestate as this might be the second time called */
+  if (ms)
+    {
+      ms->n = 0;
+      ms->min_gallop = MIN_GALLOP;
+    }
+  else
+    merge_init ();
+
   merge_getmem (1024);
+
 #ifdef INLINE_ASCENDING_SORT
   if (compare == ascending_compare)
     sort (data, nel, std::less<T> ());
@@ -1539,7 +1562,17 @@
 void
 octave_sort<T>::sort (T *data, octave_idx_type *idx, octave_idx_type nel)
 {
+  /* Re-initialize the Mergestate as this might be the second time called */
+  if (ms)
+    {
+      ms->n = 0;
+      ms->min_gallop = MIN_GALLOP;
+    }
+  else
+    merge_init ();
+
   merge_getmemi (1024);
+
 #ifdef INLINE_ASCENDING_SORT
   if (compare == ascending_compare)
     sort (data, idx, nel, std::less<T> ());
@@ -1761,6 +1794,188 @@
 }
 
 
+template <class T> template <class Comp>
+octave_idx_type 
+octave_sort<T>::lookup (const T *data, octave_idx_type nel,
+                        const T& value, Comp comp)
+{
+  return std::upper_bound (data, data + nel, value, comp) - data;
+}
+
+template <class T>
+octave_idx_type 
+octave_sort<T>::lookup (const T *data, octave_idx_type nel,
+                        const T& value)
+{
+  octave_idx_type retval = 0;
+
+#ifdef INLINE_ASCENDING_SORT
+  if (compare == ascending_compare)
+    retval = lookup (data, nel, value, std::less<T> ());
+  else
+#endif
+#ifdef INLINE_DESCENDING_SORT    
+    if (compare == descending_compare)
+      retval = lookup (data, nel, value, std::greater<T> ());
+  else
+#endif
+    if (compare)
+      retval = lookup (data, nel, value, std::ptr_fun (compare));
+
+  return retval;
+}
+
+// a unary functor that checks whether a value is outside [a,b) range
+template<class T, class Comp>
+class out_of_range_pred : public std::unary_function<T, bool>
+{
+public:
+  out_of_range_pred (const T& aa, const T& bb, Comp c) 
+    : a (aa), b (bb), comp (c) { }
+  bool operator () (const T& x) { return comp (x, a) || ! comp (x, b); }
+
+private:
+  T a, b;
+  Comp comp;
+};
+
+// a unary functor that checks whether a value is < a
+template<class T, class Comp>
+class less_than_pred : public std::unary_function<T, bool>
+{
+  typedef typename ref_param<T>::type param_type;
+public:
+  less_than_pred (param_type aa, Comp c) 
+    : a (aa), comp (c) { }
+  bool operator () (const T& x) { return comp (x, a); }
+
+private:
+  T a;
+  Comp comp;
+};
+
+// a unary functor that checks whether a value is >= a
+template<class T, class Comp>
+class greater_or_equal_pred : public std::unary_function<T, bool>
+{
+public:
+  greater_or_equal_pred (const T& aa, Comp c) 
+    : a (aa), comp (c) { }
+  bool operator () (const T& x) { return ! comp (x, a); }
+
+private:
+  T a;
+  Comp comp;
+};
+
+// conveniently constructs the above functors.
+// NOTE: with SGI extensions, this one can be written as
+// compose2 (logical_and<bool>(), bind2nd (less<T>(), a),
+//           not1 (bind2nd (less<T>(), b)))
+template<class T, class Comp>
+inline out_of_range_pred<T, Comp> 
+out_of_range (const T& a, 
+              const T& b, Comp comp)
+{
+  return out_of_range_pred<T, Comp> (a, b, comp);
+}
+
+// Note: these could be written as
+//    std::not1 (std::bind2nd (comp, *cur))
+// and
+//    std::bind2nd (comp, *(cur-1)));
+// but that doesn't work for functions with reference parameters in g++ 4.3.
+template<class T, class Comp>
+inline less_than_pred<T, Comp> 
+less_than (const T& a, Comp comp)
+{
+  return less_than_pred<T, Comp> (a, comp);
+}
+template<class T, class Comp>
+inline greater_or_equal_pred<T, Comp> 
+greater_or_equal (const T& a, Comp comp)
+{
+  return greater_or_equal_pred<T, Comp> (a, comp);
+}
+
+
+template <class T> template <class Comp>
+void 
+octave_sort<T>::lookup (const T *data, octave_idx_type nel,
+                        const T *values, octave_idx_type nvalues,
+                        octave_idx_type *idx, octave_idx_type offset, Comp comp)
+{
+  if (nel == 0)
+    // the trivial case of empty table
+    std::fill_n (idx, nvalues, offset);
+  else
+    {
+      const T *vcur = values;
+      const T *vend = values + nvalues;
+
+      const T *cur = data;
+      const T *end = data + nel;
+
+      while (vcur != vend)
+        {
+          // determine the enclosing interval for next value, trying
+          // ++cur as a special case;
+          if (cur == end || comp (*vcur, *cur))
+            cur = std::upper_bound (data, cur, *vcur, comp);
+          else
+            {
+              ++cur;
+              if (cur != end && ! comp (*vcur, *cur))
+                cur = std::upper_bound (cur + 1, end, *vcur, comp);
+            }
+
+          octave_idx_type vidx = cur - data + offset;
+          // store index of the current interval.
+          *(idx++) = vidx;
+          ++vcur;
+
+          // find first value not in current subrange
+          const T *vnew;
+          if (cur != end)
+            if (cur != data)
+              // inner interval
+              vnew = std::find_if (vcur, vend,
+                                   out_of_range (*(cur-1), *cur, comp));
+            else
+              // special case: lowermost range (-Inf, min) 
+              vnew = std::find_if (vcur, vend, greater_or_equal (*cur, comp));
+          else
+            // special case: uppermost range [max, Inf)
+            vnew = std::find_if (vcur, vend, less_than (*(cur-1), comp));
+
+          // store index of the current interval.
+          std::fill_n (idx, vnew - vcur, vidx);
+          idx += (vnew - vcur);
+          vcur = vnew;
+        }
+    }
+}
+
+template <class T>
+void 
+octave_sort<T>::lookup (const T *data, octave_idx_type nel,
+                        const T* values, octave_idx_type nvalues,
+                        octave_idx_type *idx, octave_idx_type offset)
+{
+#ifdef INLINE_ASCENDING_SORT
+  if (compare == ascending_compare)
+    lookup (data, nel, values, nvalues, idx, offset, std::less<T> ());
+  else
+#endif
+#ifdef INLINE_DESCENDING_SORT    
+    if (compare == descending_compare)
+      lookup (data, nel, values, nvalues, idx, offset, std::greater<T> ());
+  else
+#endif
+    if (compare)
+      lookup (data, nel, values, nvalues, idx, offset, std::ptr_fun (compare));
+}
+
 template <class T>
 bool 
 octave_sort<T>::ascending_compare (typename ref_param<T>::type x,
diff -r 70d06ed27c08 -r de16ebeef93d liboctave/oct-sort.h
--- a/liboctave/oct-sort.h	Thu Feb 19 07:34:15 2009 -0500
+++ b/liboctave/oct-sort.h	Thu Feb 19 15:19:59 2009 +0100
@@ -114,7 +114,7 @@
 
   octave_sort (compare_fcn_type);
   
-  ~octave_sort (void) { merge_freemem (); }
+  ~octave_sort (void); 
 
   void set_compare (compare_fcn_type comp) { compare = comp; }
 
@@ -138,6 +138,16 @@
   bool is_sorted_rows (const T *data, 
                        octave_idx_type rows, octave_idx_type cols);
 
+  // Do a binary lookup in a sorted array.
+  octave_idx_type lookup (const T *data, octave_idx_type nel,
+                          const T& value);
+
+  // Ditto, but for an array of values, specializing on long runs.
+  // Adds offset to all indices.
+  void lookup (const T *data, octave_idx_type nel,
+               const T* values, octave_idx_type nvalues,
+               octave_idx_type *idx, octave_idx_type offset = 0);
+
   static bool ascending_compare (typename ref_param<T>::type,
 				 typename ref_param<T>::type);
 
@@ -187,7 +197,7 @@
 
   compare_fcn_type compare;
   
-  MergeState ms;
+  MergeState *ms;
   
     
   template <class Comp>
@@ -277,6 +287,15 @@
   bool is_sorted_rows (const T *data, octave_idx_type rows, 
                        octave_idx_type cols, Comp comp);
 
+  template <class Comp>
+  octave_idx_type lookup (const T *data, octave_idx_type nel,
+                          const T& value, Comp comp);
+
+  template <class Comp>
+  void lookup (const T *data, octave_idx_type nel,
+               const T* values, octave_idx_type nvalues,
+               octave_idx_type *idx, octave_idx_type offset, Comp comp);
+
 };
 
 template <class T>
diff -r 70d06ed27c08 -r de16ebeef93d src/ChangeLog
--- a/src/ChangeLog	Thu Feb 19 07:34:15 2009 -0500
+++ b/src/ChangeLog	Thu Feb 19 15:19:59 2009 +0100
@@ -1,3 +1,8 @@
+2009-02-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/lookup.cc (Flookup): Use Array<T>::lookup if possible.
+	Do not compare octave_values directly. Properly check for iscellstr.
+
 2009-02-19  John W. Eaton  <jwe@octave.org>
 
 	* data.cc, graphics.cc, help.cc, lex.l, load-path.cc, parse.y:
diff -r 70d06ed27c08 -r de16ebeef93d src/DLD-FUNCTIONS/lookup.cc
--- a/src/DLD-FUNCTIONS/lookup.cc	Thu Feb 19 07:34:15 2009 -0500
+++ b/src/DLD-FUNCTIONS/lookup.cc	Thu Feb 19 15:19:59 2009 +0100
@@ -32,7 +32,6 @@
 
 #include "dNDArray.h"
 #include "CNDArray.h"
-#include "oct-lookup.h"
 
 #include "Cell.h"
 #include "defun-dld.h"
@@ -49,20 +48,6 @@
 	  || str.find (std::toupper (c)) != std::string::npos);
 }
 
-// normal ascending comparator
-static bool
-ov_str_less (const octave_value& a, const octave_value& b)
-{
-  return a.string_value () < b.string_value ();
-}
-
-// normal descending comparator
-static bool
-ov_str_greater (const octave_value& a, const octave_value& b)
-{
-  return a.string_value () > b.string_value ();
-}
-
 // case-insensitive character comparison functors
 struct icmp_char_lt : public std::binary_function<char, char, bool>
 {
@@ -76,30 +61,47 @@
     { return std::toupper (x) > std::toupper (y); }
 };
 
+// FIXME: maybe these should go elsewhere?
 // case-insensitive ascending comparator
 static bool
-ov_stri_less (const octave_value& a, const octave_value& b)
+stri_comp_lt (const std::string& a, const std::string& b)
 {
-  std::string as = a.string_value ();
-  std::string bs = b.string_value ();
-
-  return std::lexicographical_compare (as.begin (), as.end (), 
-                                       bs.begin (), bs.end (),
+  return std::lexicographical_compare (a.begin (), a.end (), 
+                                       b.begin (), b.end (),
                                        icmp_char_lt());
 }
 
 // case-insensitive descending comparator
 static bool
-ov_stri_greater (const octave_value& a, const octave_value& b)
+stri_comp_gt (const std::string& a, const std::string& b)
 {
-  std::string as = a.string_value ();
-  std::string bs = b.string_value ();
-
-  return std::lexicographical_compare (as.begin (), as.end (), 
-                                       bs.begin (), bs.end (),
+  return std::lexicographical_compare (a.begin (), a.end (), 
+                                       b.begin (), b.end (),
                                        icmp_char_gt());
 }
 
+template <class T>
+inline sortmode 
+get_sort_mode (const Array<T>& array,
+               typename octave_sort<T>::compare_fcn_type desc_comp
+               = octave_sort<T>::descending_compare)
+{
+  octave_idx_type n = array.numel ();
+  if (n > 1 && desc_comp (array (0), array (n-1)))
+    return DESCENDING;
+  else
+    return ASCENDING;
+}
+
+// FIXME: perhaps there should be octave_value::lookup?
+// The question is, how should it behave w.r.t. the second argument's type. 
+// We'd need a dispatch on two arguments. Hmmm...
+
+#define INT_ARRAY_LOOKUP(TYPE) \
+  (table.is_ ## TYPE ## _type () && y.is_ ## TYPE ## _type ()) \
+    idx = table.TYPE ## _array_value ().lookup (y.TYPE ## _array_value (), \
+                                                UNSORTED, left_inf, right_inf);
+
 DEFUN_DLD (lookup, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {@var{idx} =} lookup (@var{table}, @var{y}, @var{opt})\n\
@@ -118,14 +120,15 @@
 \n\
 The algorithm used by lookup is standard binary search, with optimizations\n\
 to speed up the case of partially ordered arrays (dense downsampling).\n\
-In particular, looking up a single entry is of binary complexity.\n\
+In particular, looking up a single entry is of logarithmic complexity\n\
+(unless a conversion occurs due to non-numeric or unequal types).\n\
 \n\
-@var{table} and @var{y} can also be a cell array of strings\n\
+@var{table} and @var{y} can also be cell arrays of strings\n\
 (or @var{y} can be a single string). In this case, string lookup\n\
 is performed using lexicographical comparison.\n\
+\n\
 If @var{opts} is specified, it shall be a string with letters indicating\n\
 additional options.\n\
-\n\
 For numeric lookup, 'l' in @var{opts} indicates that\n\
 the leftmost subinterval shall be extended to infinity (i.e. all indices\n\
 at least 1), and 'r' indicates that the rightmost subinterval shall be\n\
@@ -144,12 +147,12 @@
       return retval;
     }
 
-  octave_value argtable = args(0), argy = args(1);
-  if (argtable.ndims () > 2 || (argtable.columns () > 1 && argtable.rows () > 1))
+  octave_value table = args(0), y = args(1);
+  if (table.ndims () > 2 || (table.columns () > 1 && table.rows () > 1))
     warning ("lookup: table is not a vector");
 
-  bool num_case = argtable.is_numeric_type () && argy.is_numeric_type ();
-  bool str_case = argtable.is_cell () && (argy.is_cell () || argy.is_string ());
+  bool num_case = table.is_numeric_type () && y.is_numeric_type ();
+  bool str_case = table.is_cellstr () && (y.is_string () || y.is_cellstr ());
 
   if (num_case) 
     {
@@ -163,74 +166,42 @@
           right_inf = contains_char (opt, 'r');
         }
 
-      // in the case of a complex array, absolute values will be used for compatibility
+      // In the case of a complex array, absolute values will be used for compatibility
       // (though it's not too meaningful).
-      ArrayN<octave_idx_type> idx;
-
-      if (argtable.is_single_type () || argy.is_single_type ())
-	{
-	  FloatNDArray table = (argtable.is_complex_type ()) 
-	    ? argtable.float_complex_array_value ().abs ()
-	    : argtable.float_array_value ();
+      
+      if (table.is_complex_type ())
+        table = table.abs ();
 
-	  FloatNDArray y = (argy.is_complex_type ()) 
-	    ? argy.float_complex_array_value ().abs ()
-	    : argy.float_array_value ();
-
-	  idx = ArrayN<octave_idx_type> (y.dims ());
+      if (y.is_complex_type ())
+        y = y.abs ();
 
-	  // determine whether the array is descending. 
-	  bool desc = is_descending (table.data (), table.length ());
-	  octave_idx_type offset = left_inf ? 1 : 0;
-	  octave_idx_type size = table.length () - offset - (right_inf ? 1 : 0);
-	  if (size < 0) 
-	    size = 0;
+      Array<octave_idx_type> idx;
 
-	  if (desc)
-	    seq_lookup (table.data (), offset, size, 
-			y.data (), y.length (), idx.fortran_vec (),
-			std::greater<float> ());
-	  else
-	    seq_lookup (table.data (), offset, size, 
-			y.data (), y.length (), idx.fortran_vec (),
-			std::less<float> ());
-	}
+      // PS: I learned this from data.cc
+      if INT_ARRAY_LOOKUP (int8)
+      else if INT_ARRAY_LOOKUP (int16)
+      else if INT_ARRAY_LOOKUP (int32)
+      else if INT_ARRAY_LOOKUP (int64)
+      else if INT_ARRAY_LOOKUP (uint8)
+      else if INT_ARRAY_LOOKUP (uint16)
+      else if INT_ARRAY_LOOKUP (uint32)
+      else if INT_ARRAY_LOOKUP (uint64)
+      else if (table.is_single_type () || y.is_single_type ())
+        idx = table.float_array_value ().lookup (y.float_array_value (), 
+                                                 UNSORTED, left_inf, right_inf);
       else
-	{
-	  NDArray table = (argtable.is_complex_type ()) 
-	    ? argtable.complex_array_value ().abs ()
-	    : argtable.array_value ();
-
-	  NDArray y = (argy.is_complex_type ()) 
-	    ? argy.complex_array_value ().abs ()
-	    : argy.array_value ();
-
-	  idx = ArrayN<octave_idx_type> (y.dims ());
-
-	  // determine whether the array is descending. 
-	  bool desc = is_descending (table.data (), table.length ());
-	  octave_idx_type offset = left_inf ? 1 : 0;
-	  octave_idx_type size = table.length () - offset - (right_inf ? 1 : 0);
-	  if (size < 0) 
-	    size = 0;
-
-	  if (desc)
-	    seq_lookup (table.data (), offset, size, 
-			y.data (), y.length (), idx.fortran_vec (),
-			std::greater<double> ());
-	  else
-	    seq_lookup (table.data (), offset, size, 
-			y.data (), y.length (), idx.fortran_vec (),
-			std::less<double> ());
-	}
+        idx = table.array_value ().lookup (y.array_value (), 
+                                           UNSORTED, left_inf, right_inf);
 
       retval(0) = NDArray (idx);
+
     }
   else if (str_case)
     {
-      Cell table = argtable.cell_value ();
+      Array<std::string> str_table = table.cellstr_value ();
       
-      bool (*ov_str_comp) (const octave_value&, const octave_value&);
+      // Here we'll use octave_sort directly to avoid converting the array
+      // for case-insensitive comparison.
 
       bool icase = false;
 
@@ -241,51 +212,40 @@
           icase = contains_char (opt, 'i');
         }
 
+      sortmode mode = (icase ? get_sort_mode (str_table, stri_comp_gt)
+                       : get_sort_mode (str_table));
+
+      bool (*str_comp) (const std::string&, const std::string&);
+
       // pick the correct comparator
-      if (icase)
-        {
-          if (is_descending (table.data (), table.length (), ov_stri_less))
-            ov_str_comp = ov_stri_greater;
-          else
-            ov_str_comp = ov_stri_less;
-        }
+      if (mode == DESCENDING)
+        str_comp = icase ? stri_comp_gt : octave_sort<std::string>::descending_compare;
       else
-        {
-          if (is_descending (table.data (), table.length (), ov_str_less))
-            ov_str_comp = ov_str_greater;
-          else
-            ov_str_comp = ov_str_less;
-        }
+        str_comp = icase ? stri_comp_lt : octave_sort<std::string>::ascending_compare;
 
-
-      // query just the first cell to verify it's a string
-      if (table.is_empty () || table(0).is_string ())
+      octave_sort<std::string> lsort (str_comp);
+      if (y.is_cellstr ())
         {
-          if (argy.is_cell ())
-            {
-              Cell y = argy.cell_value ();
-              ArrayN<octave_idx_type> idx (y.dims ());
+          Array<std::string> str_y = y.cellstr_value ();
 
-
+          Array<octave_idx_type> idx (str_y.dims ());
 
-              for (int i = 0; i < y.numel (); i++)
-                  idx(i) = bin_lookup (table.data (), table.length (), y(i), 
-                                       std::ptr_fun (ov_str_comp));
+          lsort.lookup (str_table.data (), str_table.nelem (), str_y.data (),
+                        str_y.nelem (), idx.fortran_vec ());
 
-              retval(0) = NDArray (idx);
-            }
-          else
-            {
-              octave_idx_type idx;
+          retval(0) = NDArray (idx);
+        }
+      else if (y.is_string ())
+        {
+          std::string str_y = y.string_value ();
 
-              idx = bin_lookup (table.data (), table.length (), argy, 
-                                std::ptr_fun (ov_str_comp));
+          octave_idx_type idx;
 
-              retval(0) = static_cast<double> (idx);
-            }
+          lsort.lookup (str_table.data (), str_table.nelem (), &str_y,
+                        1, &idx);
+
+          retval(0) = idx;
         }
-      else
-        error("lookup: table is not a cell array of strings.");
     }
   else
     print_usage ();