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