gl/lib/randperm.c

   1 /* Generate random permutations.
   2
   3    Copyright (C) 2006-2022 Free Software Foundation, Inc.
   4
   5    This program is free software: you can redistribute it and/or modify
   6    it under the terms of the GNU General Public License as published by
   7    the Free Software Foundation, either version 3 of the License, or
   8    (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13    GNU 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, see <https://www.gnu.org/licenses/>.  */
  17
  18 /* Written by Paul Eggert.  */
  19
  20 #include <config.h>
  21
  22 #include "randperm.h"
  23
  24 #include <limits.h>
  25 #include <stdint.h>
  26 #include <stdlib.h>
  27
  28 #include "attribute.h"
  29 #include "count-leading-zeros.h"
  30 #include "hash.h"
  31 #include "verify.h"
  32 #include "xalloc.h"
  33
  34 /* Return the floor of the log base 2 of N.  If N is zero, return -1.  */
  35
  36 ATTRIBUTE_CONST static int
  37 floor_lg (size_t n)
  38 {
  39   verify (SIZE_WIDTH <= ULLONG_WIDTH);
  40   return (n == 0 ? -1
  41           : SIZE_WIDTH <= UINT_WIDTH
  42           ? UINT_WIDTH - 1 - count_leading_zeros (n)
  43           : SIZE_WIDTH <= ULONG_WIDTH
  44           ? ULONG_WIDTH - 1 - count_leading_zeros_l (n)
  45           : ULLONG_WIDTH - 1 - count_leading_zeros_ll (n));
  46 }
  47
  48 /* Return an upper bound on the number of random bytes needed to
  49    generate the first H elements of a random permutation of N
  50    elements.  H must not exceed N.  */
  51
  52 size_t
  53 randperm_bound (size_t h, size_t n)
  54 {
  55   /* Upper bound on number of bits needed to generate the first number
  56      of the permutation.  */
  57   uintmax_t lg_n = floor_lg (n) + 1;
  58
  59   /* Upper bound on number of bits needed to generated the first H elements.  */
  60   uintmax_t ar = lg_n * h;
  61
  62   /* Convert the bit count to a byte count.  */
  63   size_t bound = (ar + CHAR_BIT - 1) / CHAR_BIT;
  64
  65   return bound;
  66 }
  67
  68 /* Swap elements I and J in array V.  */
  69
  70 static void
  71 swap (size_t *v, size_t i, size_t j)
  72 {
  73   size_t t = v[i];
  74   v[i] = v[j];
  75   v[j] = t;
  76 }
  77
  78 /* Structures and functions for a sparse_map abstract data type that's
  79    used to effectively swap elements I and J in array V like swap(),
  80    but in a more memory efficient manner (when the number of permutations
  81    performed is significantly less than the size of the input).  */
  82
  83 struct sparse_ent_
  84 {
  85    size_t index;
  86    size_t val;
  87 };
  88
  89 static size_t
  90 sparse_hash_ (void const *x, size_t table_size)
  91 {
  92   struct sparse_ent_ const *ent = x;
  93   return ent->index % table_size;
  94 }
  95
  96 static bool
  97 sparse_cmp_ (void const *x, void const *y)
  98 {
  99   struct sparse_ent_ const *ent1 = x;
 100   struct sparse_ent_ const *ent2 = y;
 101   return ent1->index == ent2->index;
 102 }
 103
 104 typedef Hash_table sparse_map;
 105
 106 /* Initialize the structure for the sparse map,
 107    when a best guess as to the number of entries
 108    specified with SIZE_HINT.  */
 109
 110 static sparse_map *
 111 sparse_new (size_t size_hint)
 112 {
 113   return hash_initialize (size_hint, NULL, sparse_hash_, sparse_cmp_, free);
 114 }
 115
 116 /* Swap the values for I and J.  If a value is not already present
 117    then assume it's equal to the index.  Update the value for
 118    index I in array V.  */
 119
 120 static void
 121 sparse_swap (sparse_map *sv, size_t *v, size_t i, size_t j)
 122 {
 123   struct sparse_ent_ *v1 = hash_remove (sv, &(struct sparse_ent_) {i,0});
 124   struct sparse_ent_ *v2 = hash_remove (sv, &(struct sparse_ent_) {j,0});
 125
 126   /* FIXME: reduce the frequency of these mallocs.  */
 127   if (!v1)
 128     {
 129       v1 = xmalloc (sizeof *v1);
 130       v1->index = v1->val = i;
 131     }
 132   if (!v2)
 133     {
 134       v2 = xmalloc (sizeof *v2);
 135       v2->index = v2->val = j;
 136     }
 137
 138   size_t t = v1->val;
 139   v1->val = v2->val;
 140   v2->val = t;
 141   if (!hash_insert (sv, v1))
 142     xalloc_die ();
 143   if (!hash_insert (sv, v2))
 144     xalloc_die ();
 145
 146   v[i] = v1->val;
 147 }
 148
 149 static void
 150 sparse_free (sparse_map *sv)
 151 {
 152   hash_free (sv);
 153 }
 154
 155
 156 /* From R, allocate and return a malloc'd array of the first H elements
 157    of a random permutation of N elements.  H must not exceed N.
 158    Return NULL if H is zero.  */
 159
 160 size_t *
 161 randperm_new (struct randint_source *r, size_t h, size_t n)
 162 {
 163   size_t *v;
 164
 165   switch (h)
 166     {
 167     case 0:
 168       v = NULL;
 169       break;
 170
 171     case 1:
 172       v = xmalloc (sizeof *v);
 173       v[0] = randint_choose (r, n);
 174       break;
 175
 176     default:
 177       {
 178         /* The algorithm is essentially the same in both
 179            the sparse and non sparse case.  In the sparse case we use
 180            a hash to implement sparse storage for the set of n numbers
 181            we're shuffling.  When to use the sparse method was
 182            determined with the help of this script:
 183
 184            #!/bin/sh
 185            for n in $(seq 2 32); do
 186              for h in $(seq 2 32); do
 187                test $h -gt $n && continue
 188                for s in o n; do
 189                  test $s = o && shuf=shuf || shuf=./shuf
 190                  num=$(env time -f "$s:${h},${n} = %e,%M" \
 191                        $shuf -i0-$((2**$n-2)) -n$((2**$h-2)) | wc -l)
 192                  test $num = $((2**$h-2)) || echo "$s:${h},${n} = failed" >&2
 193                done
 194              done
 195            done
 196
 197            This showed that if sparseness = n/h, then:
 198
 199            sparseness = 128 => .125 mem used, and about same speed
 200            sparseness =  64 => .25  mem used, but 1.5 times slower
 201            sparseness =  32 => .5   mem used, but 2 times slower
 202
 203            Also the memory usage was only significant when n > 128Ki
 204         */
 205         bool sparse = (n >= (128 * 1024)) && (n / h >= 32);
 206
 207         size_t i;
 208         sparse_map *sv;
 209
 210         if (sparse)
 211           {
 212             sv = sparse_new (h * 2);
 213             if (sv == NULL)
 214               xalloc_die ();
 215             v = xnmalloc (h, sizeof *v);
 216           }
 217         else
 218           {
 219             sv = NULL; /* To placate GCC's -Wuninitialized.  */
 220             v = xnmalloc (n, sizeof *v);
 221             for (i = 0; i < n; i++)
 222               v[i] = i;
 223           }
 224
 225         for (i = 0; i < h; i++)
 226           {
 227             size_t j = i + randint_choose (r, n - i);
 228             if (sparse)
 229               sparse_swap (sv, v, i, j);
 230             else
 231               swap (v, i, j);
 232           }
 233
 234         if (sparse)
 235           sparse_free (sv);
 236         else
 237           v = xnrealloc (v, h, sizeof *v);
 238       }
 239       break;
 240     }
 241
 242   return v;
 243 }