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
bt_init(struct bt * bt,bt_compare_func * compare,const void * aux)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 *
bt_insert(struct bt * bt,struct bt_node * 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
bt_delete(struct bt * bt,struct bt_node * p)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 *
bt_first(const struct bt * bt)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 *
bt_last(const struct bt * bt)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 *
bt_find(const struct bt * bt,const struct bt_node * target)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 *
bt_find_ge(const struct bt * bt,const struct bt_node * target)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 *
bt_find_le(const struct bt * bt,const struct bt_node * target)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 *
bt_next(const struct bt * bt,const struct bt_node * p)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 *
bt_prev(const struct bt * bt,const struct bt_node * p)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 *
bt_changed(struct bt * bt,struct bt_node * p)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
bt_moved(struct bt * bt,struct bt_node * p)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
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
rebalance_subtree(struct bt * bt,struct bt_node * subtree,size_t count)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
tree_to_vine(struct bt_node ** q)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
compress(struct bt_node ** q,unsigned long count)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
vine_to_tree(struct bt_node ** q,size_t count)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
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 **
down_link(struct bt * bt,struct bt_node * p)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 *
sibling(struct bt_node * p)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
count_nodes_in_subtree(const struct bt_node * subtree)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
575 /* Arithmetic. */
576
577 /* Returns the number of high-order 0-bits in X.
578 Undefined if X is zero. */
579 static inline int
count_leading_zeros(size_t x)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
floor_log2(size_t x)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
pow_sqrt2(int x)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
calculate_h_alpha(size_t n)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
654