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