src/libpspp/bt.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 2007 Free Software Foundation, Inc.
   3
   4    This program is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU General Public License as
   6    published by the Free Software Foundation; either version 2 of the
   7    License, or (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful, but
  10    WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  17    02110-1301, USA. */
  18
  19 /* Balanced tree (BT) data structure.
  20
  21    The client should not need to be aware of the form of
  22    balancing applied to the balanced tree, as its operation is
  23    fully encapsulated.  The current implementation is a scapegoat
  24    tree.  Scapegoat trees have the advantage over many other
  25    forms of balanced trees that they do not store any additional
  26    information in each node; thus, they are slightly more
  27    space-efficient than, say, red-black or AVL trees.  Compared
  28    to splay trees, scapegoat trees provide guaranteed logarithmic
  29    worst-case search time and do not restructure the tree during
  30    a search.
  31
  32    For information on scapegoat trees, see Galperin and Rivest,
  33    "Scapegoat Trees", Proc. 4th ACM-SIAM Symposium on Discrete
  34    Algorithms, or <http://en.wikipedia.org/wiki/Scapegoat_tree>,
  35    which includes source code and links to other resources, such
  36    as the Galperin and Rivest paper.
  37
  38    One potentially tricky part of scapegoat tree design is the
  39    choice of alpha, which is a real constant that must be greater
  40    than 0.5 and less than 1.  We must be able to efficiently
  41    calculate h_alpha = floor(log(n)/log(1/alpha)) for integer n >
  42    0.  One way to do so would be to maintain a table relating
  43    h_alpha to the minimum value of n that yields that h_alpha.
  44    Then, we can track the value of h_alpha(n) in the number of
  45    nodes in the tree n as nodes are inserted and deleted.
  46
  47    This implementation uses a different approach.  We choose
  48    alpha = sqrt(2)/2 = 1/sqrt(2) ~= .707.  Then, computing
  49    h_alpha is a matter of computing a logarithm in base sqrt(2).
  50    This is easy: we simply compute twice the base-2 logarithm,
  51    then adjust upward by 1 if necessary.  See calculate_h_alpha
  52    for details. */
  53
  54 /* These library routines have no external dependencies other
  55    than the standard C library.
  56
  57    If you add routines in this file, please add a corresponding
  58    test to bt-test.c.  This test program should achieve 100%
  59    coverage of lines and branches in this code, as reported by
  60    "gcov -b". */
  61
  62 #ifdef HAVE_CONFIG_H
  63 #include <config.h>
  64 #endif
  65
  66 #include <libpspp/bt.h>
  67
  68 #include <stdbool.h>
  69 #include <stdint.h>
  70
  71 static void rebalance_subtree (struct bt *, struct bt_node *, size_t);
  72
  73 static struct bt_node **down_link (struct bt *, struct bt_node *);
  74 static inline struct bt_node *sibling (struct bt_node *p);
  75 static size_t count_nodes_in_subtree (const struct bt_node *);
  76
  77 static inline int floor_log2 (size_t);
  78 static inline int calculate_h_alpha (size_t);
  79
  80 /* Initializes BT as an empty BT that uses the given COMPARE
  81    function, passing in AUX as auxiliary data. */
  82 void
  83 bt_init (struct bt *bt,
  84          bt_compare_func *compare,
  85          const void *aux)
  86 {
  87   bt->root = NULL;
  88   bt->compare = compare;
  89   bt->aux = aux;
  90   bt->size = 0;
  91   bt->max_size = 0;
  92 }
  93
  94 /* Inserts the given NODE into BT.
  95    Returns a null pointer if successful.
  96    Returns the existing node already in BT equal to NODE, on
  97    failure. */
  98 struct bt_node *
  99 bt_insert (struct bt *bt, struct bt_node *node)
 100 {
 101   int depth = 0;
 102
 103   node->down[0] = NULL;
 104   node->down[1] = NULL;
 105
 106   if (bt->root == NULL)
 107     {
 108       bt->root = node;
 109       node->up = NULL;
 110     }
 111   else
 112     {
 113       struct bt_node *p = bt->root;
 114       for (;;)
 115         {
 116           int cmp, dir;
 117
 118           cmp = bt->compare (node, p, bt->aux);
 119           if (cmp == 0)
 120             return p;
 121           depth++;
 122
 123           dir = cmp > 0;
 124           if (p->down[dir] == NULL)
 125             {
 126               p->down[dir] = node;
 127               node->up = p;
 128               break;
 129             }
 130           p = p->down[dir];
 131         }
 132     }
 133
 134   bt->size++;
 135   if (bt->size > bt->max_size)
 136     bt->max_size = bt->size;
 137
 138   if (depth > calculate_h_alpha (bt->size))
 139     {
 140       /* We use the "alternative" method of finding a scapegoat
 141          node described by Galperin and Rivest. */
 142       struct bt_node *s = node;
 143       size_t size = 1;
 144       int i;
 145
 146       for (i = 1; ; i++)
 147         if (i < depth)
 148           {
 149             size += 1 + count_nodes_in_subtree (sibling (s));
 150             s = s->up;
 151             if (i > calculate_h_alpha (size))
 152               {
 153                 rebalance_subtree (bt, s, size);
 154                 break;
 155               }
 156           }
 157         else
 158           {
 159             rebalance_subtree (bt, bt->root, bt->size);
 160             bt->max_size = bt->size;
 161             break;
 162           }
 163     }
 164
 165   return NULL;
 166 }
 167
 168 /* Deletes P from BT. */
 169 void
 170 bt_delete (struct bt *bt, struct bt_node *p)
 171 {
 172   struct bt_node **q = down_link (bt, p);
 173   struct bt_node *r = p->down[1];
 174   if (r == NULL)
 175     {
 176       *q = p->down[0];
 177       if (*q)
 178         (*q)->up = p->up;
 179     }
 180   else if (r->down[0] == NULL)
 181     {
 182       r->down[0] = p->down[0];
 183       *q = r;
 184       r->up = p->up;
 185       if (r->down[0] != NULL)
 186         r->down[0]->up = r;
 187     }
 188   else
 189     {
 190       struct bt_node *s = r->down[0];
 191       while (s->down[0] != NULL)
 192         s = s->down[0];
 193       r = s->up;
 194       r->down[0] = s->down[1];
 195       s->down[0] = p->down[0];
 196       s->down[1] = p->down[1];
 197       *q = s;
 198       if (s->down[0] != NULL)
 199         s->down[0]->up = s;
 200       s->down[1]->up = s;
 201       s->up = p->up;
 202       if (r->down[0] != NULL)
 203         r->down[0]->up = r;
 204     }
 205   bt->size--;
 206
 207   /* We approximate .707 as .75 here.  This is conservative: it
 208      will cause us to do a little more rebalancing than strictly
 209      necessary to maintain the scapegoat tree's height
 210      invariant. */
 211   if (bt->size < bt->max_size * 3 / 4 && bt->size > 0)
 212     {
 213       rebalance_subtree (bt, bt->root, bt->size);
 214       bt->max_size = bt->size;
 215     }
 216 }
 217
 218 /* Returns the node with minimum value in BT, or a null pointer
 219    if BT is empty. */
 220 struct bt_node *
 221 bt_first (const struct bt *bt)
 222 {
 223   struct bt_node *p = bt->root;
 224   if (p != NULL)
 225     while (p->down[0] != NULL)
 226       p = p->down[0];
 227   return p;
 228 }
 229
 230 /* Returns the node with maximum value in BT, or a null pointer
 231    if BT is empty. */
 232 struct bt_node *
 233 bt_last (const struct bt *bt)
 234 {
 235   struct bt_node *p = bt->root;
 236   if (p != NULL)
 237     while (p->down[1] != NULL)
 238       p = p->down[1];
 239   return p;
 240 }
 241
 242 /* Searches BT for a node equal to TARGET.
 243    Returns the node if found, or a null pointer otherwise. */
 244 struct bt_node *
 245 bt_find (const struct bt *bt, const struct bt_node *target)
 246 {
 247   const struct bt_node *p;
 248   int cmp;
 249
 250   for (p = bt->root; p != NULL; p = p->down[cmp > 0])
 251     {
 252       cmp = bt->compare (target, p, bt->aux);
 253       if (cmp == 0)
 254         return (struct bt_node *) p;
 255     }
 256
 257   return NULL;
 258 }
 259
 260 /* Searches BT for, and returns, the first node in in-order whose
 261    value is greater than or equal to TARGET.  Returns a null
 262    pointer if all nodes in BT are less than TARGET.
 263
 264    Another way to look at the return value is that it is the node
 265    that would be returned by "bt_next (BT, TARGET)" if TARGET
 266    were inserted in BT (assuming that TARGET would not be a
 267    duplicate). */
 268 struct bt_node *
 269 bt_find_ge (const struct bt *bt, const struct bt_node *target)
 270 {
 271   const struct bt_node *p = bt->root;
 272   const struct bt_node *q = NULL;
 273   while (p != NULL)
 274     {
 275       int cmp = bt->compare (target, p, bt->aux);
 276       if (cmp > 0)
 277         p = p->down[1];
 278       else
 279         {
 280           q = p;
 281           if (cmp < 0)
 282             p = p->down[0];
 283           else
 284             break;
 285         }
 286     }
 287   return (struct bt_node *) q;
 288 }
 289
 290 /* Searches BT for, and returns, the last node in in-order whose
 291    value is less than or equal to TARGET, which should not be in
 292    BT.  Returns a null pointer if all nodes in BT are greater
 293    than TARGET.
 294
 295    Another way to look at the return value is that it is the node
 296    that would be returned by "bt_prev (BT, TARGET)" if TARGET
 297    were inserted in BT (assuming that TARGET would not be a
 298    duplicate). */
 299 struct bt_node *
 300 bt_find_le (const struct bt *bt, const struct bt_node *target)
 301 {
 302   const struct bt_node *p = bt->root;
 303   const struct bt_node *q = NULL;
 304   while (p != NULL)
 305     {
 306       int cmp = bt->compare (target, p, bt->aux);
 307       if (cmp < 0)
 308         p = p->down[0];
 309       else
 310         {
 311           q = p;
 312           if (cmp > 0)
 313             p = p->down[1];
 314           else
 315             break;
 316         }
 317     }
 318   return (struct bt_node *) q;
 319 }
 320
 321 /* Returns the node in BT following P in in-order.
 322    If P is null, returns the minimum node in BT.
 323    Returns a null pointer if P is the maximum node in BT or if P
 324    is null and BT is empty. */
 325 struct bt_node *
 326 bt_next (const struct bt *bt, const struct bt_node *p)
 327 {
 328   if (p == NULL)
 329     return bt_first (bt);
 330   else if (p->down[1] == NULL)
 331     {
 332       struct bt_node *q;
 333       for (q = p->up; ; p = q, q = q->up)
 334         if (q == NULL || p == q->down[0])
 335           return q;
 336     }
 337   else
 338     {
 339       p = p->down[1];
 340       while (p->down[0] != NULL)
 341         p = p->down[0];
 342       return (struct bt_node *) p;
 343     }
 344 }
 345
 346 /* Returns the node in BT preceding P in in-order.
 347    If P is null, returns the maximum node in BT.
 348    Returns a null pointer if P is the minimum node in BT or if P
 349    is null and BT is empty. */
 350 struct bt_node *
 351 bt_prev (const struct bt *bt, const struct bt_node *p)
 352 {
 353   if (p == NULL)
 354     return bt_last (bt);
 355   else if (p->down[0] == NULL)
 356     {
 357       struct bt_node *q;
 358       for (q = p->up; ; p = q, q = q->up)
 359         if (q == NULL || p == q->down[1])
 360           return q;
 361     }
 362   else
 363     {
 364       p = p->down[0];
 365       while (p->down[1] != NULL)
 366         p = p->down[1];
 367       return (struct bt_node *) p;
 368     }
 369 }
 370
 371 /* Moves P around in BT to compensate for its key having
 372    changed.  Returns a null pointer if successful.  If P's new
 373    value is equal to that of some other node in BT, returns the
 374    other node after removing P from BT.
 375
 376    This function is an optimization only if it is likely that P
 377    can actually retain its relative position in BT, e.g. its key
 378    has only been adjusted slightly.  Otherwise, it is more
 379    efficient to simply remove P from BT, change its key, and
 380    re-insert P.
 381
 382    It is not safe to update more than one node's key, then to
 383    call this function for each node.  Instead, update a single
 384    node's key, call this function, update another node's key, and
 385    so on.  Alternatively, remove all affected nodes from the
 386    tree, update their keys, then re-insert all of them. */
 387 struct bt_node *
 388 bt_changed (struct bt *bt, struct bt_node *p)
 389 {
 390   struct bt_node *prev = bt_prev (bt, p);
 391   struct bt_node *next = bt_next (bt, p);
 392
 393   if ((prev != NULL && bt->compare (prev, p, bt->aux) >= 0)
 394       || (next != NULL && bt->compare (p, next, bt->aux) >= 0))
 395     {
 396       bt_delete (bt, p);
 397       return bt_insert (bt, p);
 398     }
 399   return NULL;
 400  }
 401
 402 /* BT nodes may be moved around in memory as necessary, e.g. as
 403    the result of an realloc operation on a block that contains a
 404    node.  Once this is done, call this function passing node P
 405    that was moved and its BT before attempting any other
 406    operation on BT.
 407
 408    It is not safe to move more than one node, then to call this
 409    function for each node.  Instead, move a single node, call
 410    this function, move another node, and so on.  Alternatively,
 411    remove all affected nodes from the tree, move them, then
 412    re-insert all of them. */
 413 void
 414 bt_moved (struct bt *bt, struct bt_node *p)
 415 {
 416   if (p->up != NULL)
 417     {
 418       int d = p->up->down[0] == NULL || bt->compare (p, p->up, bt->aux) > 0;
 419       p->up->down[d] = p;
 420     }
 421   else
 422     bt->root = p;
 423   if (p->down[0] != NULL)
 424     p->down[0]->up = p;
 425   if (p->down[1] != NULL)
 426     p->down[1]->up = p;
 427 }
 428 \f
 429 /* Tree rebalancing.
 430
 431    This algorithm is from Q. F. Stout and B. L. Warren, "Tree
 432    Rebalancing in Optimal Time and Space", CACM 29(1986):9,
 433    pp. 902-908.  It uses O(N) time and O(1) space to rebalance a
 434    subtree that contains N nodes. */
 435
 436 static void tree_to_vine (struct bt_node **);
 437 static void vine_to_tree (struct bt_node **, size_t count);
 438
 439 /* Rebalances the subtree in BT rooted at SUBTREE, which contains
 440    exactly COUNT nodes. */
 441 static void
 442 rebalance_subtree (struct bt *bt, struct bt_node *subtree, size_t count)
 443 {
 444   struct bt_node *up = subtree->up;
 445   struct bt_node **q = down_link (bt, subtree);
 446   tree_to_vine (q);
 447   vine_to_tree (q, count);
 448   (*q)->up = up;
 449 }
 450
 451 /* Converts the subtree rooted at *Q into a vine (a binary search
 452    tree in which all the right links are null), and updates *Q to
 453    point to the new root of the subtree. */
 454 static void
 455 tree_to_vine (struct bt_node **q)
 456 {
 457   struct bt_node *p = *q;
 458   while (p != NULL)
 459     if (p->down[1] == NULL)
 460       {
 461         q = &p->down[0];
 462         p = *q;
 463       }
 464     else
 465       {
 466         struct bt_node *r = p->down[1];
 467         p->down[1] = r->down[0];
 468         r->down[0] = p;
 469         p = r;
 470         *q = r;
 471       }
 472 }
 473
 474 /* Performs a compression transformation COUNT times, starting at
 475    *Q, and updates *Q to point to the new root of the subtree. */
 476 static void
 477 compress (struct bt_node **q, unsigned long count)
 478 {
 479   while (count--)
 480     {
 481       struct bt_node *red = *q;
 482       struct bt_node *black = red->down[0];
 483
 484       *q = black;
 485       red->down[0] = black->down[1];
 486       black->down[1] = red;
 487       red->up = black;
 488       if (red->down[0] != NULL)
 489         red->down[0]->up = red;
 490       q = &black->down[0];
 491     }
 492 }
 493
 494 /* Converts the vine rooted at *Q, which contains exactly COUNT
 495    nodes, into a balanced tree, and updates *Q to point to the
 496    new root of the balanced tree. */
 497 static void
 498 vine_to_tree (struct bt_node **q, size_t count)
 499 {
 500   size_t leaf_nodes = count + 1 - ((size_t) 1 << floor_log2 (count + 1));
 501   size_t vine_nodes = count - leaf_nodes;
 502
 503   compress (q, leaf_nodes);
 504   while (vine_nodes > 1)
 505     {
 506       vine_nodes /= 2;
 507       compress (q, vine_nodes);
 508     }
 509   while ((*q)->down[0] != NULL)
 510     {
 511       (*q)->down[0]->up = *q;
 512       q = &(*q)->down[0];
 513     }
 514 }
 515 \f
 516 /* Other binary tree helper functions. */
 517
 518 /* Returns the address of the pointer that points down to P
 519    within BT. */
 520 static struct bt_node **
 521 down_link (struct bt *bt, struct bt_node *p)
 522 {
 523   struct bt_node *q = p->up;
 524   return q != NULL ? &q->down[q->down[0] != p] : &bt->root;
 525 }
 526
 527 /* Returns node P's sibling; that is, the other child of its
 528    parent.  P must not be the root. */
 529 static inline struct bt_node *
 530 sibling (struct bt_node *p)
 531 {
 532   struct bt_node *q = p->up;
 533   return q->down[q->down[0] == p];
 534 }
 535
 536 /* Returns the number of nodes in the given SUBTREE. */
 537 static size_t
 538 count_nodes_in_subtree (const struct bt_node *subtree)
 539 {
 540   /* This is an in-order traversal modified to iterate only the
 541      nodes in SUBTREE. */
 542   size_t count = 0;
 543   if (subtree != NULL)
 544     {
 545       const struct bt_node *p = subtree;
 546       while (p->down[0] != NULL)
 547         p = p->down[0];
 548       for (;;)
 549         {
 550           count++;
 551           if (p->down[1] != NULL)
 552             {
 553               p = p->down[1];
 554               while (p->down[0] != NULL)
 555                 p = p->down[0];
 556             }
 557           else
 558             {
 559               for (;;)
 560                 {
 561                   const struct bt_node *q;
 562                   if (p == subtree)
 563                     goto done;
 564                   q = p;
 565                   p = p->up;
 566                   if (p->down[0] == q)
 567                     break;
 568                 }
 569             }
 570         }
 571     }
 572  done:
 573   return count;
 574 }
 575 \f
 576 /* Arithmetic. */
 577
 578 /* Returns the number of high-order 0-bits in X.
 579    Undefined if X is zero. */
 580 static inline int
 581 count_leading_zeros (size_t x)
 582 {
 583 #if __GNUC__ >= 4
 584 #if SIZE_MAX > ULONG_MAX
 585   return __builtin_clzll (x);
 586 #elif SIZE_MAX > UINT_MAX
 587   return __builtin_clzl (x);
 588 #else
 589   return __builtin_clz (x);
 590 #endif
 591 #else
 592   /* This algorithm is from _Hacker's Delight_ section 5.3. */
 593   size_t y;
 594   int n;
 595
 596 #define COUNT_STEP(BITS)                        \
 597         y = x >> BITS;                          \
 598         if (y != 0)                             \
 599           {                                     \
 600             n -= BITS;                          \
 601             x = y;                              \
 602           }
 603
 604   n = sizeof (size_t) * CHAR_BIT;
 605 #if SIZE_MAX >> 31 >> 31 >> 2
 606   COUNT_STEP (64);
 607 #endif
 608 #if SIZE_MAX >> 31 >> 1
 609   COUNT_STEP (32);
 610 #endif
 611   COUNT_STEP (16);
 612   COUNT_STEP (8);
 613   COUNT_STEP (4);
 614   COUNT_STEP (2);
 615   y = x >> 1;
 616   return y != 0 ? n - 2 : n - x;
 617 #endif
 618 }
 619
 620 /* Returns floor(log2(x)).
 621    Undefined if X is zero. */
 622 static inline int
 623 floor_log2 (size_t x)
 624 {
 625   return sizeof (size_t) * CHAR_BIT - 1 - count_leading_zeros (x);
 626 }
 627
 628 /* Returns floor(pow(sqrt(2), x * 2 + 1)).
 629    Defined for X from 0 up to the number of bits in size_t minus
 630    1. */
 631 static inline size_t
 632 pow_sqrt2 (int x)
 633 {
 634   /* These constants are sqrt(2) multiplied by 2**63 or 2**31,
 635      respectively, and then rounded to nearest. */
 636 #if SIZE_MAX >> 31 >> 1
 637   return (UINT64_C(0xb504f333f9de6484) >> (63 - x)) + 1;
 638 #else
 639   return (0xb504f334 >> (31 - x)) + 1;
 640 #endif
 641 }
 642
 643 /* Returns floor(log(n)/log(sqrt(2))).
 644    Undefined if N is 0. */
 645 static inline int
 646 calculate_h_alpha (size_t n)
 647 {
 648   int log2 = floor_log2 (n);
 649
 650   /* The correct answer is either 2 * log2 or one more.  So we
 651      see if n >= pow(sqrt(2), 2 * log2 + 1) and if so, add 1. */
 652   return (2 * log2) + (n >= pow_sqrt2 (log2));
 653 }
 654