| 1 | /*
|
|---|
| 2 | * Copyright (c) 1983 Regents of the University of California.
|
|---|
| 3 | * All rights reserved.
|
|---|
| 4 | *
|
|---|
| 5 | * Redistribution and use in source and binary forms, with or without
|
|---|
| 6 | * modification, are permitted provided that the following conditions
|
|---|
| 7 | * are met:
|
|---|
| 8 | * 1. Redistributions of source code must retain the above copyright
|
|---|
| 9 | * notice, this list of conditions and the following disclaimer.
|
|---|
| 10 | * 2. Redistributions in binary form must reproduce the above copyright
|
|---|
| 11 | * notice, this list of conditions and the following disclaimer in the
|
|---|
| 12 | * documentation and/or other materials provided with the distribution.
|
|---|
| 13 | * 3. [rescinded 22 July 1999]
|
|---|
| 14 | * 4. Neither the name of the University nor the names of its contributors
|
|---|
| 15 | * may be used to endorse or promote products derived from this software
|
|---|
| 16 | * without specific prior written permission.
|
|---|
| 17 | *
|
|---|
| 18 | * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
|
|---|
| 19 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|---|
| 20 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|---|
| 21 | * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
|
|---|
| 22 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
|---|
| 23 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
|---|
| 24 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
|---|
| 25 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|---|
| 26 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
|---|
| 27 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
|---|
| 28 | * SUCH DAMAGE.
|
|---|
| 29 | */
|
|---|
| 30 |
|
|---|
| 31 | /*
|
|---|
| 32 | * This is derived from the Berkeley source:
|
|---|
| 33 | * @(#)random.c 5.5 (Berkeley) 7/6/88
|
|---|
| 34 | * It was reworked for the GNU C Library by Roland McGrath.
|
|---|
| 35 | */
|
|---|
| 36 |
|
|---|
| 37 | #include <errno.h>
|
|---|
| 38 |
|
|---|
| 39 | #if 0
|
|---|
| 40 |
|
|---|
| 41 | #include <ansidecl.h>
|
|---|
| 42 | #include <limits.h>
|
|---|
| 43 | #include <stddef.h>
|
|---|
| 44 | #include <stdlib.h>
|
|---|
| 45 |
|
|---|
| 46 | #else
|
|---|
| 47 |
|
|---|
| 48 | #define ULONG_MAX ((unsigned long)(~0L)) /* 0xFFFFFFFF for 32-bits */
|
|---|
| 49 | #define LONG_MAX ((long)(ULONG_MAX >> 1)) /* 0x7FFFFFFF for 32-bits*/
|
|---|
| 50 |
|
|---|
| 51 | #ifdef __STDC__
|
|---|
| 52 | # define PTR void *
|
|---|
| 53 | # ifndef NULL
|
|---|
| 54 | # define NULL (void *) 0
|
|---|
| 55 | # endif
|
|---|
| 56 | #else
|
|---|
| 57 | # define PTR char *
|
|---|
| 58 | # ifndef NULL
|
|---|
| 59 | # define NULL (void *) 0
|
|---|
| 60 | # endif
|
|---|
| 61 | #endif
|
|---|
| 62 |
|
|---|
| 63 | #endif
|
|---|
| 64 |
|
|---|
| 65 | long int random ();
|
|---|
| 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 initiallized
|
|---|
| 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 zeroeth 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 zeroeth 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 thi
|
|---|
| 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 | static int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 };
|
|---|
| 145 | static int seps[MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
|
|---|
| 146 |
|
|---|
| 147 |
|
|---|
| 148 |
|
|---|
| 149 | /* Initially, everything is set up as if from:
|
|---|
| 150 | initstate(1, randtbl, 128);
|
|---|
| 151 | Note that this initialization takes advantage of the fact that srandom
|
|---|
| 152 | advances the front and rear pointers 10*rand_deg times, and hence the
|
|---|
| 153 | rear pointer which starts at 0 will also end up at zero; thus the zeroeth
|
|---|
| 154 | element of the state information, which contains info about the current
|
|---|
| 155 | position of the rear pointer is just
|
|---|
| 156 | (MAX_TYPES * (rptr - state)) + TYPE_3 == TYPE_3. */
|
|---|
| 157 |
|
|---|
| 158 | static long int randtbl[DEG_3 + 1] =
|
|---|
| 159 | { TYPE_3,
|
|---|
| 160 | 0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342,
|
|---|
| 161 | 0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb,
|
|---|
| 162 | 0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
|
|---|
| 163 | 0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86,
|
|---|
| 164 | 0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7,
|
|---|
| 165 | 0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
|
|---|
| 166 | 0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b,
|
|---|
| 167 | 0xf5ad9d0e, 0x8999220b, 0x27fb47b9
|
|---|
| 168 | };
|
|---|
| 169 |
|
|---|
| 170 | /* FPTR and RPTR are two pointers into the state info, a front and a rear
|
|---|
| 171 | pointer. These two pointers are always rand_sep places aparts, as they
|
|---|
| 172 | cycle through the state information. (Yes, this does mean we could get
|
|---|
| 173 | away with just one pointer, but the code for random is more efficient
|
|---|
| 174 | this way). The pointers are left positioned as they would be from the call:
|
|---|
| 175 | initstate(1, randtbl, 128);
|
|---|
| 176 | (The position of the rear pointer, rptr, is really 0 (as explained above
|
|---|
| 177 | in the initialization of randtbl) because the state table pointer is set
|
|---|
| 178 | to point to randtbl[1] (as explained below).) */
|
|---|
| 179 |
|
|---|
| 180 | static long int *fptr = &randtbl[SEP_3 + 1];
|
|---|
| 181 | static long int *rptr = &randtbl[1];
|
|---|
| 182 |
|
|---|
| 183 |
|
|---|
| 184 |
|
|---|
| 185 | /* The following things are the pointer to the state information table,
|
|---|
| 186 | the type of the current generator, the degree of the current polynomial
|
|---|
| 187 | being used, and the separation between the two pointers.
|
|---|
| 188 | Note that for efficiency of random, we remember the first location of
|
|---|
| 189 | the state information, not the zeroeth. Hence it is valid to access
|
|---|
| 190 | state[-1], which is used to store the type of the R.N.G.
|
|---|
| 191 | Also, we remember the last location, since this is more efficient than
|
|---|
| 192 | indexing every time to find the address of the last element to see if
|
|---|
| 193 | the front and rear pointers have wrapped. */
|
|---|
| 194 |
|
|---|
| 195 | static long int *state = &randtbl[1];
|
|---|
| 196 |
|
|---|
| 197 | static int rand_type = TYPE_3;
|
|---|
| 198 | static int rand_deg = DEG_3;
|
|---|
| 199 | static int rand_sep = SEP_3;
|
|---|
| 200 |
|
|---|
| 201 | static long int *end_ptr = &randtbl[sizeof(randtbl) / sizeof(randtbl[0])];
|
|---|
| 202 | |
|---|
| 203 |
|
|---|
| 204 | /* Initialize the random number generator based on the given seed. If the
|
|---|
| 205 | type is the trivial no-state-information type, just remember the seed.
|
|---|
| 206 | Otherwise, initializes state[] based on the given "seed" via a linear
|
|---|
| 207 | congruential generator. Then, the pointers are set to known locations
|
|---|
| 208 | that are exactly rand_sep places apart. Lastly, it cycles the state
|
|---|
| 209 | information a given number of times to get rid of any initial dependencies
|
|---|
| 210 | introduced by the L.C.R.N.G. Note that the initialization of randtbl[]
|
|---|
| 211 | for default usage relies on values produced by this routine. */
|
|---|
| 212 | void
|
|---|
| 213 | srandom (x)
|
|---|
| 214 | unsigned int x;
|
|---|
| 215 | {
|
|---|
| 216 | state[0] = x;
|
|---|
| 217 | if (rand_type != TYPE_0)
|
|---|
| 218 | {
|
|---|
| 219 | register long int i;
|
|---|
| 220 | for (i = 1; i < rand_deg; ++i)
|
|---|
| 221 | state[i] = (1103515145 * state[i - 1]) + 12345;
|
|---|
| 222 | fptr = &state[rand_sep];
|
|---|
| 223 | rptr = &state[0];
|
|---|
| 224 | for (i = 0; i < 10 * rand_deg; ++i)
|
|---|
| 225 | random();
|
|---|
| 226 | }
|
|---|
| 227 | }
|
|---|
| 228 | |
|---|
| 229 |
|
|---|
| 230 | /* Initialize the state information in the given array of N bytes for
|
|---|
| 231 | future random number generation. Based on the number of bytes we
|
|---|
| 232 | are given, and the break values for the different R.N.G.'s, we choose
|
|---|
| 233 | the best (largest) one we can and set things up for it. srandom is
|
|---|
| 234 | then called to initialize the state information. Note that on return
|
|---|
| 235 | from srandom, we set state[-1] to be the type multiplexed with the current
|
|---|
| 236 | value of the rear pointer; this is so successive calls to initstate won't
|
|---|
| 237 | lose this information and will be able to restart with setstate.
|
|---|
| 238 | Note: The first thing we do is save the current state, if any, just like
|
|---|
| 239 | setstate so that it doesn't matter when initstate is called.
|
|---|
| 240 | Returns a pointer to the old state. */
|
|---|
| 241 | PTR
|
|---|
| 242 | initstate (seed, arg_state, n)
|
|---|
| 243 | unsigned int seed;
|
|---|
| 244 | PTR arg_state;
|
|---|
| 245 | unsigned long n;
|
|---|
| 246 | {
|
|---|
| 247 | PTR ostate = (PTR) &state[-1];
|
|---|
| 248 |
|
|---|
| 249 | if (rand_type == TYPE_0)
|
|---|
| 250 | state[-1] = rand_type;
|
|---|
| 251 | else
|
|---|
| 252 | state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
|
|---|
| 253 | if (n < BREAK_1)
|
|---|
| 254 | {
|
|---|
| 255 | if (n < BREAK_0)
|
|---|
| 256 | {
|
|---|
| 257 | errno = EINVAL;
|
|---|
| 258 | return NULL;
|
|---|
| 259 | }
|
|---|
| 260 | rand_type = TYPE_0;
|
|---|
| 261 | rand_deg = DEG_0;
|
|---|
| 262 | rand_sep = SEP_0;
|
|---|
| 263 | }
|
|---|
| 264 | else if (n < BREAK_2)
|
|---|
| 265 | {
|
|---|
| 266 | rand_type = TYPE_1;
|
|---|
| 267 | rand_deg = DEG_1;
|
|---|
| 268 | rand_sep = SEP_1;
|
|---|
| 269 | }
|
|---|
| 270 | else if (n < BREAK_3)
|
|---|
| 271 | {
|
|---|
| 272 | rand_type = TYPE_2;
|
|---|
| 273 | rand_deg = DEG_2;
|
|---|
| 274 | rand_sep = SEP_2;
|
|---|
| 275 | }
|
|---|
| 276 | else if (n < BREAK_4)
|
|---|
| 277 | {
|
|---|
| 278 | rand_type = TYPE_3;
|
|---|
| 279 | rand_deg = DEG_3;
|
|---|
| 280 | rand_sep = SEP_3;
|
|---|
| 281 | }
|
|---|
| 282 | else
|
|---|
| 283 | {
|
|---|
| 284 | rand_type = TYPE_4;
|
|---|
| 285 | rand_deg = DEG_4;
|
|---|
| 286 | rand_sep = SEP_4;
|
|---|
| 287 | }
|
|---|
| 288 |
|
|---|
| 289 | state = &((long int *) arg_state)[1]; /* First location. */
|
|---|
| 290 | /* Must set END_PTR before srandom. */
|
|---|
| 291 | end_ptr = &state[rand_deg];
|
|---|
| 292 | srandom(seed);
|
|---|
| 293 | if (rand_type == TYPE_0)
|
|---|
| 294 | state[-1] = rand_type;
|
|---|
| 295 | else
|
|---|
| 296 | state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
|
|---|
| 297 |
|
|---|
| 298 | return ostate;
|
|---|
| 299 | }
|
|---|
| 300 | |
|---|
| 301 |
|
|---|
| 302 | /* Restore the state from the given state array.
|
|---|
| 303 | Note: It is important that we also remember the locations of the pointers
|
|---|
| 304 | in the current state information, and restore the locations of the pointers
|
|---|
| 305 | from the old state information. This is done by multiplexing the pointer
|
|---|
| 306 | location into the zeroeth word of the state information. Note that due
|
|---|
| 307 | to the order in which things are done, it is OK to call setstate with the
|
|---|
| 308 | same state as the current state
|
|---|
| 309 | Returns a pointer to the old state information. */
|
|---|
| 310 |
|
|---|
| 311 | PTR
|
|---|
| 312 | setstate (arg_state)
|
|---|
| 313 | PTR arg_state;
|
|---|
| 314 | {
|
|---|
| 315 | register long int *new_state = (long int *) arg_state;
|
|---|
| 316 | register int type = new_state[0] % MAX_TYPES;
|
|---|
| 317 | register int rear = new_state[0] / MAX_TYPES;
|
|---|
| 318 | PTR ostate = (PTR) &state[-1];
|
|---|
| 319 |
|
|---|
| 320 | if (rand_type == TYPE_0)
|
|---|
| 321 | state[-1] = rand_type;
|
|---|
| 322 | else
|
|---|
| 323 | state[-1] = (MAX_TYPES * (rptr - state)) + rand_type;
|
|---|
| 324 |
|
|---|
| 325 | switch (type)
|
|---|
| 326 | {
|
|---|
| 327 | case TYPE_0:
|
|---|
| 328 | case TYPE_1:
|
|---|
| 329 | case TYPE_2:
|
|---|
| 330 | case TYPE_3:
|
|---|
| 331 | case TYPE_4:
|
|---|
| 332 | rand_type = type;
|
|---|
| 333 | rand_deg = degrees[type];
|
|---|
| 334 | rand_sep = seps[type];
|
|---|
| 335 | break;
|
|---|
| 336 | default:
|
|---|
| 337 | /* State info munged. */
|
|---|
| 338 | errno = EINVAL;
|
|---|
| 339 | return NULL;
|
|---|
| 340 | }
|
|---|
| 341 |
|
|---|
| 342 | state = &new_state[1];
|
|---|
| 343 | if (rand_type != TYPE_0)
|
|---|
| 344 | {
|
|---|
| 345 | rptr = &state[rear];
|
|---|
| 346 | fptr = &state[(rear + rand_sep) % rand_deg];
|
|---|
| 347 | }
|
|---|
| 348 | /* Set end_ptr too. */
|
|---|
| 349 | end_ptr = &state[rand_deg];
|
|---|
| 350 |
|
|---|
| 351 | return ostate;
|
|---|
| 352 | }
|
|---|
| 353 | |
|---|
| 354 |
|
|---|
| 355 | /* If we are using the trivial TYPE_0 R.N.G., just do the old linear
|
|---|
| 356 | congruential bit. Otherwise, we do our fancy trinomial stuff, which is the
|
|---|
| 357 | same in all ther other cases due to all the global variables that have been
|
|---|
| 358 | set up. The basic operation is to add the number at the rear pointer into
|
|---|
| 359 | the one at the front pointer. Then both pointers are advanced to the next
|
|---|
| 360 | location cyclically in the table. The value returned is the sum generated,
|
|---|
| 361 | reduced to 31 bits by throwing away the "least random" low bit.
|
|---|
| 362 | Note: The code takes advantage of the fact that both the front and
|
|---|
| 363 | rear pointers can't wrap on the same call by not testing the rear
|
|---|
| 364 | pointer if the front one has wrapped. Returns a 31-bit random number. */
|
|---|
| 365 |
|
|---|
| 366 | long int
|
|---|
| 367 | random ()
|
|---|
| 368 | {
|
|---|
| 369 | if (rand_type == TYPE_0)
|
|---|
| 370 | {
|
|---|
| 371 | state[0] = ((state[0] * 1103515245) + 12345) & LONG_MAX;
|
|---|
| 372 | return state[0];
|
|---|
| 373 | }
|
|---|
| 374 | else
|
|---|
| 375 | {
|
|---|
| 376 | long int i;
|
|---|
| 377 | *fptr += *rptr;
|
|---|
| 378 | /* Chucking least random bit. */
|
|---|
| 379 | i = (*fptr >> 1) & LONG_MAX;
|
|---|
| 380 | ++fptr;
|
|---|
| 381 | if (fptr >= end_ptr)
|
|---|
| 382 | {
|
|---|
| 383 | fptr = state;
|
|---|
| 384 | ++rptr;
|
|---|
| 385 | }
|
|---|
| 386 | else
|
|---|
| 387 | {
|
|---|
| 388 | ++rptr;
|
|---|
| 389 | if (rptr >= end_ptr)
|
|---|
| 390 | rptr = state;
|
|---|
| 391 | }
|
|---|
| 392 | return i;
|
|---|
| 393 | }
|
|---|
| 394 | }
|
|---|