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