src/libpspp/bt.c

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