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 #include <config.h>
  74 #include <assert.h>
  75 #include <limits.h>
  76 #include <stdlib.h>
  77 #include <string.h>
  78 #include "alloc.h"
  79 #include "algorithm.h"
  80 #include "random.h"
  81 \f
  82 /* Byte-wise swap two items of size SIZE. */
  83 #define SWAP(a, b, size)                                                      \
  84   do                                                                          \
  85     {                                                                         \
  86       register size_t __size = (size);                                        \
  87       register char *__a = (a), *__b = (b);                                   \
  88       do                                                                      \
  89         {                                                                     \
  90           char __tmp = *__a;                                                  \
  91           *__a++ = *__b;                                                      \
  92           *__b++ = __tmp;                                                     \
  93         } while (--__size > 0);                                               \
  94     } while (0)
  95
  96 /* Makes the elements in ARRAY unique, by moving up duplicates,
  97    and returns the new number of elements in the array.  Sorted
  98    arrays only.  Arguments same as for sort() above. */
  99 size_t
 100 unique (void *array, size_t count, size_t size,
 101         algo_compare_func *compare, void *aux)
 102 {
 103   char *first = array;
 104   char *last = first + size * count;
 105   char *result = array;
 106
 107   for (;;)
 108     {
 109       first += size;
 110       if (first >= last)
 111         return count;
 112
 113       if (compare (result, first, aux))
 114         {
 115           result += size;
 116           if (result != first)
 117             memcpy (result, first, size);
 118         }
 119       else
 120         count--;
 121     }
 122 }
 123
 124 /* Helper function that calls sort(), then unique(). */
 125 size_t
 126 sort_unique (void *array, size_t count, size_t size,
 127              algo_compare_func *compare, void *aux)
 128 {
 129   sort (array, count, size, compare, aux);
 130   return unique (array, count, size, compare, aux);
 131 }
 132
 133 #ifdef TEST_UNIQUE
 134 #include <stdio.h>
 135
 136 void *
 137 xmalloc (size_t size)
 138 {
 139   return malloc (size);
 140 }
 141
 142 int
 143 compare_ints (const void *a_, const void *b_, void *aux)
 144 {
 145   const int *a = a_;
 146   const int *b = b_;
 147
 148   if (*a > *b)
 149     return 1;
 150   else if (*a < *b)
 151     return 1;
 152   else
 153     return 0;
 154 }
 155
 156 void
 157 try_unique (const char *title,
 158             int *in, size_t in_cnt,
 159             size_t out_cnt)
 160 {
 161   size_t i;
 162
 163   in_cnt = unique (in, in_cnt, sizeof *in, compare_ints, NULL);
 164   if (in_cnt != out_cnt)
 165     {
 166       fprintf (stderr, "unique_test: %s: in_cnt %d, expected %d\n",
 167                title, (int) in_cnt, (int) out_cnt);
 168       return;
 169     }
 170
 171   for (i = 0; i < out_cnt; i++)
 172     {
 173       if (in[i] != i)
 174         fprintf (stderr, "unique_test: %s: idx %d = %d, expected %d\n",
 175                  title, (int) i, in[i], i);
 176     }
 177 }
 178
 179 int
 180 main (void)
 181 {
 182   int a_in[] = {0, 0, 0, 1, 2, 3, 3, 4, 5, 5};
 183   int b_in[] = {0, 1, 2, 2, 2, 3};
 184   int c_in[] = {0};
 185   int d_in;
 186
 187   try_unique ("a", a_in, sizeof a_in / sizeof *a_in, 6);
 188   try_unique ("b", b_in, sizeof b_in / sizeof *b_in, 4);
 189   try_unique ("c", c_in, sizeof c_in / sizeof *c_in, 1);
 190   try_unique ("d", &d_in, 0, 0);
 191
 192 }
 193 #endif /* TEST_UNIQUE */
 194 \f
 195 /* Reorders ARRAY, which contains COUNT elements of SIZE bytes
 196    each, so that the elements for which PREDICATE returns nonzero
 197    precede those for which PREDICATE returns zero.  AUX is
 198    passed to each predicate as auxiliary data.  Returns the
 199    number of elements for which PREDICATE returns nonzero.  Not
 200    stable. */
 201 size_t
 202 partition (void *array, size_t count, size_t size,
 203            algo_predicate_func *predicate, void *aux)
 204 {
 205   char *first = array;
 206   char *last = first + count * size;
 207
 208   for (;;)
 209     {
 210       /* Move FIRST forward to point to first element that fails
 211          PREDICATE. */
 212       for (;;)
 213         {
 214           if (first == last)
 215             return count;
 216           else if (!predicate (first, aux))
 217             break;
 218
 219           first += size;
 220         }
 221       count--;
 222
 223       /* Move LAST backward to point to last element that passes
 224          PREDICATE. */
 225       for (;;)
 226         {
 227           last -= size;
 228
 229           if (first == last)
 230             return count;
 231           else if (predicate (last, aux))
 232             break;
 233           else
 234             count--;
 235         }
 236
 237       /* By swapping FIRST and LAST we extend the starting and
 238          ending sequences that pass and fail, respectively,
 239          PREDICATE. */
 240       SWAP (first, last, size);
 241       first += size;
 242     }
 243 }
 244 \f
 245 /* A algo_random_func that uses random.h. */
 246 unsigned
 247 algo_default_random (unsigned max, void *aux unused)
 248 {
 249   return rng_get_unsigned (pspp_rng ()) % max;
 250 }
 251
 252 /* Randomly reorders ARRAY, which contains COUNT elements of SIZE
 253    bytes each.  Uses RANDOM as a source of random data, passing
 254    AUX as the auxiliary data.  RANDOM may be null to use a
 255    default random source. */
 256 void
 257 random_shuffle (void *array_, size_t count, size_t size,
 258                 algo_random_func *random, void *aux)
 259 {
 260   unsigned char *array = array_;
 261   int i;
 262
 263   if (random == NULL)
 264     random = algo_default_random;
 265
 266   for (i = 1; i < count; i++)
 267     SWAP (array + i * size, array + random (i + 1, aux) * size, size);
 268 }
 269 \f
 270 /* Copies the COUNT elements of SIZE bytes each from ARRAY to
 271    RESULT, except that elements for which PREDICATE is false are
 272    not copied.  Returns the number of elements copied.  AUX is
 273    passed to PREDICATE as auxiliary data.  */
 274 size_t
 275 copy_if (const void *array, size_t count, size_t size,
 276          void *result,
 277          algo_predicate_func *predicate, void *aux)
 278 {
 279   const unsigned char *input = array;
 280   const unsigned char *last = input + size * count;
 281   unsigned char *output = result;
 282
 283   while (input <= last)
 284     {
 285       if (predicate (input, aux))
 286         {
 287           memcpy (output, input, size);
 288           output += size;
 289         }
 290       else
 291         count--;
 292
 293       input += size;
 294     }
 295
 296   return count;
 297 }
 298
 299 /* A predicate and its auxiliary data. */
 300 struct pred_aux
 301   {
 302     algo_predicate_func *predicate;
 303     void *aux;
 304   };
 305
 306 static int
 307 not (const void *data, void *pred_aux_)
 308 {
 309   const struct pred_aux *pred_aux = pred_aux_;
 310
 311   return !pred_aux->predicate (data, pred_aux->aux);
 312 }
 313
 314 /* Copies the COUNT elements of SIZE bytes each from ARRAY to
 315    RESULT, except that elements for which PREDICATE is true are
 316    not copied.  Returns the number of elements copied.  AUX is
 317    passed to PREDICATE as auxiliary data.  */
 318 size_t
 319 remove_copy_if (const void *array, size_t count, size_t size,
 320                 void *result,
 321                 algo_predicate_func *predicate, void *aux)
 322 {
 323   struct pred_aux pred_aux;
 324   pred_aux.predicate = predicate;
 325   pred_aux.aux = aux;
 326   return copy_if (array, count, size, result, not, &pred_aux);
 327 }
 328 \f
 329 /* Searches ARRAY, which contains COUNT of SIZE bytes each, using
 330    a binary search.  Returns any element that equals VALUE, if
 331    one exists, or a null pointer otherwise.  ARRAY must ordered
 332    according to COMPARE.  AUX is passed to COMPARE as auxiliary
 333    data. */
 334 void *
 335 binary_search (const void *array, size_t count, size_t size,
 336                void *value,
 337                algo_compare_func *compare, void *aux)
 338 {
 339   assert (array != NULL);
 340   assert (count <= INT_MAX);
 341   assert (compare != NULL);
 342
 343   if (count != 0)
 344     {
 345       const unsigned char *first = array;
 346       int low = 0;
 347       int high = count - 1;
 348
 349       while (low < high)
 350         {
 351           int middle = (low + high) / 2;
 352           const unsigned char *element = first + middle * size;
 353           int cmp = compare (value, element, aux);
 354
 355           if (cmp > 0)
 356             low = middle + 1;
 357           else if (cmp < 0)
 358             high = middle - 1;
 359           else
 360             return (void *) element;
 361         }
 362     }
 363   return NULL;
 364 }
 365 \f
 366 /* Lexicographically compares ARRAY1, which contains COUNT1
 367    elements of SIZE bytes each, to ARRAY2, which contains COUNT2
 368    elements of SIZE bytes, according to COMPARE.  Returns a
 369    strcmp()-type result.  AUX is passed to COMPARE as auxiliary
 370    data. */
 371 int
 372 lexicographical_compare (const void *array1, size_t count1,
 373                          const void *array2, size_t count2,
 374                          size_t size,
 375                          algo_compare_func *compare, void *aux)
 376 {
 377   const unsigned char *first1 = array1;
 378   const unsigned char *first2 = array2;
 379   size_t min_count = count1 < count2 ? count1 : count2;
 380
 381   while (min_count > 0)
 382     {
 383       int cmp = compare (first1, first2, aux);
 384       if (cmp != 0)
 385         return cmp;
 386
 387       first1 += size;
 388       first2 += size;
 389       min_count--;
 390     }
 391
 392   return count1 < count2 ? -1 : count1 > count2;
 393 }
 394
 395 \f
 396 /* Copyright (C) 1991, 1992, 1996, 1997, 1999 Free Software Foundation, Inc.
 397    This file is part of the GNU C Library.
 398    Written by Douglas C. Schmidt (schmidt@ics.uci.edu).
 399
 400    The GNU C Library is free software; you can redistribute it and/or
 401    modify it under the terms of the GNU Lesser General Public
 402    License as published by the Free Software Foundation; either
 403    version 2.1 of the License, or (at your option) any later version.
 404
 405    The GNU C Library is distributed in the hope that it will be useful,
 406    but WITHOUT ANY WARRANTY; without even the implied warranty of
 407    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 408    Lesser General Public License for more details.
 409
 410    You should have received a copy of the GNU Lesser General Public
 411    License along with the GNU C Library; if not, write to the Free
 412    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 413    02111-1307 USA.  */
 414
 415 /* If you consider tuning this algorithm, you should consult first:
 416    Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
 417    Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */
 418
 419 #include <alloca.h>
 420 #include <limits.h>
 421 #include <stdlib.h>
 422 #include <string.h>
 423
 424 /* Discontinue quicksort algorithm when partition gets below this size.
 425    This particular magic number was chosen to work best on a Sun 4/260. */
 426 #define MAX_THRESH 4
 427
 428 /* Stack node declarations used to store unfulfilled partition obligations. */
 429 typedef struct
 430   {
 431     char *lo;
 432     char *hi;
 433   } stack_node;
 434
 435 /* The next 4 #defines implement a very fast in-line stack abstraction. */
 436 /* The stack needs log (total_elements) entries (we could even subtract
 437    log(MAX_THRESH)).  Since total_elements has type size_t, we get as
 438    upper bound for log (total_elements):
 439    bits per byte (CHAR_BIT) * sizeof(size_t).  */
 440 #define STACK_SIZE      (CHAR_BIT * sizeof(size_t))
 441 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
 442 #define POP(low, high)  ((void) (--top, (low = top->lo), (high = top->hi)))
 443 #define STACK_NOT_EMPTY (stack < top)
 444
 445
 446 /* Order size using quicksort.  This implementation incorporates
 447    four optimizations discussed in Sedgewick:
 448
 449    1. Non-recursive, using an explicit stack of pointer that store the
 450       next array partition to sort.  To save time, this maximum amount
 451       of space required to store an array of SIZE_MAX is allocated on the
 452       stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
 453       only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
 454       Pretty cheap, actually.
 455
 456    2. Chose the pivot element using a median-of-three decision tree.
 457       This reduces the probability of selecting a bad pivot value and
 458       eliminates certain extraneous comparisons.
 459
 460    3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
 461       insertion sort to order the MAX_THRESH items within each partition.
 462       This is a big win, since insertion sort is faster for small, mostly
 463       sorted array segments.
 464
 465    4. The larger of the two sub-partitions is always pushed onto the
 466       stack first, with the algorithm then concentrating on the
 467       smaller partition.  This *guarantees* no more than log (total_elems)
 468       stack size is needed (actually O(1) in this case)!  */
 469
 470 void
 471 sort (void *const pbase, size_t total_elems, size_t size,
 472       algo_compare_func *cmp, void *aux)
 473 {
 474   register char *base_ptr = (char *) pbase;
 475
 476   const size_t max_thresh = MAX_THRESH * size;
 477
 478   if (total_elems == 0)
 479     /* Avoid lossage with unsigned arithmetic below.  */
 480     return;
 481
 482   if (total_elems > MAX_THRESH)
 483     {
 484       char *lo = base_ptr;
 485       char *hi = &lo[size * (total_elems - 1)];
 486       stack_node stack[STACK_SIZE];
 487       stack_node *top = stack + 1;
 488
 489       while (STACK_NOT_EMPTY)
 490         {
 491           char *left_ptr;
 492           char *right_ptr;
 493
 494           /* Select median value from among LO, MID, and HI. Rearrange
 495              LO and HI so the three values are sorted. This lowers the
 496              probability of picking a pathological pivot value and
 497              skips a comparison for both the LEFT_PTR and RIGHT_PTR in
 498              the while loops. */
 499
 500           char *mid = lo + size * ((hi - lo) / size >> 1);
 501
 502           if ((*cmp) ((void *) mid, (void *) lo, aux) < 0)
 503             SWAP (mid, lo, size);
 504           if ((*cmp) ((void *) hi, (void *) mid, aux) < 0)
 505             SWAP (mid, hi, size);
 506           else
 507             goto jump_over;
 508           if ((*cmp) ((void *) mid, (void *) lo, aux) < 0)
 509             SWAP (mid, lo, size);
 510         jump_over:;
 511
 512           left_ptr  = lo + size;
 513           right_ptr = hi - size;
 514
 515           /* Here's the famous ``collapse the walls'' section of quicksort.
 516              Gotta like those tight inner loops!  They are the main reason
 517              that this algorithm runs much faster than others. */
 518           do
 519             {
 520               while ((*cmp) ((void *) left_ptr, (void *) mid, aux) < 0)
 521                 left_ptr += size;
 522
 523               while ((*cmp) ((void *) mid, (void *) right_ptr, aux) < 0)
 524                 right_ptr -= size;
 525
 526               if (left_ptr < right_ptr)
 527                 {
 528                   SWAP (left_ptr, right_ptr, size);
 529                   if (mid == left_ptr)
 530                     mid = right_ptr;
 531                   else if (mid == right_ptr)
 532                     mid = left_ptr;
 533                   left_ptr += size;
 534                   right_ptr -= size;
 535                 }
 536               else if (left_ptr == right_ptr)
 537                 {
 538                   left_ptr += size;
 539                   right_ptr -= size;
 540                   break;
 541                 }
 542             }
 543           while (left_ptr <= right_ptr);
 544
 545           /* Set up pointers for next iteration.  First determine whether
 546              left and right partitions are below the threshold size.  If so,
 547              ignore one or both.  Otherwise, push the larger partition's
 548              bounds on the stack and continue sorting the smaller one. */
 549
 550           if ((size_t) (right_ptr - lo) <= max_thresh)
 551             {
 552               if ((size_t) (hi - left_ptr) <= max_thresh)
 553                 /* Ignore both small partitions. */
 554                 POP (lo, hi);
 555               else
 556                 /* Ignore small left partition. */
 557                 lo = left_ptr;
 558             }
 559           else if ((size_t) (hi - left_ptr) <= max_thresh)
 560             /* Ignore small right partition. */
 561             hi = right_ptr;
 562           else if ((right_ptr - lo) > (hi - left_ptr))
 563             {
 564               /* Push larger left partition indices. */
 565               PUSH (lo, right_ptr);
 566               lo = left_ptr;
 567             }
 568           else
 569             {
 570               /* Push larger right partition indices. */
 571               PUSH (left_ptr, hi);
 572               hi = right_ptr;
 573             }
 574         }
 575     }
 576
 577   /* Once the BASE_PTR array is partially sorted by quicksort the rest
 578      is completely sorted using insertion sort, since this is efficient
 579      for partitions below MAX_THRESH size. BASE_PTR points to the beginning
 580      of the array to sort, and END_PTR points at the very last element in
 581      the array (*not* one beyond it!). */
 582
 583 #define min(x, y) ((x) < (y) ? (x) : (y))
 584
 585   {
 586     char *const end_ptr = &base_ptr[size * (total_elems - 1)];
 587     char *tmp_ptr = base_ptr;
 588     char *thresh = min(end_ptr, base_ptr + max_thresh);
 589     register char *run_ptr;
 590
 591     /* Find smallest element in first threshold and place it at the
 592        array's beginning.  This is the smallest array element,
 593        and the operation speeds up insertion sort's inner loop. */
 594
 595     for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
 596       if ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, aux) < 0)
 597         tmp_ptr = run_ptr;
 598
 599     if (tmp_ptr != base_ptr)
 600       SWAP (tmp_ptr, base_ptr, size);
 601
 602     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
 603
 604     run_ptr = base_ptr + size;
 605     while ((run_ptr += size) <= end_ptr)
 606       {
 607         tmp_ptr = run_ptr - size;
 608         while ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, aux) < 0)
 609           tmp_ptr -= size;
 610
 611         tmp_ptr += size;
 612         if (tmp_ptr != run_ptr)
 613           {
 614             char *trav;
 615
 616             trav = run_ptr + size;
 617             while (--trav >= run_ptr)
 618               {
 619                 char c = *trav;
 620                 char *hi, *lo;
 621
 622                 for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
 623                   *hi = *lo;
 624                 *hi = c;
 625               }
 626           }
 627       }
 628   }
 629 }