git.zerfleddert.de Git - micropolis/blob - src/sim/random.c

   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 #include "sim.h"
  35
  36 #include <stdio.h>
  37 #include <stdlib.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 longs; 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 longs 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 QUAD sim_random();
 121 void sim_srandom();
 122 char *sim_initstate();
 123 char *sim_setstate();
 124
 125 /*
 126  * Initially, everything is set up as if from:
 127  *
 128  *      initstate(1, &randtbl, 128);
 129  *
 130  * Note that this initialization takes advantage of the fact that srandom()
 131  * advances the front and rear pointers 10*rand_deg times, and hence the
 132  * rear pointer which starts at 0 will also end up at zero; thus the zeroeth
 133  * element of the state information, which contains info about the current
 134  * position of the rear pointer is just
 135  *
 136  *      MAX_TYPES * (rptr - state) + TYPE_3 == TYPE_3.
 137  */
 138
 139 static QUAD randtbl[DEG_3 + 1] = {
 140         TYPE_3,
 141         0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342, 0xde3b81e0, 0xdf0a6fb5,
 142         0xf103bc02, 0x48f340fb, 0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
 143         0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86, 0xda672e2a, 0x1588ca88,
 144         0xe369735d, 0x904f35f7, 0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
 145         0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b, 0xf5ad9d0e, 0x8999220b,
 146         0x27fb47b9,
 147 };
 148
 149 /*
 150  * fptr and rptr are two pointers into the state info, a front and a rear
 151  * pointer.  These two pointers are always rand_sep places aparts, as they
 152  * cycle cyclically through the state information.  (Yes, this does mean we
 153  * could get away with just one pointer, but the code for random() is more
 154  * efficient this way).  The pointers are left positioned as they would be
 155  * from the call
 156  *
 157  *      initstate(1, randtbl, 128);
 158  *
 159  * (The position of the rear pointer, rptr, is really 0 (as explained above
 160  * in the initialization of randtbl) because the state table pointer is set
 161  * to point to randtbl[1] (as explained below).
 162  */
 163 static QUAD *fptr = &randtbl[SEP_3 + 1];
 164 static QUAD *rptr = &randtbl[1];
 165
 166 /*
 167  * The following things are the pointer to the state information table, the
 168  * type of the current generator, the degree of the current polynomial being
 169  * used, and the separation between the two pointers.  Note that for efficiency
 170  * of random(), we remember the first location of the state information, not
 171  * the zeroeth.  Hence it is valid to access state[-1], which is used to
 172  * store the type of the R.N.G.  Also, we remember the last location, since
 173  * this is more efficient than indexing every time to find the address of
 174  * the last element to see if the front and rear pointers have wrapped.
 175  */
 176 static QUAD *state = &randtbl[1];
 177 static int rand_type = TYPE_3;
 178 static int rand_deg = DEG_3;
 179 static int rand_sep = SEP_3;
 180 static QUAD *end_ptr = &randtbl[DEG_3 + 1];
 181
 182 /*
 183  * srandom:
 184  *
 185  * Initialize the random number generator based on the given seed.  If the
 186  * type is the trivial no-state-information type, just remember the seed.
 187  * Otherwise, initializes state[] based on the given "seed" via a linear
 188  * congruential generator.  Then, the pointers are set to known locations
 189  * that are exactly rand_sep places apart.  Lastly, it cycles the state
 190  * information a given number of times to get rid of any initial dependencies
 191  * introduced by the L.C.R.N.G.  Note that the initialization of randtbl[]
 192  * for default usage relies on values produced by this routine.
 193  */
 194 void
 195 sim_srandom(x)
 196         unsigned int x;
 197 {
 198         register int i, j;
 199
 200         if (rand_type == TYPE_0)
 201                 state[0] = x;
 202         else {
 203                 j = 1;
 204                 state[0] = x;
 205                 for (i = 1; i < rand_deg; i++)
 206                         state[i] = 1103515245 * state[i - 1] + 12345;
 207                 fptr = &state[rand_sep];
 208                 rptr = &state[0];
 209                 for (i = 0; i < 10 * rand_deg; i++)
 210                         (void)sim_random();
 211         }
 212 }
 213
 214 /*
 215  * initstate:
 216  *
 217  * Initialize the state information in the given array of n bytes for future
 218  * random number generation.  Based on the number of bytes we are given, and
 219  * the break values for the different R.N.G.'s, we choose the best (largest)
 220  * one we can and set things up for it.  srandom() is then called to
 221  * initialize the state information.
 222  *
 223  * Note that on return from srandom(), we set state[-1] to be the type
 224  * multiplexed with the current value of the rear pointer; this is so
 225  * successive calls to initstate() won't lose this information and will be
 226  * able to restart with setstate().
 227  *
 228  * Note: the first thing we do is save the current state, if any, just like
 229  * setstate() so that it doesn't matter when initstate is called.
 230  *
 231  * Returns a pointer to the old state.
 232  */
 233 char *
 234 sim_initstate(seed, arg_state, n)
 235         unsigned int seed;              /* seed for R.N.G. */
 236         char *arg_state;                /* pointer to state array */
 237         int n;                          /* # bytes of state info */
 238 {
 239         register char *ostate = (char *)(&state[-1]);
 240
 241         if (rand_type == TYPE_0)
 242                 state[-1] = rand_type;
 243         else
 244                 state[-1] = MAX_TYPES * (rptr - state) + rand_type;
 245         if (n < BREAK_0) {
 246                 (void)fprintf(stderr,
 247                     "random: not enough state (%d bytes); ignored.\n", n);
 248                 return(0);
 249         }
 250         if (n < BREAK_1) {
 251                 rand_type = TYPE_0;
 252                 rand_deg = DEG_0;
 253                 rand_sep = SEP_0;
 254         } else if (n < BREAK_2) {
 255                 rand_type = TYPE_1;
 256                 rand_deg = DEG_1;
 257                 rand_sep = SEP_1;
 258         } else if (n < BREAK_3) {
 259                 rand_type = TYPE_2;
 260                 rand_deg = DEG_2;
 261                 rand_sep = SEP_2;
 262         } else if (n < BREAK_4) {
 263                 rand_type = TYPE_3;
 264                 rand_deg = DEG_3;
 265                 rand_sep = SEP_3;
 266         } else {
 267                 rand_type = TYPE_4;
 268                 rand_deg = DEG_4;
 269                 rand_sep = SEP_4;
 270         }
 271         state = &(((QUAD *)arg_state)[1]);      /* first location */
 272         end_ptr = &state[rand_deg];     /* must set end_ptr before srandom */
 273         sim_srandom(seed);
 274         if (rand_type == TYPE_0)
 275                 state[-1] = rand_type;
 276         else
 277                 state[-1] = MAX_TYPES*(rptr - state) + rand_type;
 278         return(ostate);
 279 }
 280
 281 /*
 282  * setstate:
 283  *
 284  * Restore the state from the given state array.
 285  *
 286  * Note: it is important that we also remember the locations of the pointers
 287  * in the current state information, and restore the locations of the pointers
 288  * from the old state information.  This is done by multiplexing the pointer
 289  * location into the zeroeth word of the state information.
 290  *
 291  * Note that due to the order in which things are done, it is OK to call
 292  * setstate() with the same state as the current state.
 293  *
 294  * Returns a pointer to the old state information.
 295  */
 296 char *
 297 sim_setstate(arg_state)
 298         char *arg_state;
 299 {
 300         register QUAD *new_state = (QUAD *)arg_state;
 301         register int type = new_state[0] % MAX_TYPES;
 302         register int rear = new_state[0] / MAX_TYPES;
 303         char *ostate = (char *)(&state[-1]);
 304
 305         if (rand_type == TYPE_0)
 306                 state[-1] = rand_type;
 307         else
 308                 state[-1] = MAX_TYPES * (rptr - state) + rand_type;
 309         switch(type) {
 310         case TYPE_0:
 311         case TYPE_1:
 312         case TYPE_2:
 313         case TYPE_3:
 314         case TYPE_4:
 315                 rand_type = type;
 316                 rand_deg = degrees[type];
 317                 rand_sep = seps[type];
 318                 break;
 319         default:
 320                 (void)fprintf(stderr,
 321                     "random: state info corrupted; not changed.\n");
 322         }
 323         state = &new_state[1];
 324         if (rand_type != TYPE_0) {
 325                 rptr = &state[rear];
 326                 fptr = &state[(rear + rand_sep) % rand_deg];
 327         }
 328         end_ptr = &state[rand_deg];             /* set end_ptr too */
 329         return(ostate);
 330 }
 331
 332 /*
 333  * random:
 334  *
 335  * If we are using the trivial TYPE_0 R.N.G., just do the old linear
 336  * congruential bit.  Otherwise, we do our fancy trinomial stuff, which is
 337  * the same in all the other cases due to all the global variables that have
 338  * been set up.  The basic operation is to add the number at the rear pointer
 339  * into the one at the front pointer.  Then both pointers are advanced to
 340  * the next location cyclically in the table.  The value returned is the sum
 341  * generated, reduced to 31 bits by throwing away the "least random" low bit.
 342  *
 343  * Note: the code takes advantage of the fact that both the front and
 344  * rear pointers can't wrap on the same call by not testing the rear
 345  * pointer if the front one has wrapped.
 346  *
 347  * Returns a 31-bit random number.
 348  */
 349 QUAD
 350 sim_random()
 351 {
 352         QUAD i;
 353
 354         if (rand_type == TYPE_0)
 355                 i = state[0] = (state[0] * 1103515245 + 12345) & 0x7fffffff;
 356         else {
 357                 *fptr += *rptr;
 358                 i = (*fptr >> 1) & 0x7fffffff;  /* chucking least random bit */
 359                 if (++fptr >= end_ptr) {
 360                         fptr = state;
 361                         ++rptr;
 362                 } else if (++rptr >= end_ptr)
 363                         rptr = state;
 364         }
 365         return(i);
 366 }