random.c 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453
  1. /*
  2. * Copyright (c) 1983, 1993
  3. * The Regents of the University of California. 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. * From:
  35. * NetBSD: random.c,v 1.19 2000/01/22 22:19:20 mycroft Exp
  36. *
  37. * Hacked gruesomely for OS/161.
  38. */
  39. #include <assert.h>
  40. #include <errno.h>
  41. #include <stdlib.h>
  42. /*
  43. * For a thread-safe libc, declare a lock for this file and change
  44. * these to be nonempty.
  45. */
  46. #define LOCKME()
  47. #define UNLOCKME()
  48. static void srandom_unlocked(unsigned long);
  49. static long random_unlocked(void);
  50. /*
  51. * random.c:
  52. *
  53. * An improved random number generation package. In addition to the standard
  54. * rand()/srand() like interface, this package also has a special state info
  55. * interface. The initstate() routine is called with a seed, an array of
  56. * bytes, and a count of how many bytes are being passed in; this array is
  57. * then initialized to contain information for random number generation with
  58. * that much state information. Good sizes for the amount of state
  59. * information are 32, 64, 128, and 256 bytes. The state can be switched by
  60. * calling the setstate() routine with the same array as was initiallized
  61. * with initstate(). By default, the package runs with 128 bytes of state
  62. * information and generates far better random numbers than a linear
  63. * congruential generator. If the amount of state information is less than
  64. * 32 bytes, a simple linear congruential R.N.G. is used.
  65. *
  66. * Internally, the state information is treated as an array of longs; the
  67. * zeroeth element of the array is the type of R.N.G. being used (small
  68. * integer); the remainder of the array is the state information for the
  69. * R.N.G. Thus, 32 bytes of state information will give 7 longs worth of
  70. * state information, which will allow a degree seven polynomial. (Note:
  71. * the zeroeth word of state information also has some other information
  72. * stored in it -- see setstate() for details).
  73. *
  74. * The random number generation technique is a linear feedback shift register
  75. * approach, employing trinomials (since there are fewer terms to sum up that
  76. * way). In this approach, the least significant bit of all the numbers in
  77. * the state table will act as a linear feedback shift register, and will
  78. * have period 2^deg - 1 (where deg is the degree of the polynomial being
  79. * used, assuming that the polynomial is irreducible and primitive). The
  80. * higher order bits will have longer periods, since their values are also
  81. * influenced by pseudo-random carries out of the lower bits. The total
  82. * period of the generator is approximately deg*(2**deg - 1); thus doubling
  83. * the amount of state information has a vast influence on the period of the
  84. * generator. Note: the deg*(2**deg - 1) is an approximation only good for
  85. * large deg, when the period of the shift register is the dominant factor.
  86. * With deg equal to seven, the period is actually much longer than the
  87. * 7*(2**7 - 1) predicted by this formula.
  88. *
  89. * Modified 28 December 1994 by Jacob S. Rosenberg.
  90. * The following changes have been made:
  91. * All references to the type u_int have been changed to unsigned long.
  92. * All references to type int have been changed to type long. Other
  93. * cleanups have been made as well. A warning for both initstate and
  94. * setstate has been inserted to the effect that on Sparc platforms
  95. * the 'arg_state' variable must be forced to begin on word boundaries.
  96. * This can be easily done by casting a long integer array to char *.
  97. * The overall logic has been left STRICTLY alone. This software was
  98. * tested on both a VAX and Sun SpacsStation with exactly the same
  99. * results. The new version and the original give IDENTICAL results.
  100. * The new version is somewhat faster than the original. As the
  101. * documentation says: "By default, the package runs with 128 bytes of
  102. * state information and generates far better random numbers than a linear
  103. * congruential generator. If the amount of state information is less than
  104. * 32 bytes, a simple linear congruential R.N.G. is used." For a buffer of
  105. * 128 bytes, this new version runs about 19 percent faster and for a 16
  106. * byte buffer it is about 5 percent faster.
  107. */
  108. /*
  109. * For each of the currently supported random number generators, we have a
  110. * break value on the amount of state information (you need at least this
  111. * many bytes of state info to support this random number generator), a degree
  112. * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
  113. * the separation between the two lower order coefficients of the trinomial.
  114. */
  115. #define TYPE_0 0 /* linear congruential */
  116. #define BREAK_0 8
  117. #define DEG_0 0
  118. #define SEP_0 0
  119. #define TYPE_1 1 /* x**7 + x**3 + 1 */
  120. #define BREAK_1 32
  121. #define DEG_1 7
  122. #define SEP_1 3
  123. #define TYPE_2 2 /* x**15 + x + 1 */
  124. #define BREAK_2 64
  125. #define DEG_2 15
  126. #define SEP_2 1
  127. #define TYPE_3 3 /* x**31 + x**3 + 1 */
  128. #define BREAK_3 128
  129. #define DEG_3 31
  130. #define SEP_3 3
  131. #define TYPE_4 4 /* x**63 + x + 1 */
  132. #define BREAK_4 256
  133. #define DEG_4 63
  134. #define SEP_4 1
  135. /*
  136. * Array versions of the above information to make code run faster --
  137. * relies on fact that TYPE_i == i.
  138. */
  139. #define MAX_TYPES 5 /* max number of types above */
  140. static const int degrees[MAX_TYPES] = { DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 };
  141. static const int seps[MAX_TYPES] = { SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
  142. /*
  143. * Initially, everything is set up as if from:
  144. *
  145. * initstate(1, &randtbl, 128);
  146. *
  147. * Note that this initialization takes advantage of the fact that srandom()
  148. * advances the front and rear pointers 10*rand_deg times, and hence the
  149. * rear pointer which starts at 0 will also end up at zero; thus the zeroeth
  150. * element of the state information, which contains info about the current
  151. * position of the rear pointer is just
  152. *
  153. * MAX_TYPES * (rptr - state) + TYPE_3 == TYPE_3.
  154. */
  155. static long randtbl[DEG_3 + 1] = {
  156. TYPE_3,
  157. (long)0x9a319039L, (long)0x32d9c024L, (long)0x9b663182L,
  158. (long)0x5da1f342L, (long)0xde3b81e0L, (long)0xdf0a6fb5L,
  159. (long)0xf103bc02L, (long)0x48f340fbL, (long)0x7449e56bL,
  160. (long)0xbeb1dbb0L, (long)0xab5c5918L, (long)0x946554fdL,
  161. (long)0x8c2e680fL, (long)0xeb3d799fL, (long)0xb11ee0b7L,
  162. (long)0x2d436b86L, (long)0xda672e2aL, (long)0x1588ca88L,
  163. (long)0xe369735dL, (long)0x904f35f7L, (long)0xd7158fd6L,
  164. (long)0x6fa6f051L, (long)0x616e6b96L, (long)0xac94efdcL,
  165. (long)0x36413f93L, (long)0xc622c298L, (long)0xf5a42ab8L,
  166. (long)0x8a88d77bL, (long)0xf5ad9d0eL, (long)0x8999220bL,
  167. (long)0x27fb47b9L,
  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 cyclically through the state information. (Yes, this does mean we
  173. * could get away with just one pointer, but the code for random() is more
  174. * efficient this way). The pointers are left positioned as they would be
  175. * from the call
  176. *
  177. * initstate(1, randtbl, 128);
  178. *
  179. * (The position of the rear pointer, rptr, is really 0 (as explained above
  180. * in the initialization of randtbl) because the state table pointer is set
  181. * to point to randtbl[1] (as explained below).
  182. */
  183. static long *fptr = &randtbl[SEP_3 + 1];
  184. static long *rptr = &randtbl[1];
  185. /*
  186. * The following things are the pointer to the state information table, the
  187. * type of the current generator, the degree of the current polynomial being
  188. * used, and the separation between the two pointers. Note that for efficiency
  189. * of random(), we remember the first location of the state information, not
  190. * the zeroeth. Hence it is valid to access state[-1], which is used to
  191. * store the type of the R.N.G. Also, we remember the last location, since
  192. * this is more efficient than indexing every time to find the address of
  193. * the last element to see if the front and rear pointers have wrapped.
  194. */
  195. static long *state = &randtbl[1];
  196. static long rand_type = TYPE_3;
  197. static int rand_deg = DEG_3;
  198. static int rand_sep = SEP_3;
  199. static long *end_ptr = &randtbl[DEG_3 + 1];
  200. /*
  201. * srandom:
  202. *
  203. * Initialize the random number generator based on the given seed. If the
  204. * type is the trivial no-state-information type, just remember the seed.
  205. * Otherwise, initializes state[] based on the given "seed" via a linear
  206. * congruential generator. Then, the pointers are set to known locations
  207. * that are exactly rand_sep places apart. Lastly, it cycles the state
  208. * information a given number of times to get rid of any initial dependencies
  209. * introduced by the L.C.R.N.G. Note that the initialization of randtbl[]
  210. * for default usage relies on values produced by this routine.
  211. */
  212. static
  213. void
  214. srandom_unlocked(unsigned long x)
  215. {
  216. int i;
  217. if (rand_type == TYPE_0)
  218. state[0] = x;
  219. else {
  220. state[0] = x;
  221. for (i = 1; i < rand_deg; i++)
  222. state[i] = 1103515245L * state[i - 1] + 12345L;
  223. fptr = &state[rand_sep];
  224. rptr = &state[0];
  225. for (i = 0; i < 10 * rand_deg; i++)
  226. (void)random_unlocked();
  227. }
  228. }
  229. void
  230. srandom(unsigned long x)
  231. {
  232. LOCKME();
  233. srandom_unlocked(x);
  234. UNLOCKME();
  235. }
  236. /*
  237. * initstate:
  238. *
  239. * Initialize the state information in the given array of n bytes for future
  240. * random number generation. Based on the number of bytes we are given, and
  241. * the break values for the different R.N.G.'s, we choose the best (largest)
  242. * one we can and set things up for it. srandom() is then called to
  243. * initialize the state information.
  244. *
  245. * Note that on return from srandom(), we set state[-1] to be the type
  246. * multiplexed with the current value of the rear pointer; this is so
  247. * successive calls to initstate() won't lose this information and will be
  248. * able to restart with setstate().
  249. *
  250. * Note: the first thing we do is save the current state, if any, just like
  251. * setstate() so that it doesn't matter when initstate is called.
  252. *
  253. * Returns a pointer to the old state.
  254. *
  255. * Note: The Sparc platform requires that arg_state begin on a long
  256. * word boundary; otherwise a bus error will occur. Even so, lint will
  257. * complain about mis-alignment, but you should disregard these messages.
  258. */
  259. char *
  260. initstate(
  261. unsigned long seed, /* seed for R.N.G. */
  262. char *arg_state, /* pointer to state array */
  263. size_t n) /* # bytes of state info */
  264. {
  265. void *ostate = (void *)(&state[-1]);
  266. long *long_arg_state;
  267. assert(arg_state != NULL);
  268. long_arg_state = (long *)(void *)arg_state;
  269. LOCKME();
  270. if (rand_type == TYPE_0)
  271. state[-1] = rand_type;
  272. else
  273. state[-1] = MAX_TYPES * (rptr - state) + rand_type;
  274. if (n < BREAK_0) {
  275. UNLOCKME();
  276. return (NULL);
  277. } else if (n < BREAK_1) {
  278. rand_type = TYPE_0;
  279. rand_deg = DEG_0;
  280. rand_sep = SEP_0;
  281. } else if (n < BREAK_2) {
  282. rand_type = TYPE_1;
  283. rand_deg = DEG_1;
  284. rand_sep = SEP_1;
  285. } else if (n < BREAK_3) {
  286. rand_type = TYPE_2;
  287. rand_deg = DEG_2;
  288. rand_sep = SEP_2;
  289. } else if (n < BREAK_4) {
  290. rand_type = TYPE_3;
  291. rand_deg = DEG_3;
  292. rand_sep = SEP_3;
  293. } else {
  294. rand_type = TYPE_4;
  295. rand_deg = DEG_4;
  296. rand_sep = SEP_4;
  297. }
  298. state = (long *) (long_arg_state + 1); /* first location */
  299. end_ptr = &state[rand_deg]; /* must set end_ptr before srandom */
  300. srandom_unlocked(seed);
  301. if (rand_type == TYPE_0)
  302. long_arg_state[0] = rand_type;
  303. else
  304. long_arg_state[0] = MAX_TYPES * (rptr - state) + rand_type;
  305. UNLOCKME();
  306. return((char *)ostate);
  307. }
  308. /*
  309. * setstate:
  310. *
  311. * Restore the state from the given state array.
  312. *
  313. * Note: it is important that we also remember the locations of the pointers
  314. * in the current state information, and restore the locations of the pointers
  315. * from the old state information. This is done by multiplexing the pointer
  316. * location into the zeroeth word of the state information.
  317. *
  318. * Note that due to the order in which things are done, it is OK to call
  319. * setstate() with the same state as the current state.
  320. *
  321. * Returns a pointer to the old state information.
  322. *
  323. * Note: The Sparc platform requires that arg_state begin on a long
  324. * word boundary; otherwise a bus error will occur. Even so, lint will
  325. * complain about mis-alignment, but you should disregard these messages.
  326. */
  327. char *
  328. setstate(char *arg_state) /* pointer to state array */
  329. {
  330. long *new_state;
  331. int type;
  332. int rear;
  333. void *ostate = (void *)(&state[-1]);
  334. assert(arg_state != NULL);
  335. new_state = (long *)(void *)arg_state;
  336. type = (int)(new_state[0] % MAX_TYPES);
  337. rear = (int)(new_state[0] / MAX_TYPES);
  338. LOCKME();
  339. if (rand_type == TYPE_0)
  340. state[-1] = rand_type;
  341. else
  342. state[-1] = MAX_TYPES * (rptr - state) + rand_type;
  343. switch(type) {
  344. case TYPE_0:
  345. case TYPE_1:
  346. case TYPE_2:
  347. case TYPE_3:
  348. case TYPE_4:
  349. rand_type = type;
  350. rand_deg = degrees[type];
  351. rand_sep = seps[type];
  352. break;
  353. default:
  354. UNLOCKME();
  355. return (NULL);
  356. }
  357. state = (long *) (new_state + 1);
  358. if (rand_type != TYPE_0) {
  359. rptr = &state[rear];
  360. fptr = &state[(rear + rand_sep) % rand_deg];
  361. }
  362. end_ptr = &state[rand_deg]; /* set end_ptr too */
  363. UNLOCKME();
  364. return((char *)ostate);
  365. }
  366. /*
  367. * random:
  368. *
  369. * If we are using the trivial TYPE_0 R.N.G., just do the old linear
  370. * congruential bit. Otherwise, we do our fancy trinomial stuff, which is
  371. * the same in all the other cases due to all the global variables that have
  372. * been set up. The basic operation is to add the number at the rear pointer
  373. * into the one at the front pointer. Then both pointers are advanced to
  374. * the next location cyclically in the table. The value returned is the sum
  375. * generated, reduced to 31 bits by throwing away the "least random" low bit.
  376. *
  377. * Note: the code takes advantage of the fact that both the front and
  378. * rear pointers can't wrap on the same call by not testing the rear
  379. * pointer if the front one has wrapped.
  380. *
  381. * Returns a 31-bit random number.
  382. */
  383. static
  384. long
  385. random_unlocked(void)
  386. {
  387. long i;
  388. long *f, *r;
  389. if (rand_type == TYPE_0) {
  390. i = state[0];
  391. state[0] = i = (i * 1103515245L + 12345L) & 0x7fffffff;
  392. } else {
  393. /*
  394. * Use local variables rather than static variables for speed.
  395. */
  396. f = fptr; r = rptr;
  397. *f += *r;
  398. /* chucking least random bit */
  399. i = ((unsigned long)*f >> 1) & 0x7fffffff;
  400. if (++f >= end_ptr) {
  401. f = state;
  402. ++r;
  403. }
  404. else if (++r >= end_ptr) {
  405. r = state;
  406. }
  407. fptr = f; rptr = r;
  408. }
  409. return(i);
  410. }
  411. long
  412. random(void)
  413. {
  414. long r;
  415. LOCKME();
  416. r = random_unlocked();
  417. UNLOCKME();
  418. return (r);
  419. }