00001 /*00002 * Copyright (c) 1983, 199300003 * The Regents of the University of California. All rights reserved.00004 *00005 * Redistribution and use in source and binary forms, with or without00006 * modification, are permitted provided that the following conditions00007 * are met:00008 * 1. Redistributions of source code must retain the above copyright00009 * notice, this list of conditions and the following disclaimer.00010 * 2. Redistributions in binary form must reproduce the above copyright00011 * notice, this list of conditions and the following disclaimer in the00012 * documentation and/or other materials provided with the distribution.00013 * 3. All advertising materials mentioning features or use of this software00014 * must display the following acknowledgement:00015 * This product includes software developed by the University of00016 * California, Berkeley and its contributors.00017 * 4. Neither the name of the University nor the names of its contributors00018 * may be used to endorse or promote products derived from this software00019 * without specific prior written permission.00020 *00021 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND00022 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE00023 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE00024 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE00025 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL00026 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS00027 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)00028 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT00029 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY00030 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF00031 * SUCH DAMAGE.00032 */00033
00034 /* This is a stripped-down version of random.c distributed with00035 FreeBSD 2.2 */00036
00037 /*00038 * random.c:00039 *00040 * An improved random number generation package. In addition to the standard00041 * rand()/srand() like interface, this package also has a special state info00042 * interface. The initstate() routine is called with a seed, an array of00043 * bytes, and a count of how many bytes are being passed in; this array is00044 * then initialized to contain information for random number generation with00045 * that much state information. Good sizes for the amount of state00046 * information are 32, 64, 128, and 256 bytes. The state can be switched by00047 * calling the setstate() routine with the same array as was initiallized00048 * with initstate(). By default, the package runs with 128 bytes of state00049 * information and generates far better random numbers than a linear00050 * congruential generator. If the amount of state information is less than00051 * 32 bytes, a simple linear congruential R.N.G. is used.00052 *00053 * Internally, the state information is treated as an array of longs; the00054 * zeroeth element of the array is the type of R.N.G. being used (small00055 * integer); the remainder of the array is the state information for the00056 * R.N.G. Thus, 32 bytes of state information will give 7 longs worth of00057 * state information, which will allow a degree seven polynomial. (Note:00058 * the zeroeth word of state information also has some other information00059 * stored in it -- see setstate() for details).00060 *00061 * The random number generation technique is a linear feedback shift register00062 * approach, employing trinomials (since there are fewer terms to sum up that00063 * way). In this approach, the least significant bit of all the numbers in00064 * the state table will act as a linear feedback shift register, and will00065 * have period 2^deg - 1 (where deg is the degree of the polynomial being00066 * used, assuming that the polynomial is irreducible and primitive). The00067 * higher order bits will have longer periods, since their values are also00068 * influenced by pseudo-random carries out of the lower bits. The total00069 * period of the generator is approximately deg*(2**deg - 1); thus doubling00070 * the amount of state information has a vast influence on the period of the00071 * generator. Note: the deg*(2**deg - 1) is an approximation only good for00072 * 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 the00074 * 7*(2**7 - 1) predicted by this formula.00075 */00076
00077 /*00078 * For each of the currently supported random number generators, we have a00079 * break value on the amount of state information (you need at least this00080 * many bytes of state info to support this random number generator), a degree00081 * for the polynomial (actually a trinomial) that the R.N.G. is based on, and00082 * the separation between the two lower order coefficients of the trinomial.00083 */00084 /* x**31 + x**3 + 1 */00085 #define DEG MZ_RANDOM_STATE_DEG00086 #define SEP 300087 00088 /*00089 * fpos and rpos are two pointers into the state info, a front and a rear00090 * pointer. These two pointers are always SEP places aparts, as they00091 * cycle cyclically through the state information. (Yes, this does mean we00092 * could get away with just one pointer, but the code for random() is more00093 * efficient this way). The pointers are left positioned as they would be00094 * from the call00095 *00096 * initstate(1, randtbl, 128);00097 *00098 * (The position of the rear pointer, rptr, is really 0 (as explained above00099 * in the initialization of randtbl) because the state table pointer is set00100 * to point to randtbl[1] (as explained below).00101 */00102 #if 000103 /* In schpriv.h: */00104 typedefstruct {
00105 Scheme_Typetype;
00106 short fpos, rpos;
00107 long state[DEG];
00108 } Scheme_Random_State;
00109 #endif00110 00111 /*00112 * sch_rand:00113 *00114 * The basic operation is to add the number at the rear pointer00115 * into the one at the front pointer. Then both pointers are advanced to00116 * the next location cyclically in the table. The value returned is the sum00117 * 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 and00120 * rear pointers can't wrap on the same call by not testing the rear00121 * pointer if the front one has wrapped.00122 *00123 * Returns a 31-bit random number.00124 */00125 longscheme_rand(unsignedlong n, Scheme_Random_State *rs)
00126 {
00127 longi;
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 } elseif (++(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 the00144 * type is the trivial no-state-information type, just remember the seed.00145 * Otherwise, initializes state[] based on the given "seed" via a linear00146 * congruential generator. Then, the pointers are set to known locations00147 * that are exactly SEP places apart. Lastly, it cycles the state00148 * information a given number of times to get rid of any initial dependencies00149 * 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 staticvoid sch_srand(longx, Scheme_Random_State *rs, int reset)
00153 {
00154 registerinti;
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) + 283600163 * 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 registerlong 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 }