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