src/libpspp/array.c

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