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