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