diff options
Diffstat (limited to 'src/lib/libc/stdlib/random.c')
| -rw-r--r-- | src/lib/libc/stdlib/random.c | 370 |
1 files changed, 370 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..79344f30f1 --- /dev/null +++ b/src/lib/libc/stdlib/random.c | |||
| @@ -0,0 +1,370 @@ | |||
| 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. All advertising materials mentioning features or use of this software | ||
| 14 | * must display the following acknowledgement: | ||
| 15 | * This product includes software developed by the University of | ||
| 16 | * California, Berkeley and its contributors. | ||
| 17 | * 4. Neither the name of the University nor the names of its contributors | ||
| 18 | * may be used to endorse or promote products derived from this software | ||
| 19 | * without specific prior written permission. | ||
| 20 | * | ||
| 21 | * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND | ||
| 22 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| 23 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
| 24 | * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE | ||
| 25 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | ||
| 26 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | ||
| 27 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | ||
| 28 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | ||
| 29 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | ||
| 30 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | ||
| 31 | * SUCH DAMAGE. | ||
| 32 | */ | ||
| 33 | |||
| 34 | #if defined(LIBC_SCCS) && !defined(lint) | ||
| 35 | static char *rcsid = "$OpenBSD: random.c,v 1.6 1998/02/07 02:16:25 millert Exp $"; | ||
| 36 | #endif /* LIBC_SCCS and not lint */ | ||
| 37 | |||
| 38 | #include <stdio.h> | ||
| 39 | #include <stdlib.h> | ||
| 40 | |||
| 41 | /* | ||
| 42 | * random.c: | ||
| 43 | * | ||
| 44 | * An improved random number generation package. In addition to the standard | ||
| 45 | * rand()/srand() like interface, this package also has a special state info | ||
| 46 | * interface. The initstate() routine is called with a seed, an array of | ||
| 47 | * bytes, and a count of how many bytes are being passed in; this array is | ||
| 48 | * then initialized to contain information for random number generation with | ||
| 49 | * that much state information. Good sizes for the amount of state | ||
| 50 | * information are 32, 64, 128, and 256 bytes. The state can be switched by | ||
| 51 | * calling the setstate() routine with the same array as was initiallized | ||
| 52 | * with initstate(). By default, the package runs with 128 bytes of state | ||
| 53 | * information and generates far better random numbers than a linear | ||
| 54 | * congruential generator. If the amount of state information is less than | ||
| 55 | * 32 bytes, a simple linear congruential R.N.G. is used. | ||
| 56 | * | ||
| 57 | * Internally, the state information is treated as an array of longs; the | ||
| 58 | * zeroeth element of the array is the type of R.N.G. being used (small | ||
| 59 | * integer); the remainder of the array is the state information for the | ||
| 60 | * R.N.G. Thus, 32 bytes of state information will give 7 longs worth of | ||
| 61 | * state information, which will allow a degree seven polynomial. (Note: | ||
| 62 | * the zeroeth word of state information also has some other information | ||
| 63 | * stored in it -- see setstate() for details). | ||
| 64 | * | ||
| 65 | * The random number generation technique is a linear feedback shift register | ||
| 66 | * approach, employing trinomials (since there are fewer terms to sum up that | ||
| 67 | * way). In this approach, the least significant bit of all the numbers in | ||
| 68 | * the state table will act as a linear feedback shift register, and will | ||
| 69 | * have period 2^deg - 1 (where deg is the degree of the polynomial being | ||
| 70 | * used, assuming that the polynomial is irreducible and primitive). The | ||
| 71 | * higher order bits will have longer periods, since their values are also | ||
| 72 | * influenced by pseudo-random carries out of the lower bits. The total | ||
| 73 | * period of the generator is approximately deg*(2**deg - 1); thus doubling | ||
| 74 | * the amount of state information has a vast influence on the period of the | ||
| 75 | * generator. Note: the deg*(2**deg - 1) is an approximation only good for | ||
| 76 | * large deg, when the period of the shift register is the dominant factor. | ||
| 77 | * With deg equal to seven, the period is actually much longer than the | ||
| 78 | * 7*(2**7 - 1) predicted by this formula. | ||
| 79 | */ | ||
| 80 | |||
| 81 | /* | ||
| 82 | * For each of the currently supported random number generators, we have a | ||
| 83 | * break value on the amount of state information (you need at least this | ||
| 84 | * many bytes of state info to support this random number generator), a degree | ||
| 85 | * for the polynomial (actually a trinomial) that the R.N.G. is based on, and | ||
| 86 | * the separation between the two lower order coefficients of the trinomial. | ||
| 87 | */ | ||
| 88 | #define TYPE_0 0 /* linear congruential */ | ||
| 89 | #define BREAK_0 8 | ||
| 90 | #define DEG_0 0 | ||
| 91 | #define SEP_0 0 | ||
| 92 | |||
| 93 | #define TYPE_1 1 /* x**7 + x**3 + 1 */ | ||
| 94 | #define BREAK_1 32 | ||
| 95 | #define DEG_1 7 | ||
| 96 | #define SEP_1 3 | ||
| 97 | |||
| 98 | #define TYPE_2 2 /* x**15 + x + 1 */ | ||
| 99 | #define BREAK_2 64 | ||
| 100 | #define DEG_2 15 | ||
| 101 | #define SEP_2 1 | ||
| 102 | |||
| 103 | #define TYPE_3 3 /* x**31 + x**3 + 1 */ | ||
| 104 | #define BREAK_3 128 | ||
| 105 | #define DEG_3 31 | ||
| 106 | #define SEP_3 3 | ||
| 107 | |||
| 108 | #define TYPE_4 4 /* x**63 + x + 1 */ | ||
| 109 | #define BREAK_4 256 | ||
| 110 | #define DEG_4 63 | ||
| 111 | #define SEP_4 1 | ||
| 112 | |||
| 113 | /* | ||
| 114 | * Array versions of the above information to make code run faster -- | ||
| 115 | * relies on fact that TYPE_i == i. | ||
| 116 | */ | ||
| 117 | #define MAX_TYPES 5 /* max number of types above */ | ||
| 118 | |||
| 119 | static int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 }; | ||
| 120 | static int seps [MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 }; | ||
| 121 | |||
| 122 | /* | ||
| 123 | * Initially, everything is set up as if from: | ||
| 124 | * | ||
| 125 | * initstate(1, &randtbl, 128); | ||
| 126 | * | ||
| 127 | * Note that this initialization takes advantage of the fact that srandom() | ||
| 128 | * advances the front and rear pointers 10*rand_deg times, and hence the | ||
| 129 | * rear pointer which starts at 0 will also end up at zero; thus the zeroeth | ||
| 130 | * element of the state information, which contains info about the current | ||
| 131 | * position of the rear pointer is just | ||
| 132 | * | ||
| 133 | * MAX_TYPES * (rptr - state) + TYPE_3 == TYPE_3. | ||
| 134 | */ | ||
| 135 | |||
| 136 | static long randtbl[DEG_3 + 1] = { | ||
| 137 | TYPE_3, | ||
| 138 | 0x991539b1, 0x16a5bce3, 0x6774a4cd, 0x3e01511e, 0x4e508aaa, 0x61048c05, | ||
| 139 | 0xf5500617, 0x846b7115, 0x6a19892c, 0x896a97af, 0xdb48f936, 0x14898454, | ||
| 140 | 0x37ffd106, 0xb58bff9c, 0x59e17104, 0xcf918a49, 0x09378c83, 0x52c7a471, | ||
| 141 | 0x8d293ea9, 0x1f4fc301, 0xc3db71be, 0x39b44e1c, 0xf8a44ef9, 0x4c8b80b1, | ||
| 142 | 0x19edc328, 0x87bf4bdd, 0xc9b240e5, 0xe9ee4b1b, 0x4382aee7, 0x535b6b41, | ||
| 143 | 0xf3bec5da, | ||
| 144 | }; | ||
| 145 | |||
| 146 | /* | ||
| 147 | * fptr and rptr are two pointers into the state info, a front and a rear | ||
| 148 | * pointer. These two pointers are always rand_sep places aparts, as they | ||
| 149 | * cycle cyclically through the state information. (Yes, this does mean we | ||
| 150 | * could get away with just one pointer, but the code for random() is more | ||
| 151 | * efficient this way). The pointers are left positioned as they would be | ||
| 152 | * from the call | ||
| 153 | * | ||
| 154 | * initstate(1, randtbl, 128); | ||
| 155 | * | ||
| 156 | * (The position of the rear pointer, rptr, is really 0 (as explained above | ||
| 157 | * in the initialization of randtbl) because the state table pointer is set | ||
| 158 | * to point to randtbl[1] (as explained below). | ||
| 159 | */ | ||
| 160 | static long *fptr = &randtbl[SEP_3 + 1]; | ||
| 161 | static long *rptr = &randtbl[1]; | ||
| 162 | |||
| 163 | /* | ||
| 164 | * The following things are the pointer to the state information table, the | ||
| 165 | * type of the current generator, the degree of the current polynomial being | ||
| 166 | * used, and the separation between the two pointers. Note that for efficiency | ||
| 167 | * of random(), we remember the first location of the state information, not | ||
| 168 | * the zeroeth. Hence it is valid to access state[-1], which is used to | ||
| 169 | * store the type of the R.N.G. Also, we remember the last location, since | ||
| 170 | * this is more efficient than indexing every time to find the address of | ||
| 171 | * the last element to see if the front and rear pointers have wrapped. | ||
| 172 | */ | ||
| 173 | static long *state = &randtbl[1]; | ||
| 174 | static int rand_type = TYPE_3; | ||
| 175 | static int rand_deg = DEG_3; | ||
| 176 | static int rand_sep = SEP_3; | ||
| 177 | static long *end_ptr = &randtbl[DEG_3 + 1]; | ||
| 178 | |||
| 179 | /* | ||
| 180 | * srandom: | ||
| 181 | * | ||
| 182 | * Initialize the random number generator based on the given seed. If the | ||
| 183 | * type is the trivial no-state-information type, just remember the seed. | ||
| 184 | * Otherwise, initializes state[] based on the given "seed" via a linear | ||
| 185 | * congruential generator. Then, the pointers are set to known locations | ||
| 186 | * that are exactly rand_sep places apart. Lastly, it cycles the state | ||
| 187 | * information a given number of times to get rid of any initial dependencies | ||
| 188 | * introduced by the L.C.R.N.G. Note that the initialization of randtbl[] | ||
| 189 | * for default usage relies on values produced by this routine. | ||
| 190 | */ | ||
| 191 | void | ||
| 192 | srandom(x) | ||
| 193 | u_int x; | ||
| 194 | { | ||
| 195 | register long int test; | ||
| 196 | register int i; | ||
| 197 | ldiv_t val; | ||
| 198 | |||
| 199 | if (rand_type == TYPE_0) | ||
| 200 | state[0] = x; | ||
| 201 | else { | ||
| 202 | state[0] = x; | ||
| 203 | for (i = 1; i < rand_deg; i++) { | ||
| 204 | /* | ||
| 205 | * Implement the following, without overflowing 31 bits: | ||
| 206 | * | ||
| 207 | * state[i] = (16807 * state[i - 1]) % 2147483647; | ||
| 208 | * | ||
| 209 | * 2^31-1 (prime) = 2147483647 = 127773*16807+2836 | ||
| 210 | */ | ||
| 211 | val = ldiv(state[i-1], 127773); | ||
| 212 | test = 16807 * val.rem - 2836 * val.quot; | ||
| 213 | state[i] = test + (test < 0 ? 2147483647 : 0); | ||
| 214 | } | ||
| 215 | fptr = &state[rand_sep]; | ||
| 216 | rptr = &state[0]; | ||
| 217 | for (i = 0; i < 10 * rand_deg; i++) | ||
| 218 | (void)random(); | ||
| 219 | } | ||
| 220 | } | ||
| 221 | |||
| 222 | /* | ||
| 223 | * initstate: | ||
| 224 | * | ||
| 225 | * Initialize the state information in the given array of n bytes for future | ||
| 226 | * random number generation. Based on the number of bytes we are given, and | ||
| 227 | * the break values for the different R.N.G.'s, we choose the best (largest) | ||
| 228 | * one we can and set things up for it. srandom() is then called to | ||
| 229 | * initialize the state information. | ||
| 230 | * | ||
| 231 | * Note that on return from srandom(), we set state[-1] to be the type | ||
| 232 | * multiplexed with the current value of the rear pointer; this is so | ||
| 233 | * successive calls to initstate() won't lose this information and will be | ||
| 234 | * able to restart with setstate(). | ||
| 235 | * | ||
| 236 | * Note: the first thing we do is save the current state, if any, just like | ||
| 237 | * setstate() so that it doesn't matter when initstate is called. | ||
| 238 | * | ||
| 239 | * Returns a pointer to the old state. | ||
| 240 | */ | ||
| 241 | char * | ||
| 242 | initstate(seed, arg_state, n) | ||
| 243 | u_int seed; /* seed for R.N.G. */ | ||
| 244 | char *arg_state; /* pointer to state array */ | ||
| 245 | size_t n; /* # bytes of state info */ | ||
| 246 | { | ||
| 247 | register char *ostate = (char *)(&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_0) | ||
| 254 | return(NULL); | ||
| 255 | if (n < BREAK_1) { | ||
| 256 | rand_type = TYPE_0; | ||
| 257 | rand_deg = DEG_0; | ||
| 258 | rand_sep = SEP_0; | ||
| 259 | } else if (n < BREAK_2) { | ||
| 260 | rand_type = TYPE_1; | ||
| 261 | rand_deg = DEG_1; | ||
| 262 | rand_sep = SEP_1; | ||
| 263 | } else if (n < BREAK_3) { | ||
| 264 | rand_type = TYPE_2; | ||
| 265 | rand_deg = DEG_2; | ||
| 266 | rand_sep = SEP_2; | ||
| 267 | } else if (n < BREAK_4) { | ||
| 268 | rand_type = TYPE_3; | ||
| 269 | rand_deg = DEG_3; | ||
| 270 | rand_sep = SEP_3; | ||
| 271 | } else { | ||
| 272 | rand_type = TYPE_4; | ||
| 273 | rand_deg = DEG_4; | ||
| 274 | rand_sep = SEP_4; | ||
| 275 | } | ||
| 276 | state = &(((long *)arg_state)[1]); /* first location */ | ||
| 277 | end_ptr = &state[rand_deg]; /* must set end_ptr before srandom */ | ||
| 278 | srandom(seed); | ||
| 279 | if (rand_type == TYPE_0) | ||
| 280 | state[-1] = rand_type; | ||
| 281 | else | ||
| 282 | state[-1] = MAX_TYPES*(rptr - state) + rand_type; | ||
| 283 | return(ostate); | ||
| 284 | } | ||
| 285 | |||
| 286 | /* | ||
| 287 | * setstate: | ||
| 288 | * | ||
| 289 | * Restore the state from the given state array. | ||
| 290 | * | ||
| 291 | * Note: it is important that we also remember the locations of the pointers | ||
| 292 | * in the current state information, and restore the locations of the pointers | ||
| 293 | * from the old state information. This is done by multiplexing the pointer | ||
| 294 | * location into the zeroeth word of the state information. | ||
| 295 | * | ||
| 296 | * Note that due to the order in which things are done, it is OK to call | ||
| 297 | * setstate() with the same state as the current state. | ||
| 298 | * | ||
| 299 | * Returns a pointer to the old state information. | ||
| 300 | */ | ||
| 301 | char * | ||
| 302 | setstate(arg_state) | ||
| 303 | const char *arg_state; | ||
| 304 | { | ||
| 305 | register long *new_state = (long *)arg_state; | ||
| 306 | register int type = new_state[0] % MAX_TYPES; | ||
| 307 | register int rear = new_state[0] / MAX_TYPES; | ||
| 308 | char *ostate = (char *)(&state[-1]); | ||
| 309 | |||
| 310 | if (rand_type == TYPE_0) | ||
| 311 | state[-1] = rand_type; | ||
| 312 | else | ||
| 313 | state[-1] = MAX_TYPES * (rptr - state) + rand_type; | ||
| 314 | switch(type) { | ||
| 315 | case TYPE_0: | ||
| 316 | case TYPE_1: | ||
| 317 | case TYPE_2: | ||
| 318 | case TYPE_3: | ||
| 319 | case TYPE_4: | ||
| 320 | rand_type = type; | ||
| 321 | rand_deg = degrees[type]; | ||
| 322 | rand_sep = seps[type]; | ||
| 323 | break; | ||
| 324 | default: | ||
| 325 | return(NULL); | ||
| 326 | } | ||
| 327 | state = &new_state[1]; | ||
| 328 | if (rand_type != TYPE_0) { | ||
| 329 | rptr = &state[rear]; | ||
| 330 | fptr = &state[(rear + rand_sep) % rand_deg]; | ||
| 331 | } | ||
| 332 | end_ptr = &state[rand_deg]; /* set end_ptr too */ | ||
| 333 | return(ostate); | ||
| 334 | } | ||
| 335 | |||
| 336 | /* | ||
| 337 | * random: | ||
| 338 | * | ||
| 339 | * If we are using the trivial TYPE_0 R.N.G., just do the old linear | ||
| 340 | * congruential bit. Otherwise, we do our fancy trinomial stuff, which is | ||
| 341 | * the same in all the other cases due to all the global variables that have | ||
| 342 | * been set up. The basic operation is to add the number at the rear pointer | ||
| 343 | * into the one at the front pointer. Then both pointers are advanced to | ||
| 344 | * the next location cyclically in the table. The value returned is the sum | ||
| 345 | * generated, reduced to 31 bits by throwing away the "least random" low bit. | ||
| 346 | * | ||
| 347 | * Note: the code takes advantage of the fact that both the front and | ||
| 348 | * rear pointers can't wrap on the same call by not testing the rear | ||
| 349 | * pointer if the front one has wrapped. | ||
| 350 | * | ||
| 351 | * Returns a 31-bit random number. | ||
| 352 | */ | ||
| 353 | long | ||
| 354 | random() | ||
| 355 | { | ||
| 356 | long i; | ||
| 357 | |||
| 358 | if (rand_type == TYPE_0) | ||
| 359 | i = state[0] = (state[0] * 1103515245 + 12345) & 0x7fffffff; | ||
| 360 | else { | ||
| 361 | *fptr += *rptr; | ||
| 362 | i = (*fptr >> 1) & 0x7fffffff; /* chucking least random bit */ | ||
| 363 | if (++fptr >= end_ptr) { | ||
| 364 | fptr = state; | ||
| 365 | ++rptr; | ||
| 366 | } else if (++rptr >= end_ptr) | ||
| 367 | rptr = state; | ||
| 368 | } | ||
| 369 | return(i); | ||
| 370 | } | ||
