lib/random_r.c

   1 /*
   2    Copyright (C) 1995, 2005, 2008-2010 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  16
  17 /*
  18    Copyright (C) 1983 Regents of the University of California.
  19    All rights reserved.
  20
  21    Redistribution and use in source and binary forms, with or without
  22    modification, are permitted provided that the following conditions
  23    are met:
  24
  25    1. Redistributions of source code must retain the above copyright
  26       notice, this list of conditions and the following disclaimer.
  27    2. Redistributions in binary form must reproduce the above copyright
  28       notice, this list of conditions and the following disclaimer in the
  29       documentation and/or other materials provided with the distribution.
  30    4. Neither the name of the University nor the names of its contributors
  31       may be used to endorse or promote products derived from this software
  32       without specific prior written permission.
  33
  34    THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  35    ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  36    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  37    ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  38    FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  39    DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  40    OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  41    HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  42    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  43    OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  44    SUCH DAMAGE.*/
  45
  46 /*
  47  * This is derived from the Berkeley source:
  48  *      @(#)random.c    5.5 (Berkeley) 7/6/88
  49  * It was reworked for the GNU C Library by Roland McGrath.
  50  * Rewritten to be reentrant by Ulrich Drepper, 1995
  51  */
  52
  53 #include <config.h>
  54
  55 /* Don't use __attribute__ __nonnull__ in this compilation unit.  Otherwise gcc
  56    optimizes away the buf == NULL, arg_state == NULL, result == NULL tests
  57    below.  */
  58 #define _GL_ARG_NONNULL(params)
  59
  60 #include <errno.h>
  61 #include <limits.h>
  62 #include <stddef.h>
  63 #include <stdlib.h>
  64 #include <inttypes.h>
  65
  66
  67 /* An improved random number generation package.  In addition to the standard
  68    rand()/srand() like interface, this package also has a special state info
  69    interface.  The initstate() routine is called with a seed, an array of
  70    bytes, and a count of how many bytes are being passed in; this array is
  71    then initialized to contain information for random number generation with
  72    that much state information.  Good sizes for the amount of state
  73    information are 32, 64, 128, and 256 bytes.  The state can be switched by
  74    calling the setstate() function with the same array as was initialized
  75    with initstate().  By default, the package runs with 128 bytes of state
  76    information and generates far better random numbers than a linear
  77    congruential generator.  If the amount of state information is less than
  78    32 bytes, a simple linear congruential R.N.G. is used.  Internally, the
  79    state information is treated as an array of longs; the zeroth element of
  80    the array is the type of R.N.G. being used (small integer); the remainder
  81    of the array is the state information for the R.N.G.  Thus, 32 bytes of
  82    state information will give 7 longs worth of state information, which will
  83    allow a degree seven polynomial.  (Note: The zeroth word of state
  84    information also has some other information stored in it; see setstate
  85    for details).  The random number generation technique is a linear feedback
  86    shift register approach, employing trinomials (since there are fewer terms
  87    to sum up that way).  In this approach, the least significant bit of all
  88    the numbers in the state table will act as a linear feedback shift register,
  89    and will have period 2^deg - 1 (where deg is the degree of the polynomial
  90    being used, assuming that the polynomial is irreducible and primitive).
  91    The higher order bits will have longer periods, since their values are
  92    also influenced by pseudo-random carries out of the lower bits.  The
  93    total period of the generator is approximately deg*(2**deg - 1); thus
  94    doubling the amount of state information has a vast influence on the
  95    period of the generator.  Note: The deg*(2**deg - 1) is an approximation
  96    only good for large deg, when the period of the shift register is the
  97    dominant factor.  With deg equal to seven, the period is actually much
  98    longer than the 7*(2**7 - 1) predicted by this formula.  */
  99
 100
 101
 102 /* For each of the currently supported random number generators, we have a
 103    break value on the amount of state information (you need at least this many
 104    bytes of state info to support this random number generator), a degree for
 105    the polynomial (actually a trinomial) that the R.N.G. is based on, and
 106    separation between the two lower order coefficients of the trinomial.  */
 107
 108 /* Linear congruential.  */
 109 #define TYPE_0          0
 110 #define BREAK_0         8
 111 #define DEG_0           0
 112 #define SEP_0           0
 113
 114 /* x**7 + x**3 + 1.  */
 115 #define TYPE_1          1
 116 #define BREAK_1         32
 117 #define DEG_1           7
 118 #define SEP_1           3
 119
 120 /* x**15 + x + 1.  */
 121 #define TYPE_2          2
 122 #define BREAK_2         64
 123 #define DEG_2           15
 124 #define SEP_2           1
 125
 126 /* x**31 + x**3 + 1.  */
 127 #define TYPE_3          3
 128 #define BREAK_3         128
 129 #define DEG_3           31
 130 #define SEP_3           3
 131
 132 /* x**63 + x + 1.  */
 133 #define TYPE_4          4
 134 #define BREAK_4         256
 135 #define DEG_4           63
 136 #define SEP_4           1
 137
 138
 139 /* Array versions of the above information to make code run faster.
 140    Relies on fact that TYPE_i == i.  */
 141
 142 #define MAX_TYPES       5       /* Max number of types above.  */
 143
 144 struct random_poly_info
 145 {
 146   int seps[MAX_TYPES];
 147   int degrees[MAX_TYPES];
 148 };
 149
 150 static const struct random_poly_info random_poly_info =
 151 {
 152   { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 },
 153   { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 }
 154 };
 155
 156 #ifndef _LIBC
 157 # define weak_alias(local, symbol)
 158 # define __set_errno(e) errno = (e)
 159 # define __srandom_r srandom_r
 160 # define __initstate_r initstate_r
 161 # define __setstate_r setstate_r
 162 # define __random_r random_r
 163 #endif
 164
 165
 166 \f
 167 /* Initialize the random number generator based on the given seed.  If the
 168    type is the trivial no-state-information type, just remember the seed.
 169    Otherwise, initializes state[] based on the given "seed" via a linear
 170    congruential generator.  Then, the pointers are set to known locations
 171    that are exactly rand_sep places apart.  Lastly, it cycles the state
 172    information a given number of times to get rid of any initial dependencies
 173    introduced by the L.C.R.N.G.  Note that the initialization of randtbl[]
 174    for default usage relies on values produced by this routine.  */
 175 int
 176 __srandom_r (unsigned int seed, struct random_data *buf)
 177 {
 178   int type;
 179   int32_t *state;
 180   long int i;
 181   long int word;
 182   int32_t *dst;
 183   int kc;
 184
 185   if (buf == NULL)
 186     goto fail;
 187   type = buf->rand_type;
 188   if ((unsigned int) type >= MAX_TYPES)
 189     goto fail;
 190
 191   state = buf->state;
 192   /* We must make sure the seed is not 0.  Take arbitrarily 1 in this case.  */
 193   if (seed == 0)
 194     seed = 1;
 195   state[0] = seed;
 196   if (type == TYPE_0)
 197     goto done;
 198
 199   dst = state;
 200   word = seed;
 201   kc = buf->rand_deg;
 202   for (i = 1; i < kc; ++i)
 203     {
 204       /* This does:
 205            state[i] = (16807 * state[i - 1]) % 2147483647;
 206          but avoids overflowing 31 bits.  */
 207       long int hi = word / 127773;
 208       long int lo = word % 127773;
 209       word = 16807 * lo - 2836 * hi;
 210       if (word < 0)
 211         word += 2147483647;
 212       *++dst = word;
 213     }
 214
 215   buf->fptr = &state[buf->rand_sep];
 216   buf->rptr = &state[0];
 217   kc *= 10;
 218   while (--kc >= 0)
 219     {
 220       int32_t discard;
 221       (void) __random_r (buf, &discard);
 222     }
 223
 224  done:
 225   return 0;
 226
 227  fail:
 228   return -1;
 229 }
 230
 231 weak_alias (__srandom_r, srandom_r)
 232 \f
 233 /* Initialize the state information in the given array of N bytes for
 234    future random number generation.  Based on the number of bytes we
 235    are given, and the break values for the different R.N.G.'s, we choose
 236    the best (largest) one we can and set things up for it.  srandom is
 237    then called to initialize the state information.  Note that on return
 238    from srandom, we set state[-1] to be the type multiplexed with the current
 239    value of the rear pointer; this is so successive calls to initstate won't
 240    lose this information and will be able to restart with setstate.
 241    Note: The first thing we do is save the current state, if any, just like
 242    setstate so that it doesn't matter when initstate is called.
 243    Returns a pointer to the old state.  */
 244 int
 245 __initstate_r (unsigned int seed, char *arg_state, size_t n,
 246                struct random_data *buf)
 247 {
 248   int32_t *old_state;
 249   int32_t *state;
 250   int type;
 251   int degree;
 252   int separation;
 253
 254   if (buf == NULL)
 255     goto fail;
 256
 257   old_state = buf->state;
 258   if (old_state != NULL)
 259     {
 260       int old_type = buf->rand_type;
 261       if (old_type == TYPE_0)
 262         old_state[-1] = TYPE_0;
 263       else
 264         old_state[-1] = (MAX_TYPES * (buf->rptr - old_state)) + old_type;
 265     }
 266
 267   if (n >= BREAK_3)
 268     type = n < BREAK_4 ? TYPE_3 : TYPE_4;
 269   else if (n < BREAK_1)
 270     {
 271       if (n < BREAK_0)
 272         {
 273           __set_errno (EINVAL);
 274           goto fail;
 275         }
 276       type = TYPE_0;
 277     }
 278   else
 279     type = n < BREAK_2 ? TYPE_1 : TYPE_2;
 280
 281   degree = random_poly_info.degrees[type];
 282   separation = random_poly_info.seps[type];
 283
 284   buf->rand_type = type;
 285   buf->rand_sep = separation;
 286   buf->rand_deg = degree;
 287   state = &((int32_t *) arg_state)[1];  /* First location.  */
 288   /* Must set END_PTR before srandom.  */
 289   buf->end_ptr = &state[degree];
 290
 291   buf->state = state;
 292
 293   __srandom_r (seed, buf);
 294
 295   state[-1] = TYPE_0;
 296   if (type != TYPE_0)
 297     state[-1] = (buf->rptr - state) * MAX_TYPES + type;
 298
 299   return 0;
 300
 301  fail:
 302   __set_errno (EINVAL);
 303   return -1;
 304 }
 305
 306 weak_alias (__initstate_r, initstate_r)
 307 \f
 308 /* Restore the state from the given state array.
 309    Note: It is important that we also remember the locations of the pointers
 310    in the current state information, and restore the locations of the pointers
 311    from the old state information.  This is done by multiplexing the pointer
 312    location into the zeroth word of the state information. Note that due
 313    to the order in which things are done, it is OK to call setstate with the
 314    same state as the current state
 315    Returns a pointer to the old state information.  */
 316 int
 317 __setstate_r (char *arg_state, struct random_data *buf)
 318 {
 319   int32_t *new_state = 1 + (int32_t *) arg_state;
 320   int type;
 321   int old_type;
 322   int32_t *old_state;
 323   int degree;
 324   int separation;
 325
 326   if (arg_state == NULL || buf == NULL)
 327     goto fail;
 328
 329   old_type = buf->rand_type;
 330   old_state = buf->state;
 331   if (old_type == TYPE_0)
 332     old_state[-1] = TYPE_0;
 333   else
 334     old_state[-1] = (MAX_TYPES * (buf->rptr - old_state)) + old_type;
 335
 336   type = new_state[-1] % MAX_TYPES;
 337   if (type < TYPE_0 || type > TYPE_4)
 338     goto fail;
 339
 340   buf->rand_deg = degree = random_poly_info.degrees[type];
 341   buf->rand_sep = separation = random_poly_info.seps[type];
 342   buf->rand_type = type;
 343
 344   if (type != TYPE_0)
 345     {
 346       int rear = new_state[-1] / MAX_TYPES;
 347       buf->rptr = &new_state[rear];
 348       buf->fptr = &new_state[(rear + separation) % degree];
 349     }
 350   buf->state = new_state;
 351   /* Set end_ptr too.  */
 352   buf->end_ptr = &new_state[degree];
 353
 354   return 0;
 355
 356  fail:
 357   __set_errno (EINVAL);
 358   return -1;
 359 }
 360
 361 weak_alias (__setstate_r, setstate_r)
 362 \f
 363 /* If we are using the trivial TYPE_0 R.N.G., just do the old linear
 364    congruential bit.  Otherwise, we do our fancy trinomial stuff, which is the
 365    same in all the other cases due to all the global variables that have been
 366    set up.  The basic operation is to add the number at the rear pointer into
 367    the one at the front pointer.  Then both pointers are advanced to the next
 368    location cyclically in the table.  The value returned is the sum generated,
 369    reduced to 31 bits by throwing away the "least random" low bit.
 370    Note: The code takes advantage of the fact that both the front and
 371    rear pointers can't wrap on the same call by not testing the rear
 372    pointer if the front one has wrapped.  Returns a 31-bit random number.  */
 373
 374 int
 375 __random_r (struct random_data *buf, int32_t *result)
 376 {
 377   int32_t *state;
 378
 379   if (buf == NULL || result == NULL)
 380     goto fail;
 381
 382   state = buf->state;
 383
 384   if (buf->rand_type == TYPE_0)
 385     {
 386       int32_t val = state[0];
 387       val = ((state[0] * 1103515245) + 12345) & 0x7fffffff;
 388       state[0] = val;
 389       *result = val;
 390     }
 391   else
 392     {
 393       int32_t *fptr = buf->fptr;
 394       int32_t *rptr = buf->rptr;
 395       int32_t *end_ptr = buf->end_ptr;
 396       int32_t val;
 397
 398       val = *fptr += *rptr;
 399       /* Chucking least random bit.  */
 400       *result = (val >> 1) & 0x7fffffff;
 401       ++fptr;
 402       if (fptr >= end_ptr)
 403         {
 404           fptr = state;
 405           ++rptr;
 406         }
 407       else
 408         {
 409           ++rptr;
 410           if (rptr >= end_ptr)
 411             rptr = state;
 412         }
 413       buf->fptr = fptr;
 414       buf->rptr = rptr;
 415     }
 416   return 0;
 417
 418  fail:
 419   __set_errno (EINVAL);
 420   return -1;
 421 }
 422
 423 weak_alias (__random_r, random_r)