src/algorithm.c

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