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