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 #ifdef HAVE_ALLOCA_H
 575 #include <alloca.h>
 576 #endif
 577
 578 #include <limits.h>
 579 #include <stdlib.h>
 580 #include <string.h>
 581
 582 /* Discontinue quicksort algorithm when partition gets below this size.
 583    This particular magic number was chosen to work best on a Sun 4/260. */
 584 #define MAX_THRESH 4
 585
 586 /* Stack node declarations used to store unfulfilled partition obligations. */
 587 typedef struct
 588   {
 589     char *lo;
 590     char *hi;
 591   } stack_node;
 592
 593 /* The next 4 #defines implement a very fast in-line stack abstraction. */
 594 /* The stack needs log (total_elements) entries (we could even subtract
 595    log(MAX_THRESH)).  Since total_elements has type size_t, we get as
 596    upper bound for log (total_elements):
 597    bits per byte (CHAR_BIT) * sizeof(size_t).  */
 598 #define STACK_SIZE      (CHAR_BIT * sizeof(size_t))
 599 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
 600 #define POP(low, high)  ((void) (--top, (low = top->lo), (high = top->hi)))
 601 #define STACK_NOT_EMPTY (stack < top)
 602
 603
 604 /* Order size using quicksort.  This implementation incorporates
 605    four optimizations discussed in Sedgewick:
 606
 607    1. Non-recursive, using an explicit stack of pointer that store the
 608       next array partition to sort.  To save time, this maximum amount
 609       of space required to store an array of SIZE_MAX is allocated on the
 610       stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
 611       only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
 612       Pretty cheap, actually.
 613
 614    2. Chose the pivot element using a median-of-three decision tree.
 615       This reduces the probability of selecting a bad pivot value and
 616       eliminates certain extraneous comparisons.
 617
 618    3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
 619       insertion sort to order the MAX_THRESH items within each partition.
 620       This is a big win, since insertion sort is faster for small, mostly
 621       sorted array segments.
 622
 623    4. The larger of the two sub-partitions is always pushed onto the
 624       stack first, with the algorithm then concentrating on the
 625       smaller partition.  This *guarantees* no more than log (total_elems)
 626       stack size is needed (actually O(1) in this case)!  */
 627
 628 void
 629 sort (void *array, size_t count, size_t size,
 630       algo_compare_func *compare, void *aux)
 631 {
 632   char *const first = array;
 633   const size_t max_thresh = MAX_THRESH * size;
 634
 635   if (count == 0)
 636     /* Avoid lossage with unsigned arithmetic below.  */
 637     return;
 638
 639   if (count > MAX_THRESH)
 640     {
 641       char *lo = first;
 642       char *hi = &lo[size * (count - 1)];
 643       stack_node stack[STACK_SIZE];
 644       stack_node *top = stack + 1;
 645
 646       while (STACK_NOT_EMPTY)
 647         {
 648           char *left_ptr;
 649           char *right_ptr;
 650
 651           /* Select median value from among LO, MID, and HI. Rearrange
 652              LO and HI so the three values are sorted. This lowers the
 653              probability of picking a pathological pivot value and
 654              skips a comparison for both the LEFT_PTR and RIGHT_PTR in
 655              the while loops. */
 656
 657           char *mid = lo + size * ((hi - lo) / size >> 1);
 658
 659           if (compare (mid, lo, aux) < 0)
 660             SWAP (mid, lo, size);
 661           if (compare (hi, mid, aux) < 0)
 662             SWAP (mid, hi, size);
 663           else
 664             goto jump_over;
 665           if (compare (mid, lo, aux) < 0)
 666             SWAP (mid, lo, size);
 667         jump_over:;
 668
 669           left_ptr  = lo + size;
 670           right_ptr = hi - size;
 671
 672           /* Here's the famous ``collapse the walls'' section of quicksort.
 673              Gotta like those tight inner loops!  They are the main reason
 674              that this algorithm runs much faster than others. */
 675           do
 676             {
 677               while (compare (left_ptr, mid, aux) < 0)
 678                 left_ptr += size;
 679
 680               while (compare (mid, right_ptr, aux) < 0)
 681                 right_ptr -= size;
 682
 683               if (left_ptr < right_ptr)
 684                 {
 685                   SWAP (left_ptr, right_ptr, size);
 686                   if (mid == left_ptr)
 687                     mid = right_ptr;
 688                   else if (mid == right_ptr)
 689                     mid = left_ptr;
 690                   left_ptr += size;
 691                   right_ptr -= size;
 692                 }
 693               else if (left_ptr == right_ptr)
 694                 {
 695                   left_ptr += size;
 696                   right_ptr -= size;
 697                   break;
 698                 }
 699             }
 700           while (left_ptr <= right_ptr);
 701
 702           /* Set up pointers for next iteration.  First determine whether
 703              left and right partitions are below the threshold size.  If so,
 704              ignore one or both.  Otherwise, push the larger partition's
 705              bounds on the stack and continue sorting the smaller one. */
 706
 707           if ((size_t) (right_ptr - lo) <= max_thresh)
 708             {
 709               if ((size_t) (hi - left_ptr) <= max_thresh)
 710                 /* Ignore both small partitions. */
 711                 POP (lo, hi);
 712               else
 713                 /* Ignore small left partition. */
 714                 lo = left_ptr;
 715             }
 716           else if ((size_t) (hi - left_ptr) <= max_thresh)
 717             /* Ignore small right partition. */
 718             hi = right_ptr;
 719           else if ((right_ptr - lo) > (hi - left_ptr))
 720             {
 721               /* Push larger left partition indices. */
 722               PUSH (lo, right_ptr);
 723               lo = left_ptr;
 724             }
 725           else
 726             {
 727               /* Push larger right partition indices. */
 728               PUSH (left_ptr, hi);
 729               hi = right_ptr;
 730             }
 731         }
 732     }
 733
 734   /* Once the FIRST array is partially sorted by quicksort the rest
 735      is completely sorted using insertion sort, since this is efficient
 736      for partitions below MAX_THRESH size. FIRST points to the beginning
 737      of the array to sort, and END_PTR points at the very last element in
 738      the array (*not* one beyond it!). */
 739
 740 #define min(x, y) ((x) < (y) ? (x) : (y))
 741
 742   {
 743     char *const end_ptr = &first[size * (count - 1)];
 744     char *tmp_ptr = first;
 745     char *thresh = min(end_ptr, first + max_thresh);
 746     register char *run_ptr;
 747
 748     /* Find smallest element in first threshold and place it at the
 749        array's beginning.  This is the smallest array element,
 750        and the operation speeds up insertion sort's inner loop. */
 751
 752     for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
 753       if (compare (run_ptr, tmp_ptr, aux) < 0)
 754         tmp_ptr = run_ptr;
 755
 756     if (tmp_ptr != first)
 757       SWAP (tmp_ptr, first, size);
 758
 759     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
 760
 761     run_ptr = first + size;
 762     while ((run_ptr += size) <= end_ptr)
 763       {
 764         tmp_ptr = run_ptr - size;
 765         while (compare (run_ptr, tmp_ptr, aux) < 0)
 766           tmp_ptr -= size;
 767
 768         tmp_ptr += size;
 769         if (tmp_ptr != run_ptr)
 770           {
 771             char *trav;
 772
 773             trav = run_ptr + size;
 774             while (--trav >= run_ptr)
 775               {
 776                 char c = *trav;
 777                 char *hi, *lo;
 778
 779                 for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
 780                   *hi = *lo;
 781                 *hi = c;
 782               }
 783           }
 784       }
 785   }
 786
 787   assert (is_sorted (array, count, size, compare, aux));
 788 }
 789
 790 /* Tests whether ARRAY, which contains COUNT elements of SIZE
 791    bytes each, is sorted in order according to COMPARE.  AUX is
 792    passed to COMPARE as auxiliary data. */
 793 int
 794 is_sorted (const void *array, size_t count, size_t size,
 795            algo_compare_func *compare, void *aux)
 796 {
 797   const unsigned char *first = array;
 798   size_t idx;
 799
 800   for (idx = 0; idx + 1 < count; idx++)
 801     if (compare (first + idx * size, first + (idx + 1) * size, aux) > 0)
 802       return 0;
 803
 804   return 1;
 805 }
 806 \f
 807 /* Computes the generalized set difference, ARRAY1 minus ARRAY2,
 808    into RESULT, and returns the number of elements written to
 809    RESULT.  If a value appears M times in ARRAY1 and N times in
 810    ARRAY2, then it will appear max(M - N, 0) in RESULT.  ARRAY1
 811    and ARRAY2 must be sorted, and RESULT is sorted and stable.
 812    ARRAY1 consists of COUNT1 elements, ARRAY2 of COUNT2 elements,
 813    each SIZE bytes.  AUX is passed to COMPARE as auxiliary
 814    data. */
 815 size_t set_difference (const void *array1, size_t count1,
 816                        const void *array2, size_t count2,
 817                        size_t size,
 818                        void *result_,
 819                        algo_compare_func *compare, void *aux)
 820 {
 821   const unsigned char *first1 = array1;
 822   const unsigned char *last1 = first1 + count1 * size;
 823   const unsigned char *first2 = array2;
 824   const unsigned char *last2 = first2 + count2 * size;
 825   unsigned char *result = result_;
 826   size_t result_count = 0;
 827
 828   while (first1 != last1 && first2 != last2)
 829     {
 830       int cmp = compare (first1, first2, aux);
 831       if (cmp < 0)
 832         {
 833           memcpy (result, first1, size);
 834           first1 += size;
 835           result += size;
 836           result_count++;
 837         }
 838       else if (cmp > 0)
 839         first2 += size;
 840       else
 841         {
 842           first1 += size;
 843           first2 += size;
 844         }
 845     }
 846
 847   while (first1 != last1)
 848     {
 849       memcpy (result, first1, size);
 850       first1 += size;
 851       result += size;
 852       result_count++;
 853     }
 854
 855   return result_count;
 856 }
 857 \f
 858 /* Finds the first pair of adjacent equal elements in ARRAY,
 859    which has COUNT elements of SIZE bytes.  Returns the first
 860    element in ARRAY such that COMPARE returns zero when it and
 861    its successor element are compared, or a null pointer if no
 862    such element exists.  AUX is passed to COMPARE as auxiliary
 863    data. */
 864 void *
 865 adjacent_find_equal (const void *array, size_t count, size_t size,
 866                      algo_compare_func *compare, void *aux)
 867 {
 868   const unsigned char *first = array;
 869   const unsigned char *last = first + count * size;
 870
 871   while (first < last && first + size < last)
 872     {
 873       if (compare (first, first + size, aux) == 0)
 874         return (void *) first;
 875       first += size;
 876     }
 877
 878   return NULL;
 879 }
 880 \f
 881 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
 882    the first COUNT - 1 elements of these form a heap, followed by
 883    a single element not part of the heap.  This function adds the
 884    final element, forming a heap of COUNT elements in ARRAY.
 885    Uses COMPARE to compare elements, passing AUX as auxiliary
 886    data. */
 887 void
 888 push_heap (void *array, size_t count, size_t size,
 889            algo_compare_func *compare, void *aux)
 890 {
 891   unsigned char *first = array;
 892   size_t i;
 893
 894   expensive_assert (count < 1 || is_heap (array, count - 1,
 895                                           size, compare, aux));
 896   for (i = count; i > 1; i /= 2)
 897     {
 898       unsigned char *parent = first + (i / 2 - 1) * size;
 899       unsigned char *element = first + (i - 1) * size;
 900       if (compare (parent, element, aux) < 0)
 901         SWAP (parent, element, size);
 902       else
 903         break;
 904     }
 905   expensive_assert (is_heap (array, count, size, compare, aux));
 906 }
 907
 908 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
 909    the children of ARRAY[idx - 1] are heaps, but ARRAY[idx - 1]
 910    may be smaller than its children.  This function fixes that,
 911    so that ARRAY[idx - 1] itself is a heap.  Uses COMPARE to
 912    compare elements, passing AUX as auxiliary data. */
 913 static void
 914 heapify (void *array, size_t count, size_t size,
 915          size_t idx,
 916          algo_compare_func *compare, void *aux)
 917 {
 918   unsigned char *first = array;
 919
 920   for (;;)
 921     {
 922       size_t left = 2 * idx;
 923       size_t right = left + 1;
 924       size_t largest = idx;
 925
 926       if (left <= count
 927           && compare (first + size * (left - 1),
 928                       first + size * (idx - 1), aux) > 0)
 929         largest = left;
 930
 931       if (right <= count
 932           && compare (first + size * (right - 1),
 933                       first + size * (largest - 1), aux) > 0)
 934         largest = right;
 935
 936       if (largest == idx)
 937         break;
 938
 939       SWAP (first + size * (idx - 1), first + size * (largest - 1), size);
 940       idx = largest;
 941     }
 942 }
 943
 944 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
 945    all COUNT elements form a heap.  This function moves the
 946    largest element in the heap to the final position in ARRAY and
 947    reforms a heap of the remaining COUNT - 1 elements at the
 948    beginning of ARRAY.  Uses COMPARE to compare elements, passing
 949    AUX as auxiliary data. */
 950 void
 951 pop_heap (void *array, size_t count, size_t size,
 952           algo_compare_func *compare, void *aux)
 953 {
 954   unsigned char *first = array;
 955
 956   expensive_assert (is_heap (array, count, size, compare, aux));
 957   SWAP (first, first + (count - 1) * size, size);
 958   heapify (first, count - 1, size, 1, compare, aux);
 959   expensive_assert (count < 1 || is_heap (array, count - 1,
 960                                           size, compare, aux));
 961 }
 962
 963 /* Turns ARRAY, which contains COUNT elements of SIZE bytes, into
 964    a heap.  Uses COMPARE to compare elements, passing AUX as
 965    auxiliary data. */
 966 void
 967 make_heap (void *array, size_t count, size_t size,
 968            algo_compare_func *compare, void *aux)
 969 {
 970   size_t idx;
 971
 972   for (idx = count / 2; idx >= 1; idx--)
 973     heapify (array, count, size, idx, compare, aux);
 974   expensive_assert (count < 1 || is_heap (array, count, size, compare, aux));
 975 }
 976
 977 /* ARRAY contains COUNT elements of SIZE bytes each.  Initially
 978    all COUNT elements form a heap.  This function turns the heap
 979    into a fully sorted array.  Uses COMPARE to compare elements,
 980    passing AUX as auxiliary data. */
 981 void
 982 sort_heap (void *array, size_t count, size_t size,
 983            algo_compare_func *compare, void *aux)
 984 {
 985   unsigned char *first = array;
 986   size_t idx;
 987
 988   expensive_assert (is_heap (array, count, size, compare, aux));
 989   for (idx = count; idx >= 2; idx--)
 990     {
 991       SWAP (first, first + (idx - 1) * size, size);
 992       heapify (array, idx - 1, size, 1, compare, aux);
 993     }
 994   expensive_assert (is_sorted (array, count, size, compare, aux));
 995 }
 996
 997 /* ARRAY contains COUNT elements of SIZE bytes each.  This
 998    function tests whether ARRAY is a heap and returns 1 if so, 0
 999    otherwise.  Uses COMPARE to compare elements, passing AUX as
1000    auxiliary data. */
1001 int
1002 is_heap (const void *array, size_t count, size_t size,
1003          algo_compare_func *compare, void *aux)
1004 {
1005   const unsigned char *first = array;
1006   size_t child;
1007
1008   for (child = 2; child <= count; child++)
1009     {
1010       size_t parent = child / 2;
1011       if (compare (first + (parent - 1) * size,
1012                    first + (child - 1) * size, aux) < 0)
1013         return 0;
1014     }
1015
1016   return 1;
1017 }
1018