209b0b0f63cce3e3c256adf8f4ef9293ac0fdd1a
[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 /* A predicate and its auxiliary data. */
397 struct pred_aux 
398   {
399     algo_predicate_func *predicate;
400     void *aux;
401   };
402
403 static int
404 not (const void *data, void *pred_aux_) 
405 {
406   const struct pred_aux *pred_aux = pred_aux_;
407
408   return !pred_aux->predicate (data, pred_aux->aux);
409 }
410
411 /* Removes elements equal to ELEMENT from ARRAY, which consists
412    of COUNT elements of SIZE bytes each.  Returns the number of
413    remaining elements.  AUX is passed to COMPARE as auxiliary
414    data. */
415 size_t
416 remove_equal (void *array, size_t count, size_t size,
417               void *element,
418               algo_compare_func *compare, void *aux) 
419 {
420   unsigned char *first = array;
421   unsigned char *last = first + count * size;
422   unsigned char *result;
423
424   for (;;)
425     {
426       if (first >= last)
427         goto done;
428       if (compare (first, element, aux) == 0)
429         break;
430
431       first += size;
432     }
433
434   result = first;
435   count--;
436   for (;;) 
437     {
438       first += size;
439       if (first >= last)
440         goto done;
441
442       if (compare (first, element, aux) == 0) 
443         {
444           count--; 
445           continue;
446         }
447       
448       memcpy (result, first, size);
449       result += size;
450     }
451
452  done:
453   assert (count_equal (array, count, size, element, compare, aux) == 0);
454   return count;
455 }
456
457 /* Copies the COUNT elements of SIZE bytes each from ARRAY to
458    RESULT, except that elements for which PREDICATE is true are
459    not copied.  Returns the number of elements copied.  AUX is
460    passed to PREDICATE as auxiliary data.  */
461 size_t 
462 remove_copy_if (const void *array, size_t count, size_t size,
463                 void *result,
464                 algo_predicate_func *predicate, void *aux) 
465 {
466   struct pred_aux pred_aux;
467   pred_aux.predicate = predicate;
468   pred_aux.aux = aux;
469   return copy_if (array, count, size, result, not, &pred_aux);
470 }
471 \f
472 /* Searches ARRAY, which contains COUNT of SIZE bytes each, using
473    a binary search.  Returns any element that equals VALUE, if
474    one exists, or a null pointer otherwise.  ARRAY must ordered
475    according to COMPARE.  AUX is passed to COMPARE as auxiliary
476    data. */
477 void *
478 binary_search (const void *array, size_t count, size_t size,
479                void *value,
480                algo_compare_func *compare, void *aux) 
481 {
482   assert (array != NULL);
483   assert (count <= INT_MAX);
484   assert (compare != NULL);
485
486   if (count != 0) 
487     {
488       const unsigned char *first = array;
489       int low = 0;
490       int high = count - 1;
491
492       while (low <= high) 
493         {
494           int middle = (low + high) / 2;
495           const unsigned char *element = first + middle * size;
496           int cmp = compare (value, element, aux);
497
498           if (cmp > 0) 
499             low = middle + 1;
500           else if (cmp < 0)
501             high = middle - 1;
502           else
503             return (void *) element;
504         }
505     }
506
507   expensive_assert (find (array, count, size, value, compare, aux) == NULL);
508   return NULL;
509 }
510 \f
511 /* Lexicographically compares ARRAY1, which contains COUNT1
512    elements of SIZE bytes each, to ARRAY2, which contains COUNT2
513    elements of SIZE bytes, according to COMPARE.  Returns a
514    strcmp()-type result.  AUX is passed to COMPARE as auxiliary
515    data. */
516 int
517 lexicographical_compare_3way (const void *array1, size_t count1,
518                               const void *array2, size_t count2,
519                               size_t size,
520                               algo_compare_func *compare, void *aux) 
521 {
522   const unsigned char *first1 = array1;
523   const unsigned char *first2 = array2;
524   size_t min_count = count1 < count2 ? count1 : count2;
525
526   while (min_count > 0)
527     {
528       int cmp = compare (first1, first2, aux);
529       if (cmp != 0)
530         return cmp;
531
532       first1 += size;
533       first2 += size;
534       min_count--;
535     }
536
537   return count1 < count2 ? -1 : count1 > count2;
538 }
539 \f
540 /* If you consider tuning this algorithm, you should consult first:
541    Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
542    Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */
543
544 #ifdef HAVE_ALLOCA_H
545 #include <alloca.h>
546 #endif
547
548 #include <limits.h>
549 #include <stdlib.h>
550 #include <string.h>
551
552 /* Discontinue quicksort algorithm when partition gets below this size.
553    This particular magic number was chosen to work best on a Sun 4/260. */
554 #define MAX_THRESH 4
555
556 /* Stack node declarations used to store unfulfilled partition obligations. */
557 typedef struct
558   {
559     char *lo;
560     char *hi;
561   } stack_node;
562
563 /* The next 4 #defines implement a very fast in-line stack abstraction. */
564 /* The stack needs log (total_elements) entries (we could even subtract
565    log(MAX_THRESH)).  Since total_elements has type size_t, we get as
566    upper bound for log (total_elements):
567    bits per byte (CHAR_BIT) * sizeof(size_t).  */
568 #define STACK_SIZE      (CHAR_BIT * sizeof(size_t))
569 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
570 #define POP(low, high)  ((void) (--top, (low = top->lo), (high = top->hi)))
571 #define STACK_NOT_EMPTY (stack < top)
572
573
574 /* Order size using quicksort.  This implementation incorporates
575    four optimizations discussed in Sedgewick:
576
577    1. Non-recursive, using an explicit stack of pointer that store the
578       next array partition to sort.  To save time, this maximum amount
579       of space required to store an array of SIZE_MAX is allocated on the
580       stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
581       only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
582       Pretty cheap, actually.
583
584    2. Chose the pivot element using a median-of-three decision tree.
585       This reduces the probability of selecting a bad pivot value and
586       eliminates certain extraneous comparisons.
587
588    3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
589       insertion sort to order the MAX_THRESH items within each partition.
590       This is a big win, since insertion sort is faster for small, mostly
591       sorted array segments.
592
593    4. The larger of the two sub-partitions is always pushed onto the
594       stack first, with the algorithm then concentrating on the
595       smaller partition.  This *guarantees* no more than log (total_elems)
596       stack size is needed (actually O(1) in this case)!  */
597
598 void
599 sort (void *array, size_t count, size_t size,
600       algo_compare_func *compare, void *aux)
601 {
602   char *const first = array;
603   const size_t max_thresh = MAX_THRESH * size;
604
605   if (count == 0)
606     /* Avoid lossage with unsigned arithmetic below.  */
607     return;
608
609   if (count > MAX_THRESH)
610     {
611       char *lo = first;
612       char *hi = &lo[size * (count - 1)];
613       stack_node stack[STACK_SIZE];
614       stack_node *top = stack + 1;
615
616       while (STACK_NOT_EMPTY)
617         {
618           char *left_ptr;
619           char *right_ptr;
620
621           /* Select median value from among LO, MID, and HI. Rearrange
622              LO and HI so the three values are sorted. This lowers the
623              probability of picking a pathological pivot value and
624              skips a comparison for both the LEFT_PTR and RIGHT_PTR in
625              the while loops. */
626
627           char *mid = lo + size * ((hi - lo) / size >> 1);
628
629           if (compare (mid, lo, aux) < 0)
630             SWAP (mid, lo, size);
631           if (compare (hi, mid, aux) < 0)
632             SWAP (mid, hi, size);
633           else
634             goto jump_over;
635           if (compare (mid, lo, aux) < 0)
636             SWAP (mid, lo, size);
637         jump_over:;
638
639           left_ptr  = lo + size;
640           right_ptr = hi - size;
641
642           /* Here's the famous ``collapse the walls'' section of quicksort.
643              Gotta like those tight inner loops!  They are the main reason
644              that this algorithm runs much faster than others. */
645           do
646             {
647               while (compare (left_ptr, mid, aux) < 0)
648                 left_ptr += size;
649
650               while (compare (mid, right_ptr, aux) < 0)
651                 right_ptr -= size;
652
653               if (left_ptr < right_ptr)
654                 {
655                   SWAP (left_ptr, right_ptr, size);
656                   if (mid == left_ptr)
657                     mid = right_ptr;
658                   else if (mid == right_ptr)
659                     mid = left_ptr;
660                   left_ptr += size;
661                   right_ptr -= size;
662                 }
663               else if (left_ptr == right_ptr)
664                 {
665                   left_ptr += size;
666                   right_ptr -= size;
667                   break;
668                 }
669             }
670           while (left_ptr <= right_ptr);
671
672           /* Set up pointers for next iteration.  First determine whether
673              left and right partitions are below the threshold size.  If so,
674              ignore one or both.  Otherwise, push the larger partition's
675              bounds on the stack and continue sorting the smaller one. */
676
677           if ((size_t) (right_ptr - lo) <= max_thresh)
678             {
679               if ((size_t) (hi - left_ptr) <= max_thresh)
680                 /* Ignore both small partitions. */
681                 POP (lo, hi);
682               else
683                 /* Ignore small left partition. */
684                 lo = left_ptr;
685             }
686           else if ((size_t) (hi - left_ptr) <= max_thresh)
687             /* Ignore small right partition. */
688             hi = right_ptr;
689           else if ((right_ptr - lo) > (hi - left_ptr))
690             {
691               /* Push larger left partition indices. */
692               PUSH (lo, right_ptr);
693               lo = left_ptr;
694             }
695           else
696             {
697               /* Push larger right partition indices. */
698               PUSH (left_ptr, hi);
699               hi = right_ptr;
700             }
701         }
702     }
703
704   /* Once the FIRST array is partially sorted by quicksort the rest
705      is completely sorted using insertion sort, since this is efficient
706      for partitions below MAX_THRESH size. FIRST points to the beginning
707      of the array to sort, and END_PTR points at the very last element in
708      the array (*not* one beyond it!). */
709
710 #define min(x, y) ((x) < (y) ? (x) : (y))
711
712   {
713     char *const end_ptr = &first[size * (count - 1)];
714     char *tmp_ptr = first;
715     char *thresh = min(end_ptr, first + max_thresh);
716     register char *run_ptr;
717
718     /* Find smallest element in first threshold and place it at the
719        array's beginning.  This is the smallest array element,
720        and the operation speeds up insertion sort's inner loop. */
721
722     for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
723       if (compare (run_ptr, tmp_ptr, aux) < 0)
724         tmp_ptr = run_ptr;
725
726     if (tmp_ptr != first)
727       SWAP (tmp_ptr, first, size);
728
729     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
730
731     run_ptr = first + size;
732     while ((run_ptr += size) <= end_ptr)
733       {
734         tmp_ptr = run_ptr - size;
735         while (compare (run_ptr, tmp_ptr, aux) < 0)
736           tmp_ptr -= size;
737
738         tmp_ptr += size;
739         if (tmp_ptr != run_ptr)
740           {
741             char *trav;
742
743             trav = run_ptr + size;
744             while (--trav >= run_ptr)
745               {
746                 char c = *trav;
747                 char *hi, *lo;
748
749                 for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
750                   *hi = *lo;
751                 *hi = c;
752               }
753           }
754       }
755   }
756
757   assert (is_sorted (array, count, size, compare, aux));
758 }
759
760 /* Tests whether ARRAY, which contains COUNT elements of SIZE
761    bytes each, is sorted in order according to COMPARE.  AUX is
762    passed to COMPARE as auxiliary data. */
763 int
764 is_sorted (const void *array, size_t count, size_t size,
765            algo_compare_func *compare, void *aux) 
766 {
767   const unsigned char *first = array;
768   size_t idx;
769       
770   for (idx = 0; idx + 1 < count; idx++)
771     if (compare (first + idx * size, first + (idx + 1) * size, aux) > 0)
772       return 0; 
773   
774   return 1;
775 }
776 \f
777 /* Computes the generalized set difference, ARRAY1 minus ARRAY2,
778    into RESULT, and returns the number of elements written to
779    RESULT.  If a value appears M times in ARRAY1 and N times in
780    ARRAY2, then it will appear max(M - N, 0) in RESULT.  ARRAY1
781    and ARRAY2 must be sorted, and RESULT is sorted and stable.
782    ARRAY1 consists of COUNT1 elements, ARRAY2 of COUNT2 elements,
783    each SIZE bytes.  AUX is passed to COMPARE as auxiliary
784    data. */
785 size_t set_difference (const void *array1, size_t count1,
786                        const void *array2, size_t count2,
787                        size_t size,
788                        void *result_,
789                        algo_compare_func *compare, void *aux) 
790 {
791   const unsigned char *first1 = array1;
792   const unsigned char *last1 = first1 + count1 * size;
793   const unsigned char *first2 = array2;
794   const unsigned char *last2 = first2 + count2 * size;
795   unsigned char *result = result_;
796   size_t result_count = 0;
797   
798   while (first1 != last1 && first2 != last2) 
799     {
800       int cmp = compare (first1, first2, aux);
801       if (cmp < 0)
802         {
803           memcpy (result, first1, size);
804           first1 += size;
805           result += size;
806           result_count++;
807         }
808       else if (cmp > 0)
809         first2 += size;
810       else
811         {
812           first1 += size;
813           first2 += size;
814         }
815     }
816
817   while (first1 != last1) 
818     {
819       memcpy (result, first1, size);
820       first1 += size;
821       result += size;
822       result_count++;
823     }
824
825   return result_count;
826 }
827 \f
828 /* Finds the first pair of adjacent equal elements in ARRAY,
829    which has COUNT elements of SIZE bytes.  Returns the first
830    element in ARRAY such that COMPARE returns zero when it and
831    its successor element are compared, or a null pointer if no
832    such element exists.  AUX is passed to COMPARE as auxiliary
833    data. */
834 void *
835 adjacent_find_equal (const void *array, size_t count, size_t size,
836                      algo_compare_func *compare, void *aux) 
837 {
838   const unsigned char *first = array;
839   const unsigned char *last = first + count * size;
840
841   while (first < last && first + size < last) 
842     {
843       if (compare (first, first + size, aux) == 0)
844         return (void *) first;
845       first += size;
846     }
847
848   return NULL;
849 }
850 \f
851 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
852    the first COUNT - 1 elements of these form a heap, followed by
853    a single element not part of the heap.  This function adds the
854    final element, forming a heap of COUNT elements in ARRAY.
855    Uses COMPARE to compare elements, passing AUX as auxiliary
856    data. */
857 void
858 push_heap (void *array, size_t count, size_t size,
859            algo_compare_func *compare, void *aux) 
860 {
861   unsigned char *first = array;
862   size_t i;
863   
864   expensive_assert (count < 1 || is_heap (array, count - 1,
865                                           size, compare, aux));
866   for (i = count; i > 1; i /= 2) 
867     {
868       unsigned char *parent = first + (i / 2 - 1) * size;
869       unsigned char *element = first + (i - 1) * size;
870       if (compare (parent, element, aux) < 0)
871         SWAP (parent, element, size);
872       else
873         break; 
874     }
875   expensive_assert (is_heap (array, count, size, compare, aux));
876 }
877
878 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
879    the children of ARRAY[idx - 1] are heaps, but ARRAY[idx - 1]
880    may be smaller than its children.  This function fixes that,
881    so that ARRAY[idx - 1] itself is a heap.  Uses COMPARE to
882    compare elements, passing AUX as auxiliary data. */
883 static void
884 heapify (void *array, size_t count, size_t size,
885          size_t idx,
886          algo_compare_func *compare, void *aux) 
887 {
888   unsigned char *first = array;
889   
890   for (;;) 
891     {
892       size_t left = 2 * idx;
893       size_t right = left + 1;
894       size_t largest = idx;
895
896       if (left <= count
897           && compare (first + size * (left - 1),
898                       first + size * (idx - 1), aux) > 0)
899         largest = left;
900
901       if (right <= count
902           && compare (first + size * (right - 1),
903                       first + size * (largest - 1), aux) > 0)
904         largest = right;
905
906       if (largest == idx)
907         break;
908
909       SWAP (first + size * (idx - 1), first + size * (largest - 1), size);
910       idx = largest;
911     }
912 }
913
914 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
915    all COUNT elements form a heap.  This function moves the
916    largest element in the heap to the final position in ARRAY and
917    reforms a heap of the remaining COUNT - 1 elements at the
918    beginning of ARRAY.  Uses COMPARE to compare elements, passing
919    AUX as auxiliary data. */
920 void
921 pop_heap (void *array, size_t count, size_t size,
922           algo_compare_func *compare, void *aux) 
923 {
924   unsigned char *first = array;
925
926   expensive_assert (is_heap (array, count, size, compare, aux));
927   SWAP (first, first + (count - 1) * size, size);
928   heapify (first, count - 1, size, 1, compare, aux);
929   expensive_assert (count < 1 || is_heap (array, count - 1,
930                                           size, compare, aux));
931 }
932
933 /* Turns ARRAY, which contains COUNT elements of SIZE bytes, into
934    a heap.  Uses COMPARE to compare elements, passing AUX as
935    auxiliary data. */
936 void
937 make_heap (void *array, size_t count, size_t size,
938            algo_compare_func *compare, void *aux) 
939 {
940   size_t idx;
941   
942   for (idx = count / 2; idx >= 1; idx--)
943     heapify (array, count, size, idx, compare, aux);
944   expensive_assert (count < 1 || is_heap (array, count, size, compare, aux));
945 }
946
947 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
948    all COUNT elements form a heap.  This function turns the heap
949    into a fully sorted array.  Uses COMPARE to compare elements,
950    passing AUX as auxiliary data. */
951 void
952 sort_heap (void *array, size_t count, size_t size,
953            algo_compare_func *compare, void *aux) 
954 {
955   unsigned char *first = array;
956   size_t idx;
957
958   expensive_assert (is_heap (array, count, size, compare, aux));
959   for (idx = count; idx >= 2; idx--)
960     {
961       SWAP (first, first + (idx - 1) * size, size);
962       heapify (array, idx - 1, size, 1, compare, aux);
963     }
964   expensive_assert (is_sorted (array, count, size, compare, aux));
965 }
966
967 /* ARRAY contains COUNT elements of SIZE bytes each.  This
968    function tests whether ARRAY is a heap and returns 1 if so, 0
969    otherwise.  Uses COMPARE to compare elements, passing AUX as
970    auxiliary data. */
971 int
972 is_heap (const void *array, size_t count, size_t size,
973          algo_compare_func *compare, void *aux) 
974 {
975   const unsigned char *first = array;
976   size_t child;
977   
978   for (child = 2; child <= count; child++)
979     {
980       size_t parent = child / 2;
981       if (compare (first + (parent - 1) * size,
982                    first + (child - 1) * size, aux) < 0)
983         return 0;
984     }
985
986   return 1;
987 }
988