src/libpspp/array.c

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