octave-lyh: liboctave/oct-sort.cc comparison

comparison liboctave/oct-sort.cc @ 9362:2ebf3ca62add

use a smarter algorithm for default lookup

author	Jaroslav Hajek <highegg@gmail.com>
date	Thu, 18 Jun 2009 15:06:35 +0200
parents	9fd5c56ce57a
children	da465c405d84

comparison

equal deleted inserted replaced

-:a00e219c402d
+:2ebf3ca62add
 retval = is_sorted_rows (data, rows, cols, compare);
 return retval;
 }
+// A helper function (just to avoid repeating expressions).
+template <class T, class Comp>
+inline octave_idx_type
+lookup_impl (const T *data, octave_idx_type lo,
+octave_idx_type hi, const T& val, Comp comp)
+{
+while (lo < hi)
+{
+octave_idx_type mid = lo + ((hi-lo) >> 1);
+if (comp (val, data[mid]))
+hi = mid;
+else
+lo = mid + 1;
+}
+return lo;
+}
 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;
+return lookup_impl (data, 0, nel, value, comp);
 }
 template <class T>
 octave_idx_type
 octave_sort<T>::lookup (const T *data, octave_idx_type nel,
 retval = lookup (data, nel, value, std::ptr_fun (compare));
 return retval;
 }
-// a unary functor that checks whether a value is outside [a,b) range
+// A recursively splitting lookup of a sorted array in another sorted array.
-template<class T, class Comp>
+template <class T> template <class Comp>
-class out_of_range_pred : public std::unary_function<T, bool>
+void
-{
+octave_sort<T>::lookup_merge (const T *data, octave_idx_type lo, octave_idx_type hi,
-public:
+const T *values, octave_idx_type nvalues,
-out_of_range_pred (const T& aa, const T& bb, Comp c)
+octave_idx_type *idx, Comp comp)
-: a (aa), b (bb), comp (c) { }
+{
-bool operator () (const T& x) { return comp (x, a) || ! comp (x, b); }
+// Fake entry.
+begin:
-private:
-T a, b;
+if (hi - lo <= nvalues + 16)
-Comp comp;
+{
-};
+// Do a linear merge.
+octave_idx_type i = lo, j = 0;
-// a unary functor that checks whether a value is < a
+while (j != nvalues && i < hi)
-template<class T, class Comp>
+{
-class less_than_pred : public std::unary_function<T, bool>
+if (comp (values[j], data[i]))
-{
+idx[j++] = i;
-typedef typename ref_param<T>::type param_type;
+else
-public:
+i++;
-less_than_pred (param_type aa, Comp c)
+}
-: a (aa), comp (c) { }
+while (j != nvalues)
-bool operator () (const T& x) { return comp (x, a); }
+idx[j++] = i;
+}
-private:
+else
-T a;
+{
-Comp comp;
+// Use the ordering to split the table; look-up recursively.
-};
+switch (nvalues)
+{
-// a unary functor that checks whether a value is >= a
+case 0:
-template<class T, class Comp>
+break;
-class greater_or_equal_pred : public std::unary_function<T, bool>
+case 1:
-{
+idx[0] = lookup_impl (data,     lo,     hi, values[0], comp);
-public:
+break;
-greater_or_equal_pred (const T& aa, Comp c)
+case 2:
-: a (aa), comp (c) { }
+idx[0] = lookup_impl (data,     lo,     hi, values[0], comp);
-bool operator () (const T& x) { return ! comp (x, a); }
+idx[1] = lookup_impl (data, idx[0],     hi, values[1], comp);
+break;
-private:
+case 3:
-T a;
+idx[1] = lookup_impl (data,     lo,     hi, values[1], comp);
-Comp comp;
+idx[0] = lookup_impl (data,     lo, idx[1], values[0], comp);
-};
+idx[2] = lookup_impl (data, idx[1],     hi, values[2], comp);
+break;
-// conveniently constructs the above functors.
+case 4:
-// NOTE: with SGI extensions, this one can be written as
+idx[2] = lookup_impl (data,     lo,     hi, values[2], comp);
-// compose2 (logical_and<bool>(), bind2nd (less<T>(), a),
+idx[0] = lookup_impl (data,     lo, idx[2], values[0], comp);
-//           not1 (bind2nd (less<T>(), b)))
+idx[1] = lookup_impl (data, idx[0], idx[2], values[1], comp);
-template<class T, class Comp>
+idx[3] = lookup_impl (data, idx[2],     hi, values[3], comp);
-inline out_of_range_pred<T, Comp>
+break;
-out_of_range (const T& a,
+case 5:
-const T& b, Comp comp)
+idx[2] = lookup_impl (data,     lo,     hi, values[2], comp);
-{
+idx[0] = lookup_impl (data,     lo, idx[2], values[0], comp);
-return out_of_range_pred<T, Comp> (a, b, comp);
+idx[1] = lookup_impl (data, idx[0], idx[2], values[1], comp);
-}
+idx[3] = lookup_impl (data, idx[2],     hi, values[3], comp);
+idx[4] = lookup_impl (data, idx[3],     hi, values[4], comp);
-// Note: these could be written as
+break;
-//    std::not1 (std::bind2nd (comp, *cur))
+default:
-// and
+{
-//    std::bind2nd (comp, *(cur-1)));
+// This is a bit tricky. We do not want a normal recursion
-// but that doesn't work for functions with reference parameters in g++ 4.3.
+// because we split by idx[k], and there's no guaranteed descend
-template<class T, class Comp>
+// ratio; hence the worst-case scenario would be a linear stack
-inline less_than_pred<T, Comp>
+// growth, which is dangerous. Instead, we recursively run the
-less_than (const T& a, Comp comp)
+// smaller half, and simulate a tail call for the rest using
-{
+// goto.
-return less_than_pred<T, Comp> (a, comp);
+octave_idx_type k = nvalues / 2;
-}
+idx[k] = lookup_impl (data, lo, hi, values[k], comp);
-template<class T, class Comp>
+if (idx[k] - lo <= hi - idx[k])
-inline greater_or_equal_pred<T, Comp>
+{
-greater_or_equal (const T& a, Comp comp)
+// The smaller portion is run recursively.
-{
+lookup_merge (data, idx[k], k, values, k, idx, comp);
-return greater_or_equal_pred<T, Comp> (a, comp);
+// Simulate a tail call.
-}
+lo = idx[k];
+values += k; nvalues -= k; idx += k;
+goto begin;
+}
+else
+{
+// The smaller portion is run recursively.
+lookup_merge (data, idx[k], hi,
+values + k, nvalues - k, idx, comp);
+// Simulate a tail call.
+hi = idx[k];
+nvalues = k;
+goto begin;
+}
+break;
+}
+}
+}
+}
 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,
 if (nel == 0)
 // the trivial case of empty table
 std::fill_n (idx, nvalues, offset);
 else
 {
-const T *vcur = values;
+octave_idx_type vlo = 0;
-const T *vend = values + nvalues;
+int nbadruns = 0;
+while (vlo != nvalues && nbadruns < 16)
-const T *cur = data;
-const T *end = data + nel;
-while (vcur != vend)
 {
-// determine the enclosing interval for next value, trying
+octave_idx_type vhi;
-// ++cur as a special case;
-if (cur == end || comp (*vcur, *cur))
+// Determine a sorted run.
-cur = std::upper_bound (data, cur, *vcur, comp);
+for (vhi = vlo + 1; vhi != nvalues; vhi++)
-else
 {
-++cur;
+if (! comp (values[vhi-1], values[vhi]))
-if (cur != end && ! comp (*vcur, *cur))
+break;
-cur = std::upper_bound (cur + 1, end, *vcur, comp);
 }
-octave_idx_type vidx = cur - data + offset;
+// Do a recursive lookup.
-// store index of the current interval.
+lookup_merge (data - offset, offset, nel + offset,
-*(idx++) = vidx;
+values + vlo, vhi - vlo, idx + vlo, comp);
-++vcur;
+if (vhi - vlo <= 2)
-// find first value not in current subrange
+nbadruns++;
-const T *vnew;
+else if (vhi - vlo >= 16)
-if (cur != end)
+nbadruns = 0;
-if (cur != data)
-// inner interval
+vlo = vhi;
-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;
 }
+// The optimistic loop terminated, do the rest normally.
+for (; vlo != nvalues; vlo++)
+idx[vlo] = lookup_impl (data, 0, nel, values[vlo], comp) + offset;
 }
 }
 template <class T>
 void
 void
 octave_sort<T>::lookupm (const T *data, octave_idx_type nel,
 const T *values, octave_idx_type nvalues,
 octave_idx_type *idx, Comp comp)
 {
-const T *end = data + nel;
 for (octave_idx_type i = 0; i < nvalues; i++)
 {
-const T *ptr = std::lower_bound (data, end, values[i], comp);
+octave_idx_type j = lookup_impl (data, 0, nel, values[i], comp);
-if (ptr != end && ! comp (values[i], *ptr))
+idx[i] = (j != 0 && comp (data[j-1], values[i])) ? -1 : j-1;
-idx[i] = ptr - data;
-else
-idx[i] = -1;
 }
 }
 template <class T>
 void
 void
 octave_sort<T>::lookupb (const T *data, octave_idx_type nel,
 const T *values, octave_idx_type nvalues,
 bool *match, Comp comp)
 {
-const T *end = data + nel;
 for (octave_idx_type i = 0; i < nvalues; i++)
-match[i] = std::binary_search (data, end, values[i], comp);
+{
+octave_idx_type j = lookup_impl (data, 0, nel, values[i], comp);
+match[i] = (j != 0 && ! comp (data[j-1], values[i]));
+}
 }
 template <class T>
 void
 octave_sort<T>::lookupb (const T *data, octave_idx_type nel,

Mercurial > hg > octave-lyh

comparison liboctave/oct-sort.cc @ 9362:2ebf3ca62add