f4a1ef6d953843591f9f77881160981e7f69bdd6
[pspp-builds.git] / src / libpspp / array.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 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 /* Copyright (C) 2001 Free Software Foundation, Inc.
20
21    This file is part of the GNU ISO C++ Library.  This library is free
22    software; you can redistribute it and/or modify it under the
23    terms of the GNU General Public License as published by the
24    Free Software Foundation; either version 2, or (at your option)
25    any later version.
26
27    This library is distributed in the hope that it will be useful,
28    but WITHOUT ANY WARRANTY; without even the implied warranty of
29    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30    GNU General Public License for more details.
31
32    You should have received a copy of the GNU General Public License along
33    with this library; see the file COPYING.  If not, write to the Free
34    Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
35    USA.
36
37    As a special exception, you may use this file as part of a free software
38    library without restriction.  Specifically, if other files instantiate
39    templates or use macros or inline functions from this file, or you compile
40    this file and link it with other files to produce an executable, this
41    file does not by itself cause the resulting executable to be covered by
42    the GNU General Public License.  This exception does not however
43    invalidate any other reasons why the executable file might be covered by
44    the GNU General Public License. */
45
46 /*
47  *
48  * Copyright (c) 1994
49  * Hewlett-Packard Company
50  *
51  * Permission to use, copy, modify, distribute and sell this software
52  * and its documentation for any purpose is hereby granted without fee,
53  * provided that the above copyright notice appear in all copies and
54  * that both that copyright notice and this permission notice appear
55  * in supporting documentation.  Hewlett-Packard Company makes no
56  * representations about the suitability of this software for any
57  * purpose.  It is provided "as is" without express or implied warranty.
58  *
59  *
60  * Copyright (c) 1996
61  * Silicon Graphics Computer Systems, Inc.
62  *
63  * Permission to use, copy, modify, distribute and sell this software
64  * and its documentation for any purpose is hereby granted without fee,
65  * provided that the above copyright notice appear in all copies and
66  * that both that copyright notice and this permission notice appear
67  * in supporting documentation.  Silicon Graphics makes no
68  * representations about the suitability of this software for any
69  * purpose.  It is provided "as is" without express or implied warranty.
70  */
71
72 /* Copyright (C) 1991, 1992, 1996, 1997, 1999 Free Software Foundation, Inc.
73    This file is part of the GNU C Library.
74    Written by Douglas C. Schmidt (schmidt@ics.uci.edu).
75
76    The GNU C Library is free software; you can redistribute it and/or
77    modify it under the terms of the GNU Lesser General Public
78    License as published by the Free Software Foundation; either
79    version 2.1 of the License, or (at your option) any later version.
80
81    The GNU C Library is distributed in the hope that it will be useful,
82    but WITHOUT ANY WARRANTY; without even the implied warranty of
83    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
84    Lesser General Public License for more details.
85
86    You should have received a copy of the GNU Lesser General Public
87    License along with the GNU C Library; if not, write to the Free
88    Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
89    02110-1301 USA.  */
90
91 #include <config.h>
92 #include "array.h"
93 #include <limits.h>
94 #include <stdlib.h>
95 #include <string.h>
96 #include "alloc.h"
97 #include <libpspp/assertion.h>
98
99 #include "message.h"
100
101 #include "minmax.h"
102 \f
103 /* Finds an element in ARRAY, which contains COUNT elements of
104    SIZE bytes each, using COMPARE for comparisons.  Returns the
105    first element in ARRAY that matches TARGET, or a null pointer
106    on failure.  AUX is passed to each comparison as auxiliary
107    data. */
108 void *
109 find (const void *array, size_t count, size_t size,
110       const void *target,
111       algo_compare_func *compare, const void *aux)
112 {
113   const char *element = array;
114
115   while (count-- > 0)
116     {
117       if (compare (target, element, aux) == 0)
118         return (void *) element;
119
120       element += size;
121     }
122
123   return NULL;
124 }
125
126 /* Counts and return the number of elements in ARRAY, which
127    contains COUNT elements of SIZE bytes each, which are equal to
128    ELEMENT as compared with COMPARE.  AUX is passed as auxiliary
129    data to COMPARE. */
130 size_t
131 count_equal (const void *array, size_t count, size_t size,
132              const void *element,
133              algo_compare_func *compare, const void *aux)
134 {
135   const char *first = array;
136   size_t equal_cnt = 0;
137
138   while (count-- > 0)
139     {
140       if (compare (element, first, aux) == 0)
141         equal_cnt++;
142
143       first += size;
144     }
145
146   return equal_cnt;
147 }
148
149 /* Counts and return the number of elements in ARRAY, which
150    contains COUNT elements of SIZE bytes each, for which
151    PREDICATE returns true.  AUX is passed as auxiliary data to
152    PREDICATE. */
153 size_t
154 count_if (const void *array, size_t count, size_t size,
155           algo_predicate_func *predicate, const void *aux)
156 {
157   const char *first = array;
158   size_t true_cnt = 0;
159
160   while (count-- > 0)
161     {
162       if (predicate (first, aux) != 0)
163         true_cnt++;
164
165       first += size;
166     }
167
168   return true_cnt;
169 }
170 \f
171 /* Byte-wise swap two items of size SIZE. */
172 #define SWAP(a, b, size)                        \
173   do                                            \
174     {                                           \
175       register size_t __size = (size);          \
176       register char *__a = (a), *__b = (b);     \
177       do                                        \
178         {                                       \
179           char __tmp = *__a;                    \
180           *__a++ = *__b;                        \
181           *__b++ = __tmp;                       \
182         } while (--__size > 0);                 \
183     } while (0)
184
185 /* Makes the elements in ARRAY unique, by moving up duplicates,
186    and returns the new number of elements in the array.  Sorted
187    arrays only.  Arguments same as for sort() above. */
188 size_t
189 unique (void *array, size_t count, size_t size,
190         algo_compare_func *compare, const void *aux)
191 {
192   char *first = array;
193   char *last = first + size * count;
194   char *result = array;
195
196   for (;;)
197     {
198       first += size;
199       if (first >= last)
200         {
201           assert (adjacent_find_equal (array, count,
202                                        size, compare, aux) == NULL);
203           return count;
204         }
205
206       if (compare (result, first, aux))
207         {
208           result += size;
209           if (result != first)
210             memcpy (result, first, size);
211         }
212       else
213         count--;
214     }
215 }
216
217 /* Helper function that calls sort(), then unique(). */
218 size_t
219 sort_unique (void *array, size_t count, size_t size,
220              algo_compare_func *compare, const void *aux)
221 {
222   sort (array, count, size, compare, aux);
223   return unique (array, count, size, compare, aux);
224 }
225 \f
226 /* Reorders ARRAY, which contains COUNT elements of SIZE bytes
227    each, so that the elements for which PREDICATE returns true
228    precede those for which PREDICATE returns zero.  AUX is
229    passed to each predicate as auxiliary data.  Returns the
230    number of elements for which PREDICATE returns true.  Not
231    stable. */
232 size_t
233 partition (void *array, size_t count, size_t size,
234            algo_predicate_func *predicate, const void *aux)
235 {
236   size_t true_cnt = count;
237   char *first = array;
238   char *last = first + true_cnt * size;
239
240   for (;;)
241     {
242       /* Move FIRST forward to point to first element that fails
243          PREDICATE. */
244       for (;;)
245         {
246           if (first == last)
247             goto done;
248           else if (!predicate (first, aux))
249             break;
250
251           first += size;
252         }
253       true_cnt--;
254
255       /* Move LAST backward to point to last element that passes
256          PREDICATE. */
257       for (;;)
258         {
259           last -= size;
260
261           if (first == last)
262             goto done;
263           else if (predicate (last, aux))
264             break;
265           else
266             true_cnt--;
267         }
268
269       /* By swapping FIRST and LAST we extend the starting and
270          ending sequences that pass and fail, respectively,
271          PREDICATE. */
272       SWAP (first, last, size);
273       first += size;
274     }
275
276  done:
277   assert (is_partitioned (array, count, size, true_cnt, predicate, aux));
278   return true_cnt;
279 }
280
281 /* Checks whether ARRAY, which contains COUNT elements of SIZE
282    bytes each, is partitioned such that PREDICATE returns true
283    for the first TRUE_CNT elements and zero for the remaining
284    elements.  AUX is passed as auxiliary data to PREDICATE. */
285 bool
286 is_partitioned (const void *array, size_t count, size_t size,
287                 size_t true_cnt,
288                 algo_predicate_func *predicate, const void *aux)
289 {
290   const char *first = array;
291   size_t idx;
292
293   assert (true_cnt <= count);
294   for (idx = 0; idx < true_cnt; idx++)
295     if (predicate (first + idx * size, aux) == 0)
296       return false;
297   for (idx = true_cnt; idx < count; idx++)
298     if (predicate (first + idx * size, aux) != 0)
299       return false;
300   return true;
301 }
302 \f
303 /* Copies the COUNT elements of SIZE bytes each from ARRAY to
304    RESULT, except that elements for which PREDICATE is false are
305    not copied.  Returns the number of elements copied.  AUX is
306    passed to PREDICATE as auxiliary data.  */
307 size_t
308 copy_if (const void *array, size_t count, size_t size,
309          void *result,
310          algo_predicate_func *predicate, const void *aux)
311 {
312   const char *input = array;
313   const char *last = input + size * count;
314   char *output = result;
315   size_t nonzero_cnt = 0;
316
317   while (input < last)
318     {
319       if (predicate (input, aux))
320         {
321           memcpy (output, input, size);
322           output += size;
323           nonzero_cnt++;
324         }
325
326       input += size;
327     }
328
329   assert (nonzero_cnt == count_if (array, count, size, predicate, aux));
330   assert (nonzero_cnt == count_if (result, nonzero_cnt, size, predicate, aux));
331
332   return nonzero_cnt;
333 }
334
335 /* Removes N elements starting at IDX from ARRAY, which consists
336    of COUNT elements of SIZE bytes each, by shifting the elements
337    following them, if any, into its position. */
338 void
339 remove_range (void *array_, size_t count, size_t size,
340               size_t idx, size_t n)
341 {
342   char *array = array_;
343
344   assert (array != NULL);
345   assert (idx <= count);
346   assert (idx + n <= count);
347
348   if (idx + n < count)
349     memmove (array + idx * size, array + (idx + n) * size,
350              size * (count - idx - n));
351 }
352
353 /* Removes element IDX from ARRAY, which consists of COUNT
354    elements of SIZE bytes each, by shifting the elements
355    following it, if any, into its position. */
356 void
357 remove_element (void *array, size_t count, size_t size,
358                 size_t idx)
359 {
360   remove_range (array, count, size, idx, 1);
361 }
362
363 /* Makes room for N elements starting at IDX in ARRAY, which
364    initially consists of COUNT elements of SIZE bytes each, by
365    shifting elements IDX...COUNT (exclusive) to the right by N
366    positions. */
367 void
368 insert_range (void *array_, size_t count, size_t size,
369               size_t idx, size_t n)
370 {
371   char *array = array_;
372
373   assert (idx <= count);
374   memmove (array + (idx + n) * size, array + idx * size, (count - idx) * size);
375 }
376
377 /* Makes room for a new element at IDX in ARRAY, which initially
378    consists of COUNT elements of SIZE bytes each, by shifting
379    elements IDX...COUNT (exclusive) to the right by one
380    positions. */
381 void
382 insert_element (void *array, size_t count, size_t size,
383                 size_t idx)
384 {
385   insert_range (array, count, size, idx, 1);
386 }
387
388 /* Moves an element in ARRAY, which consists of COUNT elements of
389    SIZE bytes each, from OLD_IDX to NEW_IDX, shifting around
390    other elements as needed.  Runs in O(abs(OLD_IDX - NEW_IDX))
391    time. */
392 void
393 move_element (void *array_, size_t count, size_t size,
394               size_t old_idx, size_t new_idx)
395 {
396   assert (array_ != NULL || count == 0);
397   assert (old_idx < count);
398   assert (new_idx < count);
399
400   if (old_idx != new_idx)
401     {
402       char *array = array_;
403       char *element = xmalloc (size);
404       char *new = array + new_idx * size;
405       char *old = array + old_idx * size;
406
407       memcpy (element, old, size);
408       if (new < old)
409         memmove (new + size, new, (old_idx - new_idx) * size);
410       else
411         memmove (old, old + size, (new_idx - old_idx) * size);
412       memcpy (new, element, size);
413
414       free (element);
415     }
416 }
417
418 /* Moves N elements in ARRAY starting at OLD_IDX, which consists
419    of COUNT elements of SIZE bytes each, so that they now start
420    at NEW_IDX, shifting around other elements as needed. */
421 void
422 move_range (void *array_, size_t count, size_t size,
423             size_t old_idx, size_t new_idx, size_t n)
424 {
425   assert (array_ != NULL || count == 0);
426   assert (n <= count);
427   assert (old_idx + n <= count);
428   assert (new_idx + n <= count);
429
430   if (old_idx != new_idx && n > 0)
431     {
432       char *array = array_;
433       char *range = xmalloc (size * n);
434       char *new = array + new_idx * size;
435       char *old = array + old_idx * size;
436
437       memcpy (range, old, size * n);
438       if (new < old)
439         memmove (new + size * n, new, (old_idx - new_idx) * size);
440       else
441         memmove (old, old + size * n, (new_idx - old_idx) * size);
442       memcpy (new, range, size * n);
443
444       free (range);
445     }
446 }
447
448 /* A predicate and its auxiliary data. */
449 struct pred_aux
450   {
451     algo_predicate_func *predicate;
452     const void *aux;
453   };
454
455 static bool
456 not (const void *data, const void *pred_aux_)
457 {
458   const struct pred_aux *pred_aux = pred_aux_;
459
460   return !pred_aux->predicate (data, pred_aux->aux);
461 }
462
463 /* Removes elements equal to ELEMENT from ARRAY, which consists
464    of COUNT elements of SIZE bytes each.  Returns the number of
465    remaining elements.  AUX is passed to COMPARE as auxiliary
466    data. */
467 size_t
468 remove_equal (void *array, size_t count, size_t size,
469               void *element,
470               algo_compare_func *compare, const void *aux)
471 {
472   char *first = array;
473   char *last = first + count * size;
474   char *result;
475
476   for (;;)
477     {
478       if (first >= last)
479         goto done;
480       if (compare (first, element, aux) == 0)
481         break;
482
483       first += size;
484     }
485
486   result = first;
487   count--;
488   for (;;)
489     {
490       first += size;
491       if (first >= last)
492         goto done;
493
494       if (compare (first, element, aux) == 0)
495         {
496           count--;
497           continue;
498         }
499
500       memcpy (result, first, size);
501       result += size;
502     }
503
504  done:
505   assert (count_equal (array, count, size, element, compare, aux) == 0);
506   return count;
507 }
508
509 /* Copies the COUNT elements of SIZE bytes each from ARRAY to
510    RESULT, except that elements for which PREDICATE is true are
511    not copied.  Returns the number of elements copied.  AUX is
512    passed to PREDICATE as auxiliary data.  */
513 size_t
514 remove_copy_if (const void *array, size_t count, size_t size,
515                 void *result,
516                 algo_predicate_func *predicate, const void *aux)
517 {
518   struct pred_aux pred_aux;
519   pred_aux.predicate = predicate;
520   pred_aux.aux = aux;
521   return copy_if (array, count, size, result, not, &pred_aux);
522 }
523 \f
524 /* Searches ARRAY, which contains COUNT of SIZE bytes each, using
525    a binary search.  Returns any element that equals VALUE, if
526    one exists, or a null pointer otherwise.  ARRAY must ordered
527    according to COMPARE.  AUX is passed to COMPARE as auxiliary
528    data. */
529 void *
530 binary_search (const void *array, size_t count, size_t size,
531                void *value,
532                algo_compare_func *compare, const void *aux)
533 {
534   assert (array != NULL || count == 0);
535   assert (count <= INT_MAX);
536   assert (compare != NULL);
537
538   if (count != 0)
539     {
540       const char *first = array;
541       int low = 0;
542       int high = count - 1;
543
544       while (low <= high)
545         {
546           int middle = (low + high) / 2;
547           const char *element = first + middle * size;
548           int cmp = compare (value, element, aux);
549
550           if (cmp > 0)
551             low = middle + 1;
552           else if (cmp < 0)
553             high = middle - 1;
554           else
555             return (void *) element;
556         }
557     }
558
559   expensive_assert (find (array, count, size, value, compare, aux) == NULL);
560   return NULL;
561 }
562 \f
563 /* Lexicographically compares ARRAY1, which contains COUNT1
564    elements of SIZE bytes each, to ARRAY2, which contains COUNT2
565    elements of SIZE bytes, according to COMPARE.  Returns a
566    strcmp()-type result.  AUX is passed to COMPARE as auxiliary
567    data. */
568 int
569 lexicographical_compare_3way (const void *array1, size_t count1,
570                               const void *array2, size_t count2,
571                               size_t size,
572                               algo_compare_func *compare, const void *aux)
573 {
574   const char *first1 = array1;
575   const char *first2 = array2;
576   size_t min_count = count1 < count2 ? count1 : count2;
577
578   while (min_count > 0)
579     {
580       int cmp = compare (first1, first2, aux);
581       if (cmp != 0)
582         return cmp;
583
584       first1 += size;
585       first2 += size;
586       min_count--;
587     }
588
589   return count1 < count2 ? -1 : count1 > count2;
590 }
591 \f
592 /* If you consider tuning this algorithm, you should consult first:
593    Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
594    Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */
595
596 #include <limits.h>
597 #include <stdlib.h>
598 #include <string.h>
599
600 /* Discontinue quicksort algorithm when partition gets below this size.
601    This particular magic number was chosen to work best on a Sun 4/260. */
602 #define MAX_THRESH 4
603
604 /* Stack node declarations used to store unfulfilled partition obligations. */
605 typedef struct
606   {
607     char *lo;
608     char *hi;
609   } stack_node;
610
611 /* The next 4 #defines implement a very fast in-line stack abstraction. */
612 /* The stack needs log (total_elements) entries (we could even subtract
613    log(MAX_THRESH)).  Since total_elements has type size_t, we get as
614    upper bound for log (total_elements):
615    bits per byte (CHAR_BIT) * sizeof(size_t).  */
616 #define STACK_SIZE      (CHAR_BIT * sizeof(size_t))
617 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
618 #define POP(low, high)  ((void) (--top, (low = top->lo), (high = top->hi)))
619 #define STACK_NOT_EMPTY (stack < top)
620
621
622 /* Order size using quicksort.  This implementation incorporates
623    four optimizations discussed in Sedgewick:
624
625    1. Non-recursive, using an explicit stack of pointer that store the
626       next array partition to sort.  To save time, this maximum amount
627       of space required to store an array of SIZE_MAX is allocated on the
628       stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
629       only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
630       Pretty cheap, actually.
631
632    2. Chose the pivot element using a median-of-three decision tree.
633       This reduces the probability of selecting a bad pivot value and
634       eliminates certain extraneous comparisons.
635
636    3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
637       insertion sort to order the MAX_THRESH items within each partition.
638       This is a big win, since insertion sort is faster for small, mostly
639       sorted array segments.
640
641    4. The larger of the two sub-partitions is always pushed onto the
642       stack first, with the algorithm then concentrating on the
643       smaller partition.  This *guarantees* no more than log (total_elems)
644       stack size is needed (actually O(1) in this case)!  */
645
646 void
647 sort (void *array, size_t count, size_t size,
648             algo_compare_func *compare, const void *aux)
649 {
650   char *const first = array;
651   const size_t max_thresh = MAX_THRESH * size;
652
653   if (count == 0)
654     /* Avoid lossage with unsigned arithmetic below.  */
655     return;
656
657   if (count > MAX_THRESH)
658     {
659       char *lo = first;
660       char *hi = &lo[size * (count - 1)];
661       stack_node stack[STACK_SIZE];
662       stack_node *top = stack + 1;
663
664       while (STACK_NOT_EMPTY)
665         {
666           char *left_ptr;
667           char *right_ptr;
668
669           /* Select median value from among LO, MID, and HI. Rearrange
670              LO and HI so the three values are sorted. This lowers the
671              probability of picking a pathological pivot value and
672              skips a comparison for both the LEFT_PTR and RIGHT_PTR in
673              the while loops. */
674
675           char *mid = lo + size * ((hi - lo) / size >> 1);
676
677           if (compare (mid, lo, aux) < 0)
678             SWAP (mid, lo, size);
679           if (compare (hi, mid, aux) < 0)
680             SWAP (mid, hi, size);
681           else
682             goto jump_over;
683           if (compare (mid, lo, aux) < 0)
684             SWAP (mid, lo, size);
685         jump_over:;
686
687           left_ptr  = lo + size;
688           right_ptr = hi - size;
689
690           /* Here's the famous ``collapse the walls'' section of quicksort.
691              Gotta like those tight inner loops!  They are the main reason
692              that this algorithm runs much faster than others. */
693           do
694             {
695               while (compare (left_ptr, mid, aux) < 0)
696                 left_ptr += size;
697
698               while (compare (mid, right_ptr, aux) < 0)
699                 right_ptr -= size;
700
701               if (left_ptr < right_ptr)
702                 {
703                   SWAP (left_ptr, right_ptr, size);
704                   if (mid == left_ptr)
705                     mid = right_ptr;
706                   else if (mid == right_ptr)
707                     mid = left_ptr;
708                   left_ptr += size;
709                   right_ptr -= size;
710                 }
711               else if (left_ptr == right_ptr)
712                 {
713                   left_ptr += size;
714                   right_ptr -= size;
715                   break;
716                 }
717             }
718           while (left_ptr <= right_ptr);
719
720           /* Set up pointers for next iteration.  First determine whether
721              left and right partitions are below the threshold size.  If so,
722              ignore one or both.  Otherwise, push the larger partition's
723              bounds on the stack and continue sorting the smaller one. */
724
725           if ((size_t) (right_ptr - lo) <= max_thresh)
726             {
727               if ((size_t) (hi - left_ptr) <= max_thresh)
728                 /* Ignore both small partitions. */
729                 POP (lo, hi);
730               else
731                 /* Ignore small left partition. */
732                 lo = left_ptr;
733             }
734           else if ((size_t) (hi - left_ptr) <= max_thresh)
735             /* Ignore small right partition. */
736             hi = right_ptr;
737           else if ((right_ptr - lo) > (hi - left_ptr))
738             {
739               /* Push larger left partition indices. */
740               PUSH (lo, right_ptr);
741               lo = left_ptr;
742             }
743           else
744             {
745               /* Push larger right partition indices. */
746               PUSH (left_ptr, hi);
747               hi = right_ptr;
748             }
749         }
750     }
751
752   /* Once the FIRST array is partially sorted by quicksort the rest
753      is completely sorted using insertion sort, since this is efficient
754      for partitions below MAX_THRESH size. FIRST points to the beginning
755      of the array to sort, and END_PTR points at the very last element in
756      the array (*not* one beyond it!). */
757
758   {
759     char *const end_ptr = &first[size * (count - 1)];
760     char *tmp_ptr = first;
761     char *thresh = MIN (end_ptr, first + max_thresh);
762     register char *run_ptr;
763
764     /* Find smallest element in first threshold and place it at the
765        array's beginning.  This is the smallest array element,
766        and the operation speeds up insertion sort's inner loop. */
767
768     for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
769       if (compare (run_ptr, tmp_ptr, aux) < 0)
770         tmp_ptr = run_ptr;
771
772     if (tmp_ptr != first)
773       SWAP (tmp_ptr, first, size);
774
775     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
776
777     run_ptr = first + size;
778     while ((run_ptr += size) <= end_ptr)
779       {
780         tmp_ptr = run_ptr - size;
781         while (compare (run_ptr, tmp_ptr, aux) < 0)
782           tmp_ptr -= size;
783
784         tmp_ptr += size;
785         if (tmp_ptr != run_ptr)
786           {
787             char *trav;
788
789             trav = run_ptr + size;
790             while (--trav >= run_ptr)
791               {
792                 char c = *trav;
793                 char *hi, *lo;
794
795                 for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
796                   *hi = *lo;
797                 *hi = c;
798               }
799           }
800       }
801   }
802
803   assert (is_sorted (array, count, size, compare, aux));
804 }
805
806 /* Tests whether ARRAY, which contains COUNT elements of SIZE
807    bytes each, is sorted in order according to COMPARE.  AUX is
808    passed to COMPARE as auxiliary data. */
809 bool
810 is_sorted (const void *array, size_t count, size_t size,
811            algo_compare_func *compare, const void *aux)
812 {
813   const char *first = array;
814   size_t idx;
815
816   for (idx = 0; idx + 1 < count; idx++)
817     if (compare (first + idx * size, first + (idx + 1) * size, aux) > 0)
818       return false;
819
820   return true;
821 }
822 \f
823 /* Computes the generalized set difference, ARRAY1 minus ARRAY2,
824    into RESULT, and returns the number of elements written to
825    RESULT.  If a value appears M times in ARRAY1 and N times in
826    ARRAY2, then it will appear max(M - N, 0) in RESULT.  ARRAY1
827    and ARRAY2 must be sorted, and RESULT is sorted and stable.
828    ARRAY1 consists of COUNT1 elements, ARRAY2 of COUNT2 elements,
829    each SIZE bytes.  AUX is passed to COMPARE as auxiliary
830    data. */
831 size_t set_difference (const void *array1, size_t count1,
832                        const void *array2, size_t count2,
833                        size_t size,
834                        void *result_,
835                        algo_compare_func *compare, const void *aux)
836 {
837   const char *first1 = array1;
838   const char *last1 = first1 + count1 * size;
839   const char *first2 = array2;
840   const char *last2 = first2 + count2 * size;
841   char *result = result_;
842   size_t result_count = 0;
843
844   while (first1 != last1 && first2 != last2)
845     {
846       int cmp = compare (first1, first2, aux);
847       if (cmp < 0)
848         {
849           memcpy (result, first1, size);
850           first1 += size;
851           result += size;
852           result_count++;
853         }
854       else if (cmp > 0)
855         first2 += size;
856       else
857         {
858           first1 += size;
859           first2 += size;
860         }
861     }
862
863   while (first1 != last1)
864     {
865       memcpy (result, first1, size);
866       first1 += size;
867       result += size;
868       result_count++;
869     }
870
871   return result_count;
872 }
873 \f
874 /* Finds the first pair of adjacent equal elements in ARRAY,
875    which has COUNT elements of SIZE bytes.  Returns the first
876    element in ARRAY such that COMPARE returns zero when it and
877    its successor element are compared, or a null pointer if no
878    such element exists.  AUX is passed to COMPARE as auxiliary
879    data. */
880 void *
881 adjacent_find_equal (const void *array, size_t count, size_t size,
882                      algo_compare_func *compare, const void *aux)
883 {
884   const char *first = array;
885   const char *last = first + count * size;
886
887   while (first < last && first + size < last)
888     {
889       if (compare (first, first + size, aux) == 0)
890         return (void *) first;
891       first += size;
892     }
893
894   return NULL;
895 }
896 \f
897 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
898    the first COUNT - 1 elements of these form a heap, followed by
899    a single element not part of the heap.  This function adds the
900    final element, forming a heap of COUNT elements in ARRAY.
901    Uses COMPARE to compare elements, passing AUX as auxiliary
902    data. */
903 void
904 push_heap (void *array, size_t count, size_t size,
905            algo_compare_func *compare, const void *aux)
906 {
907   char *first = array;
908   size_t i;
909
910   expensive_assert (count < 1 || is_heap (array, count - 1,
911                                           size, compare, aux));
912   for (i = count; i > 1; i /= 2)
913     {
914       char *parent = first + (i / 2 - 1) * size;
915       char *element = first + (i - 1) * size;
916       if (compare (parent, element, aux) < 0)
917         SWAP (parent, element, size);
918       else
919         break;
920     }
921   expensive_assert (is_heap (array, count, size, compare, aux));
922 }
923
924 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
925    the children of ARRAY[idx - 1] are heaps, but ARRAY[idx - 1]
926    may be smaller than its children.  This function fixes that,
927    so that ARRAY[idx - 1] itself is a heap.  Uses COMPARE to
928    compare elements, passing AUX as auxiliary data. */
929 static void
930 heapify (void *array, size_t count, size_t size,
931          size_t idx,
932          algo_compare_func *compare, const void *aux)
933 {
934   char *first = array;
935
936   for (;;)
937     {
938       size_t left = 2 * idx;
939       size_t right = left + 1;
940       size_t largest = idx;
941
942       if (left <= count
943           && compare (first + size * (left - 1),
944                       first + size * (idx - 1), aux) > 0)
945         largest = left;
946
947       if (right <= count
948           && compare (first + size * (right - 1),
949                       first + size * (largest - 1), aux) > 0)
950         largest = right;
951
952       if (largest == idx)
953         break;
954
955       SWAP (first + size * (idx - 1), first + size * (largest - 1), size);
956       idx = largest;
957     }
958 }
959
960 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
961    all COUNT elements form a heap.  This function moves the
962    largest element in the heap to the final position in ARRAY and
963    reforms a heap of the remaining COUNT - 1 elements at the
964    beginning of ARRAY.  Uses COMPARE to compare elements, passing
965    AUX as auxiliary data. */
966 void
967 pop_heap (void *array, size_t count, size_t size,
968           algo_compare_func *compare, const void *aux)
969 {
970   char *first = array;
971
972   expensive_assert (is_heap (array, count, size, compare, aux));
973   SWAP (first, first + (count - 1) * size, size);
974   heapify (first, count - 1, size, 1, compare, aux);
975   expensive_assert (count < 1 || is_heap (array, count - 1,
976                                           size, compare, aux));
977 }
978
979 /* Turns ARRAY, which contains COUNT elements of SIZE bytes, into
980    a heap.  Uses COMPARE to compare elements, passing AUX as
981    auxiliary data. */
982 void
983 make_heap (void *array, size_t count, size_t size,
984            algo_compare_func *compare, const void *aux)
985 {
986   size_t idx;
987
988   for (idx = count / 2; idx >= 1; idx--)
989     heapify (array, count, size, idx, compare, aux);
990   expensive_assert (count < 1 || is_heap (array, count, size, compare, aux));
991 }
992
993 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
994    all COUNT elements form a heap.  This function turns the heap
995    into a fully sorted array.  Uses COMPARE to compare elements,
996    passing AUX as auxiliary data. */
997 void
998 sort_heap (void *array, size_t count, size_t size,
999            algo_compare_func *compare, const void *aux)
1000 {
1001   char *first = array;
1002   size_t idx;
1003
1004   expensive_assert (is_heap (array, count, size, compare, aux));
1005   for (idx = count; idx >= 2; idx--)
1006     {
1007       SWAP (first, first + (idx - 1) * size, size);
1008       heapify (array, idx - 1, size, 1, compare, aux);
1009     }
1010   expensive_assert (is_sorted (array, count, size, compare, aux));
1011 }
1012
1013 /* ARRAY contains COUNT elements of SIZE bytes each.  This
1014    function tests whether ARRAY is a heap and returns true if so,
1015    false otherwise.  Uses COMPARE to compare elements, passing
1016    AUX as auxiliary data. */
1017 bool
1018 is_heap (const void *array, size_t count, size_t size,
1019          algo_compare_func *compare, const void *aux)
1020 {
1021   const char *first = array;
1022   size_t child;
1023
1024   for (child = 2; child <= count; child++)
1025     {
1026       size_t parent = child / 2;
1027       if (compare (first + (parent - 1) * size,
1028                    first + (child - 1) * size, aux) < 0)
1029         return false;
1030     }
1031
1032   return true;
1033 }
1034