src/algorithm.c

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