diff options
Diffstat (limited to 'src/lib/libc/stdlib/random.c')
| -rw-r--r-- | src/lib/libc/stdlib/random.c | 395 |
1 files changed, 395 insertions, 0 deletions
diff --git a/src/lib/libc/stdlib/random.c b/src/lib/libc/stdlib/random.c new file mode 100644 index 0000000000..48e892042b --- /dev/null +++ b/src/lib/libc/stdlib/random.c | |||
| @@ -0,0 +1,395 @@ | |||
| 1 | /* $OpenBSD: random.c,v 1.15 2005/11/30 07:51:02 otto Exp $ */ | ||
| 2 | /* | ||
| 3 | * Copyright (c) 1983 Regents of the University of California. | ||
| 4 | * All rights reserved. | ||
| 5 | * | ||
| 6 | * Redistribution and use in source and binary forms, with or without | ||
| 7 | * modification, are permitted provided that the following conditions | ||
| 8 | * are met: | ||
| 9 | * 1. Redistributions of source code must retain the above copyright | ||
| 10 | * notice, this list of conditions and the following disclaimer. | ||
| 11 | * 2. Redistributions in binary form must reproduce the above copyright | ||
| 12 | * notice, this list of conditions and the following disclaimer in the | ||
| 13 | * documentation and/or other materials provided with the distribution. | ||
| 14 | * 3. 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 | #include <sys/param.h> | ||
| 32 | #include <sys/sysctl.h> | ||
| 33 | #include <sys/time.h> | ||
| 34 | #include <fcntl.h> | ||
| 35 | #include <stdio.h> | ||
| 36 | #include <stdlib.h> | ||
| 37 | #include <unistd.h> | ||
| 38 | |||
| 39 | /* | ||
| 40 | * random.c: | ||
| 41 | * | ||
| 42 | * An improved random number generation package. In addition to the standard | ||
| 43 | * rand()/srand() like interface, this package also has a special state info | ||
| 44 | * interface. The initstate() routine is called with a seed, an array of | ||
| 45 | * bytes, and a count of how many bytes are being passed in; this array is | ||
| 46 | * then initialized to contain information for random number generation with | ||
| 47 | * that much state information. Good sizes for the amount of state | ||
| 48 | * information are 32, 64, 128, and 256 bytes. The state can be switched by | ||
| 49 | * calling the setstate() routine with the same array as was initiallized | ||
| 50 | * with initstate(). By default, the package runs with 128 bytes of state | ||
| 51 | * information and generates far better random numbers than a linear | ||
| 52 | * congruential generator. If the amount of state information is less than | ||
| 53 | * 32 bytes, a simple linear congruential R.N.G. is used. | ||
| 54 | * | ||
| 55 | * Internally, the state information is treated as an array of int32_t; the | ||
| 56 | * zeroeth element of the array is the type of R.N.G. being used (small | ||
| 57 | * integer); the remainder of the array is the state information for the | ||
| 58 | * R.N.G. Thus, 32 bytes of state information will give 7 int32_ts worth of | ||
| 59 | * state information, which will allow a degree seven polynomial. (Note: | ||
| 60 | * the zeroeth word of state information also has some other information | ||
| 61 | * stored in it -- see setstate() for details). | ||
| 62 | * | ||
| 63 | * The random number generation technique is a linear feedback shift register | ||
| 64 | * approach, employing trinomials (since there are fewer terms to sum up that | ||
| 65 | * way). In this approach, the least significant bit of all the numbers in | ||
| 66 | * the state table will act as a linear feedback shift register, and will | ||
| 67 | * have period 2^deg - 1 (where deg is the degree of the polynomial being | ||
| 68 | * used, assuming that the polynomial is irreducible and primitive). The | ||
| 69 | * higher order bits will have longer periods, since their values are also | ||
| 70 | * influenced by pseudo-random carries out of the lower bits. The total | ||
| 71 | * period of the generator is approximately deg*(2**deg - 1); thus doubling | ||
| 72 | * the amount of state information has a vast influence on the period of the | ||
| 73 | * generator. Note: the deg*(2**deg - 1) is an approximation only good for | ||
| 74 | * large deg, when the period of the shift register is the dominant factor. | ||
| 75 | * With deg equal to seven, the period is actually much longer than the | ||
| 76 | * 7*(2**7 - 1) predicted by this formula. | ||
| 77 | */ | ||
| 78 | |||
| 79 | /* | ||
| 80 | * For each of the currently supported random number generators, we have a | ||
| 81 | * break value on the amount of state information (you need at least this | ||
| 82 | * many bytes of state info to support this random number generator), a degree | ||
| 83 | * for the polynomial (actually a trinomial) that the R.N.G. is based on, and | ||
| 84 | * the separation between the two lower order coefficients of the trinomial. | ||
| 85 | */ | ||
| 86 | #define TYPE_0 0 /* linear congruential */ | ||
| 87 | #define BREAK_0 8 | ||
| 88 | #define DEG_0 0 | ||
| 89 | #define SEP_0 0 | ||
| 90 | |||
| 91 | #define TYPE_1 1 /* x**7 + x**3 + 1 */ | ||
| 92 | #define BREAK_1 32 | ||
| 93 | #define DEG_1 7 | ||
| 94 | #define SEP_1 3 | ||
| 95 | |||
| 96 | #define TYPE_2 2 /* x**15 + x + 1 */ | ||
| 97 | #define BREAK_2 64 | ||
| 98 | #define DEG_2 15 | ||
| 99 | #define SEP_2 1 | ||
| 100 | |||
| 101 | #define TYPE_3 3 /* x**31 + x**3 + 1 */ | ||
| 102 | #define BREAK_3 128 | ||
| 103 | #define DEG_3 31 | ||
| 104 | #define SEP_3 3 | ||
| 105 | |||
| 106 | #define TYPE_4 4 /* x**63 + x + 1 */ | ||
| 107 | #define BREAK_4 256 | ||
| 108 | #define DEG_4 63 | ||
| 109 | #define SEP_4 1 | ||
| 110 | |||
| 111 | /* | ||
| 112 | * Array versions of the above information to make code run faster -- | ||
| 113 | * relies on fact that TYPE_i == i. | ||
| 114 | */ | ||
| 115 | #define MAX_TYPES 5 /* max number of types above */ | ||
| 116 | |||
| 117 | static int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 }; | ||
| 118 | static int seps [MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 }; | ||
| 119 | |||
| 120 | /* | ||
| 121 | * Initially, everything is set up as if from: | ||
| 122 | * | ||
| 123 | * initstate(1, &randtbl, 128); | ||
| 124 | * | ||
| 125 | * Note that this initialization takes advantage of the fact that srandom() | ||
| 126 | * advances the front and rear pointers 10*rand_deg times, and hence the | ||
| 127 | * rear pointer which starts at 0 will also end up at zero; thus the zeroeth | ||
| 128 | * element of the state information, which contains info about the current | ||
| 129 | * position of the rear pointer is just | ||
| 130 | * | ||
| 131 | * MAX_TYPES * (rptr - state) + TYPE_3 == TYPE_3. | ||
| 132 | */ | ||
| 133 | |||
| 134 | static int32_t randtbl[DEG_3 + 1] = { | ||
| 135 | TYPE_3, | ||
| 136 | 0x991539b1, 0x16a5bce3, 0x6774a4cd, 0x3e01511e, 0x4e508aaa, 0x61048c05, | ||
| 137 | 0xf5500617, 0x846b7115, 0x6a19892c, 0x896a97af, 0xdb48f936, 0x14898454, | ||
| 138 | 0x37ffd106, 0xb58bff9c, 0x59e17104, 0xcf918a49, 0x09378c83, 0x52c7a471, | ||
| 139 | 0x8d293ea9, 0x1f4fc301, 0xc3db71be, 0x39b44e1c, 0xf8a44ef9, 0x4c8b80b1, | ||
| 140 | 0x19edc328, 0x87bf4bdd, 0xc9b240e5, 0xe9ee4b1b, 0x4382aee7, 0x535b6b41, | ||
| 141 | 0xf3bec5da, | ||
| 142 | }; | ||
| 143 | |||
| 144 | /* | ||
| 145 | * fptr and rptr are two pointers into the state info, a front and a rear | ||
| 146 | * pointer. These two pointers are always rand_sep places aparts, as they | ||
| 147 | * cycle cyclically through the state information. (Yes, this does mean we | ||
| 148 | * could get away with just one pointer, but the code for random() is more | ||
| 149 | * efficient this way). The pointers are left positioned as they would be | ||
| 150 | * from the call | ||
| 151 | * | ||
| 152 | * initstate(1, randtbl, 128); | ||
| 153 | * | ||
| 154 | * (The position of the rear pointer, rptr, is really 0 (as explained above | ||
| 155 | * in the initialization of randtbl) because the state table pointer is set | ||
| 156 | * to point to randtbl[1] (as explained below). | ||
| 157 | */ | ||
| 158 | static int32_t *fptr = &randtbl[SEP_3 + 1]; | ||
| 159 | static int32_t *rptr = &randtbl[1]; | ||
| 160 | |||
| 161 | /* | ||
| 162 | * The following things are the pointer to the state information table, the | ||
| 163 | * type of the current generator, the degree of the current polynomial being | ||
| 164 | * used, and the separation between the two pointers. Note that for efficiency | ||
| 165 | * of random(), we remember the first location of the state information, not | ||
| 166 | * the zeroeth. Hence it is valid to access state[-1], which is used to | ||
| 167 | * store the type of the R.N.G. Also, we remember the last location, since | ||
| 168 | * this is more efficient than indexing every time to find the address of | ||
| 169 | * the last element to see if the front and rear pointers have wrapped. | ||
| 170 | */ | ||
| 171 | static int32_t *state = &randtbl[1]; | ||
| 172 | static int32_t *end_ptr = &randtbl[DEG_3 + 1]; | ||
| 173 | static int rand_type = TYPE_3; | ||
| 174 | static int rand_deg = DEG_3; | ||
| 175 | static int rand_sep = SEP_3; | ||
| 176 | |||
| 177 | /* | ||
| 178 | * srandom: | ||
| 179 | * | ||
| 180 | * Initialize the random number generator based on the given seed. If the | ||
| 181 | * type is the trivial no-state-information type, just remember the seed. | ||
| 182 | * Otherwise, initializes state[] based on the given "seed" via a linear | ||
| 183 | * congruential generator. Then, the pointers are set to known locations | ||
| 184 | * that are exactly rand_sep places apart. Lastly, it cycles the state | ||
| 185 | * information a given number of times to get rid of any initial dependencies | ||
| 186 | * introduced by the L.C.R.N.G. Note that the initialization of randtbl[] | ||
| 187 | * for default usage relies on values produced by this routine. | ||
| 188 | */ | ||
| 189 | void | ||
| 190 | srandom(unsigned int x) | ||
| 191 | { | ||
| 192 | int i; | ||
| 193 | int32_t test; | ||
| 194 | div_t val; | ||
| 195 | |||
| 196 | if (rand_type == TYPE_0) | ||
| 197 | state[0] = x; | ||
| 198 | else { | ||
| 199 | state[0] = x; | ||
| 200 | for (i = 1; i < rand_deg; i++) { | ||
| 201 | /* | ||
| 202 | * Implement the following, without overflowing 31 bits: | ||
| 203 | * | ||
| 204 | * state[i] = (16807 * state[i - 1]) % 2147483647; | ||
| 205 | * | ||
| 206 | * 2^31-1 (prime) = 2147483647 = 127773*16807+2836 | ||
| 207 | */ | ||
| 208 | val = div(state[i-1], 127773); | ||
| 209 | test = 16807 * val.rem - 2836 * val.quot; | ||
| 210 | state[i] = test + (test < 0 ? 2147483647 : 0); | ||
| 211 | } | ||
| 212 | fptr = &state[rand_sep]; | ||
| 213 | rptr = &state[0]; | ||
| 214 | for (i = 0; i < 10 * rand_deg; i++) | ||
| 215 | (void)random(); | ||
| 216 | } | ||
| 217 | } | ||
| 218 | |||
| 219 | /* | ||
| 220 | * srandomdev: | ||
| 221 | * | ||
| 222 | * Many programs choose the seed value in a totally predictable manner. | ||
| 223 | * This often causes problems. We seed the generator using random | ||
| 224 | * data from the kernel. | ||
| 225 | * Note that this particular seeding procedure can generate states | ||
| 226 | * which are impossible to reproduce by calling srandom() with any | ||
| 227 | * value, since the succeeding terms in the state buffer are no longer | ||
| 228 | * derived from the LC algorithm applied to a fixed seed. | ||
| 229 | */ | ||
| 230 | void | ||
| 231 | srandomdev(void) | ||
| 232 | { | ||
| 233 | int mib[2]; | ||
| 234 | size_t len; | ||
| 235 | |||
| 236 | if (rand_type == TYPE_0) | ||
| 237 | len = sizeof(state[0]); | ||
| 238 | else | ||
| 239 | len = rand_deg * sizeof(state[0]); | ||
| 240 | |||
| 241 | mib[0] = CTL_KERN; | ||
| 242 | mib[1] = KERN_ARND; | ||
| 243 | sysctl(mib, 2, state, &len, NULL, 0); | ||
| 244 | |||
| 245 | if (rand_type != TYPE_0) { | ||
| 246 | fptr = &state[rand_sep]; | ||
| 247 | rptr = &state[0]; | ||
| 248 | } | ||
| 249 | } | ||
| 250 | |||
| 251 | /* | ||
| 252 | * initstate: | ||
| 253 | * | ||
| 254 | * Initialize the state information in the given array of n bytes for future | ||
| 255 | * random number generation. Based on the number of bytes we are given, and | ||
| 256 | * the break values for the different R.N.G.'s, we choose the best (largest) | ||
| 257 | * one we can and set things up for it. srandom() is then called to | ||
| 258 | * initialize the state information. | ||
| 259 | * | ||
| 260 | * Note that on return from srandom(), we set state[-1] to be the type | ||
| 261 | * multiplexed with the current value of the rear pointer; this is so | ||
| 262 | * successive calls to initstate() won't lose this information and will be | ||
| 263 | * able to restart with setstate(). | ||
| 264 | * | ||
| 265 | * Note: the first thing we do is save the current state, if any, just like | ||
| 266 | * setstate() so that it doesn't matter when initstate is called. | ||
| 267 | * | ||
| 268 | * Returns a pointer to the old state. | ||
| 269 | */ | ||
| 270 | char * | ||
| 271 | initstate(u_int seed, char *arg_state, size_t n) | ||
| 272 | { | ||
| 273 | char *ostate = (char *)(&state[-1]); | ||
| 274 | |||
| 275 | if (rand_type == TYPE_0) | ||
| 276 | state[-1] = rand_type; | ||
| 277 | else | ||
| 278 | state[-1] = MAX_TYPES * (rptr - state) + rand_type; | ||
| 279 | if (n < BREAK_0) | ||
| 280 | return(NULL); | ||
| 281 | if (n < BREAK_1) { | ||
| 282 | rand_type = TYPE_0; | ||
| 283 | rand_deg = DEG_0; | ||
| 284 | rand_sep = SEP_0; | ||
| 285 | } else if (n < BREAK_2) { | ||
| 286 | rand_type = TYPE_1; | ||
| 287 | rand_deg = DEG_1; | ||
| 288 | rand_sep = SEP_1; | ||
| 289 | } else if (n < BREAK_3) { | ||
| 290 | rand_type = TYPE_2; | ||
| 291 | rand_deg = DEG_2; | ||
| 292 | rand_sep = SEP_2; | ||
| 293 | } else if (n < BREAK_4) { | ||
| 294 | rand_type = TYPE_3; | ||
| 295 | rand_deg = DEG_3; | ||
| 296 | rand_sep = SEP_3; | ||
| 297 | } else { | ||
| 298 | rand_type = TYPE_4; | ||
| 299 | rand_deg = DEG_4; | ||
| 300 | rand_sep = SEP_4; | ||
| 301 | } | ||
| 302 | state = &(((int32_t *)arg_state)[1]); /* first location */ | ||
| 303 | end_ptr = &state[rand_deg]; /* must set end_ptr before srandom */ | ||
| 304 | srandom(seed); | ||
| 305 | if (rand_type == TYPE_0) | ||
| 306 | state[-1] = rand_type; | ||
| 307 | else | ||
| 308 | state[-1] = MAX_TYPES*(rptr - state) + rand_type; | ||
| 309 | return(ostate); | ||
| 310 | } | ||
| 311 | |||
| 312 | /* | ||
| 313 | * setstate: | ||
| 314 | * | ||
| 315 | * Restore the state from the given state array. | ||
| 316 | * | ||
| 317 | * Note: it is important that we also remember the locations of the pointers | ||
| 318 | * in the current state information, and restore the locations of the pointers | ||
| 319 | * from the old state information. This is done by multiplexing the pointer | ||
| 320 | * location into the zeroeth word of the state information. | ||
| 321 | * | ||
| 322 | * Note that due to the order in which things are done, it is OK to call | ||
| 323 | * setstate() with the same state as the current state. | ||
| 324 | * | ||
| 325 | * Returns a pointer to the old state information. | ||
| 326 | */ | ||
| 327 | char * | ||
| 328 | setstate(const char *arg_state) | ||
| 329 | { | ||
| 330 | int32_t *new_state = (int32_t *)arg_state; | ||
| 331 | int32_t type = new_state[0] % MAX_TYPES; | ||
| 332 | int32_t rear = new_state[0] / MAX_TYPES; | ||
| 333 | char *ostate = (char *)(&state[-1]); | ||
| 334 | |||
| 335 | if (rand_type == TYPE_0) | ||
| 336 | state[-1] = rand_type; | ||
| 337 | else | ||
| 338 | state[-1] = MAX_TYPES * (rptr - state) + rand_type; | ||
| 339 | switch(type) { | ||
| 340 | case TYPE_0: | ||
| 341 | case TYPE_1: | ||
| 342 | case TYPE_2: | ||
| 343 | case TYPE_3: | ||
| 344 | case TYPE_4: | ||
| 345 | rand_type = type; | ||
| 346 | rand_deg = degrees[type]; | ||
| 347 | rand_sep = seps[type]; | ||
| 348 | break; | ||
| 349 | default: | ||
| 350 | return(NULL); | ||
| 351 | } | ||
| 352 | state = &new_state[1]; | ||
| 353 | if (rand_type != TYPE_0) { | ||
| 354 | rptr = &state[rear]; | ||
| 355 | fptr = &state[(rear + rand_sep) % rand_deg]; | ||
| 356 | } | ||
| 357 | end_ptr = &state[rand_deg]; /* set end_ptr too */ | ||
| 358 | return(ostate); | ||
| 359 | } | ||
| 360 | |||
| 361 | /* | ||
| 362 | * random: | ||
| 363 | * | ||
| 364 | * If we are using the trivial TYPE_0 R.N.G., just do the old linear | ||
| 365 | * congruential bit. Otherwise, we do our fancy trinomial stuff, which is | ||
| 366 | * the same in all the other cases due to all the global variables that have | ||
| 367 | * been set up. The basic operation is to add the number at the rear pointer | ||
| 368 | * into the one at the front pointer. Then both pointers are advanced to | ||
| 369 | * the next location cyclically in the table. The value returned is the sum | ||
| 370 | * generated, reduced to 31 bits by throwing away the "least random" low bit. | ||
| 371 | * | ||
| 372 | * Note: the code takes advantage of the fact that both the front and | ||
| 373 | * rear pointers can't wrap on the same call by not testing the rear | ||
| 374 | * pointer if the front one has wrapped. | ||
| 375 | * | ||
| 376 | * Returns a 31-bit random number. | ||
| 377 | */ | ||
| 378 | long | ||
| 379 | random(void) | ||
| 380 | { | ||
| 381 | int32_t i; | ||
| 382 | |||
| 383 | if (rand_type == TYPE_0) | ||
| 384 | i = state[0] = (state[0] * 1103515245 + 12345) & 0x7fffffff; | ||
| 385 | else { | ||
| 386 | *fptr += *rptr; | ||
| 387 | i = (*fptr >> 1) & 0x7fffffff; /* chucking least random bit */ | ||
| 388 | if (++fptr >= end_ptr) { | ||
| 389 | fptr = state; | ||
| 390 | ++rptr; | ||
| 391 | } else if (++rptr >= end_ptr) | ||
| 392 | rptr = state; | ||
| 393 | } | ||
| 394 | return((long)i); | ||
| 395 | } | ||
