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