gui: Fix behavior of print preview in psppire-output-window.
[pspp] / 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