Back to index

plt-scheme  4.2.1
random.inc
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1983, 1993
00003  *     The Regents of the University of California.  All rights reserved.
00004  *
00005  * Redistribution and use in source and binary forms, with or without
00006  * modification, are permitted provided that the following conditions
00007  * are met:
00008  * 1. Redistributions of source code must retain the above copyright
00009  *    notice, this list of conditions and the following disclaimer.
00010  * 2. Redistributions in binary form must reproduce the above copyright
00011  *    notice, this list of conditions and the following disclaimer in the
00012  *    documentation and/or other materials provided with the distribution.
00013  * 3. All advertising materials mentioning features or use of this software
00014  *    must display the following acknowledgement:
00015  *     This product includes software developed by the University of
00016  *     California, Berkeley and its contributors.
00017  * 4. Neither the name of the University nor the names of its contributors
00018  *    may be used to endorse or promote products derived from this software
00019  *    without specific prior written permission.
00020  *
00021  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
00022  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00023  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00024  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
00025  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
00026  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
00027  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00028  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00029  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
00030  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
00031  * SUCH DAMAGE.
00032  */
00033 
00034 /* This is a stripped-down version of random.c distributed with
00035    FreeBSD 2.2 */
00036 
00037 /*
00038  * random.c:
00039  *
00040  * An improved random number generation package.  In addition to the standard
00041  * rand()/srand() like interface, this package also has a special state info
00042  * interface.  The initstate() routine is called with a seed, an array of
00043  * bytes, and a count of how many bytes are being passed in; this array is
00044  * then initialized to contain information for random number generation with
00045  * that much state information.  Good sizes for the amount of state
00046  * information are 32, 64, 128, and 256 bytes.  The state can be switched by
00047  * calling the setstate() routine with the same array as was initiallized
00048  * with initstate().  By default, the package runs with 128 bytes of state
00049  * information and generates far better random numbers than a linear
00050  * congruential generator.  If the amount of state information is less than
00051  * 32 bytes, a simple linear congruential R.N.G. is used.
00052  *
00053  * Internally, the state information is treated as an array of longs; the
00054  * zeroeth element of the array is the type of R.N.G. being used (small
00055  * integer); the remainder of the array is the state information for the
00056  * R.N.G.  Thus, 32 bytes of state information will give 7 longs worth of
00057  * state information, which will allow a degree seven polynomial.  (Note:
00058  * the zeroeth word of state information also has some other information
00059  * stored in it -- see setstate() for details).
00060  *
00061  * The random number generation technique is a linear feedback shift register
00062  * approach, employing trinomials (since there are fewer terms to sum up that
00063  * way).  In this approach, the least significant bit of all the numbers in
00064  * the state table will act as a linear feedback shift register, and will
00065  * have period 2^deg - 1 (where deg is the degree of the polynomial being
00066  * used, assuming that the polynomial is irreducible and primitive).  The
00067  * higher order bits will have longer periods, since their values are also
00068  * influenced by pseudo-random carries out of the lower bits.  The total
00069  * period of the generator is approximately deg*(2**deg - 1); thus doubling
00070  * the amount of state information has a vast influence on the period of the
00071  * generator.  Note: the deg*(2**deg - 1) is an approximation only good for
00072  * large deg, when the period of the shift register is the dominant factor.
00073  * With deg equal to seven, the period is actually much longer than the
00074  * 7*(2**7 - 1) predicted by this formula.
00075  */
00076 
00077 /*
00078  * For each of the currently supported random number generators, we have a
00079  * break value on the amount of state information (you need at least this
00080  * many bytes of state info to support this random number generator), a degree
00081  * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
00082  * the separation between the two lower order coefficients of the trinomial.
00083  */
00084 /* x**31 + x**3 + 1 */
00085 #define       DEG           MZ_RANDOM_STATE_DEG
00086 #define       SEP           3
00087 
00088 /*
00089  * fpos and rpos are two pointers into the state info, a front and a rear
00090  * pointer.  These two pointers are always SEP places aparts, as they
00091  * cycle cyclically through the state information.  (Yes, this does mean we
00092  * could get away with just one pointer, but the code for random() is more
00093  * efficient this way).  The pointers are left positioned as they would be
00094  * from the call
00095  *
00096  *     initstate(1, randtbl, 128);
00097  *
00098  * (The position of the rear pointer, rptr, is really 0 (as explained above
00099  * in the initialization of randtbl) because the state table pointer is set
00100  * to point to randtbl[1] (as explained below).
00101  */
00102 #if 0
00103 /* In schpriv.h: */
00104 typedef struct {
00105   Scheme_Type type;
00106   short fpos, rpos;
00107   long state[DEG];
00108 } Scheme_Random_State;
00109 #endif
00110 
00111 /*
00112  * sch_rand:
00113  *
00114  * The basic operation is to add the number at the rear pointer
00115  * into the one at the front pointer.  Then both pointers are advanced to
00116  * the next location cyclically in the table.  The value returned is the sum
00117  * generated, reduced to 31 bits by throwing away the "least random" low bit.
00118  *
00119  * Note: the code takes advantage of the fact that both the front and
00120  * rear pointers can't wrap on the same call by not testing the rear
00121  * pointer if the front one has wrapped.
00122  *
00123  * Returns a 31-bit random number.
00124  */
00125 long scheme_rand(unsigned long n, Scheme_Random_State *rs)
00126 {
00127   long i;
00128   
00129   rs->state[rs->fpos] += rs->state[rs->rpos];
00130   i = (rs->state[rs->fpos] >> 1) & 0x7fffffff; /* chucking least random bit */
00131   if (++(rs->fpos) >= DEG) {
00132     rs->fpos = 0;
00133     rs->rpos++;
00134   } else if (++(rs->rpos) >= DEG)
00135     rs->rpos = 0;
00136 
00137   return i % n;
00138 }
00139 
00140 /*
00141  * sch_srand:
00142  *
00143  * Initialize the random number generator based on the given seed.  If the
00144  * type is the trivial no-state-information type, just remember the seed.
00145  * Otherwise, initializes state[] based on the given "seed" via a linear
00146  * congruential generator.  Then, the pointers are set to known locations
00147  * that are exactly SEP places apart.  Lastly, it cycles the state
00148  * information a given number of times to get rid of any initial dependencies
00149  * introduced by the L.C.R.N.G.  Note that the initialization of randtbl[]
00150  * for default usage relies on values produced by this routine.
00151  */
00152 static void sch_srand(long x, Scheme_Random_State *rs, int reset)
00153 {
00154   register int i;
00155   long *state = rs->state;
00156 
00157   state[0] = x & 0x7fffffff;
00158   for (i = 1; i < DEG; i++) {
00159     /*
00160      * Compute v = (7^5 * v) mod (2^31 - 1)
00161      * wihout overflowing 31 bits:
00162      *      (2^31 - 1) = 127773 * (7^5) + 2836
00163      * From "Random number generators: good ones are hard to find",
00164      * Park and Miller, Communications of the ACM, vol. 31, no. 10,
00165      * October 1988, p. 1195.
00166      */
00167     register long hi, lo, v;
00168   
00169     v = state[i - 1];
00170     hi = v / 127773;
00171     lo = v % 127773;
00172     v = 16807 * lo - 2836 * hi;
00173     if (v <= 0)
00174       v += 0x7fffffff;
00175     state[i] = v;
00176   }
00177   rs->fpos = SEP;
00178   rs->rpos = 0;
00179   for (i = 0; i < 10 * DEG; i++) {
00180     (void)scheme_rand(rs);
00181 }
00182 }
00183 
00184 Scheme_Object *scheme_make_random_state(long seed)
00185 {
00186   Scheme_Random_State *rs;
00187 
00188   rs = MALLOC_ONE_TAGGED(Scheme_Random_State);
00189   rs->so.type = scheme_random_state_type;
00190 
00191   sch_srand(seed, rs);
00192 
00193   return (Scheme_Object *)rs;
00194 }