use red-black tree in cdirect instead of repeatedly resorting ... convex_hull...

use red-black tree in cdirect instead of repeatedly resorting ... convex_hull is still O(N), though, and needs to be changed to a dynamic-hull algorithm to reduce the overall complexity below O(N^2), ugh

darcs-hash:20070827144928-c8de0-1d7b7b1b20a5f1ca81fc7e9b31a79ee2bd4479c8.gz

cdirect/Makefile.am

 AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
 
 noinst_LTLIBRARIES = libcdirect.la
-libcdirect_la_SOURCES = cdirect.c cdirect.h
+libcdirect_la_SOURCES = cdirect.c cdirect.h redblack.c redblack.h
+
+noinst_PROGRAMS = redblack_test
+redblack_test_SOURCES = redblack_test.c
+redblack_test_LDADD = libcdirect.la
 
 EXTRA_DIST = 

cdirect/cdirect.c

@@ -5,6 +5,7 @@
 #include "nlopt-util.h"
 #include "nlopt.h"
 #include "cdirect.h"
+#include "redblack.h"
 #include "config.h"
 
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
@@ -13,13 +14,13 @@
 /***************************************************************************/
 /* basic data structure:
  *
- * a hyper-rectangle is stored as an array of length 2n+2, where [0]
- * is the value of the function at the center, [1] is the "size"
- * measure of the rectangle, [2..n+1] are the coordinates of the
- * center, and [n+2..2n+1] are the widths of the sides.
+ * a hyper-rectangle is stored as an array of length L = 2n+2, where [0]
+ * is the value (f) of the function at the center, [1] is the "size"
+ * measure (d) of the rectangle, [2..n+1] are the coordinates of the
+ * center (c), and [n+2..2n+1] are the widths of the sides (w).
  *
  * a list of rectangles is just an array of N hyper-rectangles
- * stored as an N x (2n+1) in row-major order.  Generally,
+ * stored as an N x L in row-major order.  Generally,
  * we allocate more than we need, allocating Na hyper-rectangles.
  *
  * n > 0 always
@@ -30,6 +31,9 @@
 /* parameters of the search algorithm and various information that
    needs to be passed around */
 typedef struct {
+     int n; /* dimension */
+     int L; /* RECT_LEN(n) */
+     double *rects; /* the hyper-rectangles */
      double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
      int which_diam; /* which measure of hyper-rectangle diam to use:
 			0 = Jones, 1 = Gablonsky */
@@ -39,11 +43,14 @@ typedef struct {
 		       2: Jones Encyc. Opt.: pick random longest side */
      const double *lb, *ub;
      nlopt_stopping *stop; /* stopping criteria */
-     int n;
      nlopt_func f; void *f_data;
      double *work; /* workspace, of length >= 2*n */
      int *iwork, iwork_len; /* workspace, of length iwork_len >= n */
      double fmin, *xmin; /* minimum so far */
+     
+     /* red-black tree of hyperrect indices, sorted by (d,f) in
+	lexographical order */
+     rb_tree rtree;
 } params;
 
 /***************************************************************************/
@@ -115,18 +122,18 @@ static double function_eval(const double *x, params *p) {
 
 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 
-/* divide rectangle idiv in the list rects */
-static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
-				params *p)
+/* divide rectangle idiv in the list p->rects */
+static nlopt_result divide_rect(int *N, int *Na, int idiv, params *p)
 {
      int i;
      const const int n = p->n;
-     const int L = RECT_LEN(n);
-     double *r = *rects;
+     const int L = p->L;
+     double *r = p->rects;
      double *c = r + L*idiv + 2; /* center of rect to divide */
      double *w = c + n; /* widths of rect to divide */
      double wmax = w[0];
      int imax = 0, nlongest = 0;
+     rb_node *node;
 
      for (i = 1; i < n; ++i)
 	  if (w[i] > wmax)
@@ -154,15 +161,20 @@ static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
 	  }
 	  sort_fv(n, fv, isort);
 	  ALLOC_RECTS(n, Na, r, (*N)+2*nlongest); 
-	  *rects = r; c = r + L*idiv + 2; w = c + n;
+	  p->rects = r; c = r + L*idiv + 2; w = c + n;
 	  for (i = 0; i < nlongest; ++i) {
 	       int k;
+	       if (!(node = rb_tree_find(&p->rtree, idiv)))
+		    return NLOPT_FAILURE;
 	       w[isort[i]] *= THIRD;
 	       r[L*idiv + 1] = rect_diameter(n, w, p);
+	       rb_tree_resort(&p->rtree, node);
 	       for (k = 0; k <= 1; ++k) {
 		    memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 		    r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
 		    r[L*(*N)] = fv[2*isort[i]+k];
+		    if (!rb_tree_insert(&p->rtree, *N))
+			 return NLOPT_FAILURE;
 		    ++(*N);
 	       }
 	  }
@@ -181,13 +193,18 @@ static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
 	  else
 	       i = imax; /* trisect longest side */
 	  ALLOC_RECTS(n, Na, r, (*N)+2);
-          *rects = r; c = r + L*idiv + 2; w = c + n;
+          p->rects = r; c = r + L*idiv + 2; w = c + n;
+	  if (!(node = rb_tree_find(&p->rtree, idiv)))
+	       return NLOPT_FAILURE;
 	  w[i] *= THIRD;
 	  r[L*idiv + 1] = rect_diameter(n, w, p);
+	  rb_tree_resort(&p->rtree, node);
 	  for (k = 0; k <= 1; ++k) {
 	       memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 	       r[L*(*N) + 2 + i] += w[i] * (2*k-1);
 	       FUNCTION_EVAL(r[L*(*N)], r + L*(*N) + 2, p);
+	       if (!rb_tree_insert(&p->rtree, *N))
+		    return NLOPT_FAILURE;
 	       ++(*N);
 	  }
      }
@@ -198,72 +215,50 @@ static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
 /* O(N log N) convex hull algorithm, used later to find the potentially
    optimal points */
 
-/* sort ihull by xy in lexographic order by x,y */
-static int s_qsort = 1; static double *xy_qsort = 0;
-static int sort_xy_compare(const void *a_, const void *b_)
-{
-     int a = *((const int *) a_), b = *((const int *) b_);
-     double xa = xy_qsort[a*s_qsort+1], xb = xy_qsort[b*s_qsort+1];
-     if (xa < xb) return -1;
-     else if (xb < xa) return +1;
-     else {
-	  double ya = xy_qsort[a*s_qsort], yb = xy_qsort[b*s_qsort];
-	  if (ya < yb) return -1;
-	  else if (ya > yb) return +1;
-	  else return 0;
-     }
-}
-static void sort_xy(int N, double *xy, int s, int *isort)
-{
-     int i;
-
-     for (i = 0; i < N; ++i) isort[i] = i;
-     s_qsort = s; xy_qsort = xy;
-     qsort(isort, (unsigned) N, sizeof(int), sort_xy_compare);
-     xy_qsort = 0;
-}
-
 /* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
    0 <= i < N and s >= 2.
 
    The return value is the number of points in the hull, with indices
-   stored in ihull.  ihull and is should point to arrays of length >= N.
-
-   Note that we don't allow redundant points along the same line in the
-   hull, similar to Gablonsky's version of DIRECT and differing from
-   Jones'. */
-static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
+   stored in ihull.  ihull should point to arrays of length >= N.
+   rb_tree should be a red-black tree of indices (keys == i) sorted
+   in lexographic order by (xy[s*i+1], xy[s*i]).
+*/
+static int convex_hull(int N, double *xy, int s, int *ihull, rb_tree *t)
 {
-     int minmin; /* first index (smallest y) with min x */
-     int minmax; /* last index (largest y) with min x */
-     int maxmin; /* first index (smallest y) with max x */
-     int maxmax; /* last index (largest y) with max x */
-     int i, nhull = 0;
+     int nhull;
      double minslope;
-     double xmin, xmax;
+     double xmin, xmax, yminmin, ymaxmin;
+     rb_node *n, *nmax;
 
+     if (N == 0) return 0;
+     
      /* Monotone chain algorithm [Andrew, 1979]. */
 
-     sort_xy(N, xy, s, is);
+     n = rb_tree_min(t);
+     nmax = rb_tree_max(t);
 
-     xmin = xy[s*is[minmin=0]+1]; xmax = xy[s*is[maxmax=N-1]+1];
+     xmin = xy[s*(n->k)+1];
+     yminmin = xy[s*(n->k)];
+     xmax = xy[s*(nmax->k)+1];
 
-     if (xmin == xmax) { /* degenerate case */
-	  ihull[nhull++] = is[minmin];
-	  return nhull;
-     }
+     ihull[nhull = 1] = n->k;
+     if (xmin == xmax) return nhull;
 
-     for (minmax = minmin; minmax+1 < N && xy[s*is[minmax+1]+1]==xmin; 
-	  ++minmax);
-     for (maxmin = maxmax; maxmin-1>=0 && xy[s*is[maxmin-1]+1]==xmax;
-	  --maxmin);
+     /* set nmax = min mode with x == xmax */
+     while (xy[s*(nmax->k)+1] == xmax)
+	  nmax = rb_tree_pred(t, nmax); /* non-NULL since xmin != xmax */
+     nmax = rb_tree_succ(t, nmax);
 
-     minslope = (xy[s*is[maxmin]] - xy[s*is[minmin]]) / (xmax - xmin);
+     ymaxmin = xy[s*(nmax->k)];
+     minslope = (ymaxmin - yminmin) / (xmax - xmin);
 
-     ihull[nhull++] = is[minmin];
-     for (i = minmax + 1; i < maxmin; ++i) {
-	  int k = is[i];
-	  if (xy[s*k] > xy[s*is[minmin]] + (xy[s*k+1] - xmin) * minslope)
+     /* set n = first node with x != xmin */
+     while (xy[s*(n->k)+1] == xmin)
+	  n = rb_tree_succ(t, n); /* non-NULL since xmin != xmax */
+
+     for (; n != nmax; n = rb_tree_succ(t, n)) { 
+	  int k = n->k;
+	  if (xy[s*k] > yminmin + (xy[s*k+1] - xmin) * minslope)
 	       continue;
 	  /* remove points until we are making a "left turn" to k */
 	  while (nhull > 1) {
@@ -276,7 +271,7 @@ static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
 	  }
 	  ihull[nhull++] = k;
      }
-     ihull[nhull++] = is[maxmin];
+     ihull[nhull++] = nmax->k;
      return nhull;
 }
 
@@ -292,23 +287,22 @@ static int small(double *w, params *p)
      return 1;
 }
 
-static nlopt_result divide_good_rects(int *N, int *Na, double **rects, 
-				      params *p)
+static nlopt_result divide_good_rects(int *N, int *Na, params *p)
 {
      const int n = p->n;
-     const int L = RECT_LEN(n);
+     const int L = p->L;
      int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
-     double *r = *rects;
+     double *r = p->rects;
      double magic_eps = p->magic_eps;
 
-     if (p->iwork_len < n + 2*(*N)) {
-	  p->iwork_len = p->iwork_len + n + 2*(*N);
+     if (p->iwork_len < *N) {
+	  p->iwork_len = p->iwork_len + *N;
 	  p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
 	  if (!p->iwork)
 	       return NLOPT_OUT_OF_MEMORY;
      }
      ihull = p->iwork;
-     nhull = convex_hull(*N, r, L, ihull, ihull + *N);
+     nhull = convex_hull(*N, r, L, ihull, &p->rtree);
  divisions:
      for (i = 0; i < nhull; ++i) {
 	  double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
@@ -324,8 +318,8 @@ static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
 	       /* "potentially optimal" rectangle, so subdivide */
 	       divided_some = 1;
 	       nlopt_result ret;
-	       ret = divide_rect(N, Na, rects, ihull[i], p);
-	       r = *rects; /* may have grown */
+	       ret = divide_rect(N, Na, ihull[i], p);
+	       r = p->rects; /* may have grown */
 	       if (ret != NLOPT_SUCCESS) return ret;
 	       xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
 	  }
@@ -341,7 +335,7 @@ static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
 	       for (i = 1; i < *N; ++i)
 		    if (r[L*i+1] > wmax)
 			 wmax = r[L*(imax=i)+1];
-	       return divide_rect(N, Na, rects, imax, p);
+	       return divide_rect(N, Na, imax, p);
 	  }
      }
      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
@@ -349,6 +343,23 @@ static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
 
 /***************************************************************************/
 
+/* lexographic sort order (d,f) of hyper-rects, for red-black tree */
+static int hyperrect_compare(int i, int j, void *p_)
+{
+     params *p = (params *) p_;
+     int L = p->L;
+     double *r = p->rects;
+     double di = r[i*L+1], dj = r[j*L+1], fi, fj;
+     if (di < dj) return -1;
+     if (dj < di) return +1;
+     fi = r[i*L]; fj = r[j*L];
+     if (fi < fj) return -1;
+     if (fj < fi) return +1;
+     return 0;
+}
+
+/***************************************************************************/
+
 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 			      const double *lb, const double *ub,
 			      double *x,
@@ -357,7 +368,6 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 			      double magic_eps, int which_alg)
 {
      params p;
-     double *rects;
      int Na = 100, N = 1, i, x_center = 1;
      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 
@@ -367,38 +377,45 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
      p.lb = lb; p.ub = ub;
      p.stop = stop;
      p.n = n;
+     p.L = RECT_LEN(n);
      p.f = f;
      p.f_data = f_data;
      p.xmin = x;
      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
      p.work = 0;
      p.iwork = 0;
-     rects = 0;
+     p.rects = 0;
+
+     if (!rb_tree_init(&p.rtree, hyperrect_compare, &p)) goto done;
      p.work = (double *) malloc(sizeof(double) * 2*n);
      if (!p.work) goto done;
-     rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
-     if (!rects) goto done;
-     p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = 2*Na + n));
+     p.rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
+     if (!p.rects) goto done;
+     p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = Na + n));
      if (!p.iwork) goto done;
 
      for (i = 0; i < n; ++i) {
-	  rects[2+i] = 0.5 * (lb[i] + ub[i]);
+	  p.rects[2+i] = 0.5 * (lb[i] + ub[i]);
 	  x_center = x_center
-	       && (fabs(rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
-	  rects[2+n+i] = ub[i] - lb[i];
+	       && (fabs(p.rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
+	  p.rects[2+n+i] = ub[i] - lb[i];
      }
-     rects[1] = rect_diameter(n, rects+2+n, &p);
+     p.rects[1] = rect_diameter(n, p.rects+2+n, &p);
      if (x_center)
-	  rects[0] = p.fmin; /* avoid computing f(center) twice */
+	  p.rects[0] = p.fmin; /* avoid computing f(center) twice */
      else
-	  rects[0] = function_eval(rects+2, &p);
+	  p.rects[0] = function_eval(p.rects+2, &p);
+     if (!rb_tree_insert(&p.rtree, 0)) {
+	  ret = NLOPT_FAILURE;
+	  goto done;
+     }
 
-     ret = divide_rect(&N, &Na, &rects, 0, &p);
+     ret = divide_rect(&N, &Na, 0, &p);
      if (ret != NLOPT_SUCCESS) goto done;
 
      while (1) {
 	  double fmin0 = p.fmin;
-	  ret = divide_good_rects(&N, &Na, &rects, &p);
+	  ret = divide_good_rects(&N, &Na, &p);
 	  if (ret != NLOPT_SUCCESS) goto done;
 	  if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 	       ret = NLOPT_FTOL_REACHED;
@@ -407,8 +424,9 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
      }
 
  done:
+     rb_tree_destroy(&p.rtree);
      free(p.iwork);
-     free(rects);
+     free(p.rects);
      free(p.work);
 	      
      *fmin = p.fmin;

cdirect/redblack.c

+/* simple implementation of red-black trees optimized for use with DIRECT */
+
+#include <stddef.h>
+#include <stdlib.h>
+#include "redblack.h"
+
+int rb_tree_init(rb_tree *t, rb_compare compare, void *c_data) {
+     t->compare = compare; t->c_data = c_data;
+     t->nil.c = BLACK; t->nil.l = t->nil.r = t->nil.p = &t->nil;
+     t->root = &t->nil;
+     t->N = 0;
+     t->Nalloc = 100; /* allocate some space to start with */
+     t->nodes = (rb_node*) malloc(sizeof(rb_node) * t->Nalloc);
+     return t->nodes != NULL;
+}
+
+void rb_tree_destroy(rb_tree *t)
+{
+     t->root = 0; t->N = 0; t->Nalloc = 0;
+     free(t->nodes); t->nodes = 0;
+}
+
+/* in our application, we can optimize memory allocation because
+   we never delete two nodes in a row (we always add a node after
+   deleting)... or rather, we never delete but the value of
+   the key sometimes changes.  ... this means we can just
+   allocate a linear, exponentially growing stack (nodes) of
+   nodes, and don't have to worry about holes in the stack ...
+   otherwise, alloc1 should be replaced by an implementation that
+   malloc's each node separately */
+static rb_node *alloc1(rb_tree *t, int k)
+{
+     rb_node *nil = &t->nil;
+     rb_node *n;
+     if (t->Nalloc == t->N) { /* grow allocation */
+	  rb_node *old_nodes = t->nodes;
+	  ptrdiff_t change;
+	  int i;
+	  t->Nalloc = 2*t->Nalloc + 1;
+	  t->nodes = (rb_node*) realloc(t->nodes, sizeof(rb_node) * t->Nalloc);
+	  if (!t->nodes) return NULL;
+	  change = t->nodes - old_nodes;
+	  if (t->root != nil) t->root += change;
+	  for (i = 0; i < t->N; ++i) { /* shift all pointers, ugh */
+	       if (t->nodes[i].p != nil) t->nodes[i].p += change;
+	       if (t->nodes[i].r != nil) t->nodes[i].r += change;
+	       if (t->nodes[i].l != nil) t->nodes[i].l += change;
+	  }
+     }
+     n = t->nodes + t->N++;
+     n->k = k;
+     n->p = n->l = n->r = nil;
+     return n;
+}
+
+static void rotate_left(rb_node *p, rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_node *n = p->r; /* must be non-NULL */
+     p->r = n->l;
+     n->l = p;
+     if (p->p != nil) {
+	  if (p == p->p->l) p->p->l = n;
+	  else p->p->r = n;
+     }
+     else
+	  t->root = n;
+     n->p = p->p;
+     p->p = n;
+     if (p->r != nil) p->r->p = p;
+}
+
+static void rotate_right(rb_node *p, rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_node *n = p->l; /* must be non-NULL */
+     p->l = n->r;
+     n->r = p;
+     if (p->p != nil) {
+	  if (p == p->p->l) p->p->l = n;
+	  else p->p->r = n;
+     }
+     else
+	  t->root = n;
+     n->p = p->p;
+     p->p = n;
+     if (p->l != nil) p->l->p = p;
+}
+
+static void insert_node(rb_tree *t, rb_node *n)
+{
+     rb_node *nil = &t->nil;
+     rb_compare compare = t->compare;
+     void *c_data = t->c_data;
+     int k = n->k;
+     rb_node *p = t->root;
+     n->c = RED;
+     if (p == nil) {
+	  t->root = n;
+	  n->c = BLACK;
+	  return;
+     }
+     /* insert (RED) node into tree */
+     while (1) {
+	  if (compare(k, p->k, c_data) <= 0) { /* k <= p->k */
+	       if (p->l != nil)
+		    p = p->l;
+	       else {
+		    p->l = n;
+		    n->p = p;
+		    break;
+	       }
+	  }
+	  else {
+	       if (p->r != nil)
+		    p = p->r;
+	       else {
+		    p->r = n;
+		    n->p = p;
+		    break;
+	       }
+	  }
+     }
+ fixtree:
+     if (n->p->c == RED) { /* red cannot have red child */
+	  rb_node *u = p == p->p->l ? p->p->r : p->p->l;
+	  if (u != nil && u->c == RED) {
+	       p->c = u->c = BLACK;
+	       n = p->p;
+	       if ((p = n->p) != nil) {
+		    n->c = RED;
+		    goto fixtree;
+	       }
+	  }
+	  else {
+	       if (n == p->r && p == p->p->l) {
+		    rotate_left(p, t);
+		    p = n; n = n->l;
+	       }
+	       else if (n == p->l && p == p->p->r) {
+		    rotate_right(p, t);
+		    p = n; n = n->r;
+	       }
+	       p->c = BLACK;
+	       p->p->c = RED;
+	       if (n == p->l && p == p->p->l)
+		    rotate_right(p->p, t);
+	       else if (n == p->r && p == p->p->r)
+		    rotate_left(p->p, t);
+	  }
+	      
+     }
+}
+
+int rb_tree_insert(rb_tree *t, int k)
+{
+     rb_node *n = alloc1(t, k);
+     if (!n) return 0;
+     insert_node(t, n);
+     return 1;
+}
+
+static int check_node(rb_node *n, int *nblack, rb_node *nil)
+{
+     int nbl, nbr;
+     if (n == nil) { *nblack = 0; return 1; }
+     if (n->r != nil && n->r->p != n) return 0;
+     if (n->l != nil && n->l->p != n) return 0;
+     if (n->c == RED) {
+	  if (n->r != nil && n->r->c == RED) return 0;
+	  if (n->l != nil && n->l->c == RED) return 0;
+     }
+     if (!(check_node(n->r, &nbl, nil) && check_node(n->l, &nbr, nil))) 
+	  return 0;
+     if (nbl != nbr) return 0;
+     *nblack = nbl + n->c == BLACK;
+     return 1;
+}
+int rb_tree_check(rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     int nblack;
+     if (nil->c != BLACK) return 0;
+     if (t->root == nil) return 1;
+     if (t->root->c != BLACK) return 0;
+     return check_node(t->root, &nblack, nil);
+}
+
+rb_node *rb_tree_find(rb_tree *t, int k)
+{
+     rb_node *nil = &t->nil;
+     rb_compare compare = t->compare;
+     void *c_data = t->c_data;
+     rb_node *p = t->root;
+     while (p != nil) {
+	  int comp = compare(k, p->k, c_data);
+	  if (!comp) return p;
+	  p = comp <= 0 ? p->l : p->r;
+     }
+     return NULL;
+}
+
+rb_node *rb_tree_min(rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_node *n = t->root;
+     while (n != nil && n->l != nil)
+	  n = n->l;
+     return(n == nil ? NULL : n);
+}
+
+rb_node *rb_tree_max(rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_node *n = t->root;
+     while (n != nil && n->r != nil)
+	  n = n->r;
+     return(n == nil ? NULL : n);
+}
+
+rb_node *rb_tree_succ(rb_tree *t, rb_node *n)
+{
+     rb_node *nil = &t->nil;
+     if (n->r == nil) {
+	  rb_node *prev;
+	  do {
+	       prev = n;
+	       n = n->p;
+	  } while (prev == n->r && n != nil);
+	  return n == nil ? NULL : n;
+     }
+     else {
+	  n = n->r;
+	  while (n->l != nil)
+	       n = n->l;
+	  return n;
+     }
+}
+
+rb_node *rb_tree_pred(rb_tree *t, rb_node *n)
+{
+     rb_node *nil = &t->nil;
+     if (n->l == nil) {
+	  rb_node *prev;
+	  do {
+	       prev = n;
+	       n = n->p;
+	  } while (prev == n->l && n != nil);
+	  return n == nil ? NULL : n;
+     }
+     else {
+	  n = n->l;
+	  while (n->r != nil)
+	       n = n->r;
+	  return n;
+     }
+}
+
+rb_node *rb_tree_remove(rb_tree *t, rb_node *n)
+{
+     rb_node *nil = &t->nil;
+     rb_node *m;
+     if (n->l != nil && n->r != nil) {
+	  rb_node *lmax = n->l;
+	  while (lmax->r != nil) lmax = lmax->r;
+	  n->k = lmax->k;
+	  n = lmax;
+     }
+     m = n->l != nil? n->l : n->r;
+     if (n->p != nil) {
+	  if (n->p->r == n) n->p->r = m;
+	  else n->p->l = m;
+     }
+     else
+	  t->root = m;
+     m->p = n->p;
+     if (n->c == BLACK) {
+	  if (m->c == RED)
+	       m->c = BLACK;
+	  else {
+	  deleteblack:
+	       if (m->p != nil) {
+		    rb_node *s = m == m->p->l ? m->p->r : m->p->l;
+		    if (s->c == RED) {
+			 m->p->c = RED;
+			 s->c = BLACK;
+			 if (m == m->p->l) rotate_left(m->p, t);
+			 else rotate_right(m->p, t);
+			 s = m == m->p->l ? m->p->r : m->p->l;
+		    }
+		    if (m->p->c == BLACK && s->c == BLACK
+			&& s->l->c == BLACK && s->r->c == BLACK) {
+			 if (s != nil) s->c = RED;
+			 m = m->p;
+			 goto deleteblack;
+		    }
+		    else if (m->p->c == RED && s->c == BLACK &&
+			     s->l->c == BLACK && s->r->c == BLACK) {
+			 if (s != nil) s->c = RED;
+			 m->p->c = BLACK;
+		    }
+		    else {
+			 if (m == m->p->l && s->c == BLACK &&
+			     s->l->c == RED && s->r->c == BLACK) {
+			      s->c = RED;
+			      s->l->c = BLACK;
+			      rotate_right(s, t);
+			      s = m == m->p->l ? m->p->r : m->p->l;
+			 }
+			 else if (m == m->p->r && s->c == BLACK &&
+				  s->r->c == RED && s->l->c == BLACK) {
+			      s->c = RED;
+			      s->r->c = BLACK;
+			      rotate_left(s, t);
+			      s = m == m->p->l ? m->p->r : m->p->l;
+			 }
+			 s->c = m->p->c;
+			 m->p->c = BLACK;
+			 if (m == m->p->l) {
+			      s->r->c = BLACK;
+			      rotate_left(m->p, t);
+			 }
+			 else {
+			      s->l->c = BLACK;
+			      rotate_right(m->p, t);
+			 }
+		    }
+	       }
+	  }
+     }
+     return n; /* the node that was deleted may be different from initial n */
+}
+
+rb_node *rb_tree_resort(rb_tree *t, rb_node *n)
+{
+     int k = n->k;
+     n = rb_tree_remove(t, n);
+     n->p = n->l = n->r = &t->nil;
+     n->k = k; /* n may have changed during remove */
+     insert_node(t, n);
+     return n;
+}

cdirect/redblack.h

+#ifndef REDBLACK_H
+#define REDBLACK_H
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+typedef enum { RED, BLACK } rb_color;
+typedef struct rb_node_s {
+     struct rb_node_s *p, *r, *l; /* parent, right, left */
+     int k; /* key/data ... for DIRECT, an index into our hyperrect array */
+     rb_color c;
+} rb_node;
+
+typedef int (*rb_compare)(int k1, int k2, void *c_data);
+
+typedef struct {
+     rb_compare compare; void *c_data;
+     rb_node *root;
+     int N; /* number of nodes */
+
+     /* in our application, we can optimize memory allocation because
+	we never delete two nodes in a row (we always add a node after
+	deleting)... or rather, we never delete but the value of
+        the key sometimes changes.  ... this means we can just
+	allocate a linear, exponentially growing stack (nodes) of
+	nodes, and don't have to worry about holes in the stack */
+     rb_node *nodes; /* allocated data of nodes, in some order */
+     int Nalloc; /* number of allocated nodes */
+     rb_node nil; /* explicit node for NULL nodes, for convenience */
+} rb_tree;
+
+extern int rb_tree_init(rb_tree *t, rb_compare compare, void *c_data);
+extern void rb_tree_destroy(rb_tree *t);
+extern int rb_tree_insert(rb_tree *t, int k);
+extern int rb_tree_check(rb_tree *t);
+extern rb_node *rb_tree_find(rb_tree *t, int k);
+extern rb_node *rb_tree_resort(rb_tree *t, rb_node *n);
+extern rb_node *rb_tree_min(rb_tree *t);
+extern rb_node *rb_tree_max(rb_tree *t);
+extern rb_node *rb_tree_succ(rb_tree *t, rb_node *n);
+extern rb_node *rb_tree_pred(rb_tree *t, rb_node *n);
+
+/* To change a key, use rb_tree_find+resort.  Removing a node
+   currently wastes memory unless you change the allocation scheme
+   in redblack.c */
+extern rb_node *rb_tree_remove(rb_tree *t, rb_node *n);
+
+#ifdef __cplusplus
+}  /* extern "C" */
+#endif /* __cplusplus */
+
+#endif

cdirect/redblack_test.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include "redblack.h"
+
+static int comp(int k1, int k2, void *dummy)
+{
+     (void) dummy;
+     return k1 - k2;
+}
+
+int main(int argc, char **argv)
+{
+     int N, M;
+     int *k;
+     rb_tree t;
+     rb_node *n;
+     int i, j;
+
+     if (argc != 2) {
+	  fprintf(stderr, "Usage: redblack_test Ntest\n");
+	  return 1;
+     }
+
+     N = atoi(argv[1]);
+     k = (int *) malloc(N * sizeof(int));
+     if (!rb_tree_init(&t, comp, NULL)) {
+	  fprintf(stderr, "error in rb_tree_init\n");
+	  return 1;
+     }
+
+     srand((unsigned) time(NULL));
+     for (i = 0; i < N; ++i) {
+	  if (!rb_tree_insert(&t, k[i] = rand() % N)) {
+	       fprintf(stderr, "error in rb_tree_insert\n");
+	       return 1;
+	  }
+	  if (!rb_tree_check(&t)) {
+	       fprintf(stderr, "rb_tree_check_failed after insert!\n");
+	       return 1;
+	  }
+     }
+     
+     for (i = 0; i < N; ++i)
+	  if (!rb_tree_find(&t, k[i])) {
+	       fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+	       return 1;
+	  }
+     
+     n = rb_tree_min(&t);
+     for (i = 0; i < N; ++i) {
+	  if (!n) {
+	       fprintf(stderr, "not enough successors %d\n!", i);
+	       return 1;
+	  }
+	  printf("%d: %d\n", i, n->k);
+	  n = rb_tree_succ(&t, n);
+     }
+     if (n) {
+	  fprintf(stderr, "too many successors!\n");
+	  return 1;
+     }
+     
+     n = rb_tree_max(&t);
+     for (i = 0; i < N; ++i) {
+	  if (!n) {
+	       fprintf(stderr, "not enough predecessors %d\n!", i);
+	       return 1;
+	  }
+	  printf("%d: %d\n", i, n->k);
+	  n = rb_tree_pred(&t, n);
+     }
+     if (n) {
+	  fprintf(stderr, "too many predecessors!\n");
+	  return 1;
+     }
+     
+     for (M = N; M > 0; --M) {
+	  int knew;
+	  j = rand() % M;
+	  for (i = 0; i < N; ++i)
+	       if (k[i] >= 0)
+		    if (j-- == 0)
+			 break;
+	  if (i >= N) abort();
+	  if (!(n = rb_tree_find(&t, k[i]))) {
+               fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+               return 1;
+          }
+	  n->k = knew;
+	  if (!rb_tree_resort(&t, n)) {
+	       fprintf(stderr, "error in rb_tree_resort\n");
+	       return 1;
+	  }
+	  if (!rb_tree_check(&t)) {
+	       fprintf(stderr, "rb_tree_check_failed after change!\n");
+	       return 1;
+	  }
+	  k[i] = -1 - knew;
+     }
+
+     for (i = 0; i < N; ++i)
+	  k[i] = -1 - k[i];
+     
+     for (M = N; M > 0; --M) {
+	  j = rand() % M;
+	  for (i = 0; i < N; ++i)
+	       if (k[i] >= 0)
+		    if (j-- == 0)
+			 break;
+	  if (i >= N) abort();
+	  if (!(n = rb_tree_find(&t, k[i]))) {
+	       fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+	       return 1;
+	  }
+	  rb_tree_remove(&t, n);
+	  if (!rb_tree_check(&t)) {
+	       fprintf(stderr, "rb_tree_check_failed after remove!\n");
+	       return 1;
+	  }
+	  k[i] = -1 - k[i];
+     }
+     
+     rb_tree_destroy(&t);
+     free(k);
+     return 0;
+}