lib/randperm.c

   1 /* Generate random permutations.
   2
   3    Copyright (C) 2006, 2007 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 <http://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
  26 #include "xalloc.h"
  27
  28 /* Return the ceiling of the log base 2 of N.  If N is zero, return
  29    an unspecified value.  */
  30
  31 static size_t
  32 ceil_lg (size_t n)
  33 {
  34   size_t b = 0;
  35   for (n--; n != 0; n /= 2)
  36     b++;
  37   return b;
  38 }
  39
  40 /* Return an upper bound on the number of random bytes needed to
  41    generate the first H elements of a random permutation of N
  42    elements.  H must not exceed N.  */
  43
  44 size_t
  45 randperm_bound (size_t h, size_t n)
  46 {
  47   /* Upper bound on number of bits needed to generate the first number
  48      of the permutation.  */
  49   size_t lg_n = ceil_lg (n);
  50
  51   /* Upper bound on number of bits needed to generated the first H elements.  */
  52   size_t ar = lg_n * h;
  53
  54   /* Convert the bit count to a byte count.  */
  55   size_t bound = (ar + CHAR_BIT - 1) / CHAR_BIT;
  56
  57   return bound;
  58 }
  59
  60 /* From R, allocate and return a malloc'd array of the first H elements
  61    of a random permutation of N elements.  H must not exceed N.
  62    Return NULL if H is zero.  */
  63
  64 size_t *
  65 randperm_new (struct randint_source *r, size_t h, size_t n)
  66 {
  67   size_t *v;
  68
  69   switch (h)
  70     {
  71     case 0:
  72       v = NULL;
  73       break;
  74
  75     case 1:
  76       v = xmalloc (sizeof *v);
  77       v[0] = randint_choose (r, n);
  78       break;
  79
  80     default:
  81       {
  82         size_t i;
  83
  84         v = xnmalloc (n, sizeof *v);
  85         for (i = 0; i < n; i++)
  86           v[i] = i;
  87
  88         for (i = 0; i < h; i++)
  89           {
  90             size_t j = i + randint_choose (r, n - i);
  91             size_t t = v[i];
  92             v[i] = v[j];
  93             v[j] = t;
  94           }
  95
  96         v = xnrealloc (v, h, sizeof *v);
  97       }
  98       break;
  99     }
 100
 101   return v;
 102 }