#ifndef lint static char Rcs_Id[] = "$Id: mtwist.c,v 1.1 2004/03/23 09:30:04 geoff Exp geoff $"; #endif /* * C library functions for generating pseudorandom numbers using the * Mersenne Twist algorithm. See M. Matsumoto and T. Nishimura, * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform * Pseudo-Random Number Generator", ACM Transactions on Modeling and * Computer Simulation, Vol. 8, No. 1, January 1998, pp 3--30. * * The Web page on the Mersenne Twist algorithm is at: * * http://www.math.keio.ac.jp/~matumoto/emt.html * * These functions were written by Geoffrey H. Kuenning, Claremont, CA. * * IMPORTANT NOTE: the Makefile must define two machine-specific * variables to get optimum features and performance: * * MT_NO_INLINE should be defined if the compiler doesn't support * the "inline" keyword. * MT_NO_LONGLONG should be defined if the compiler doesn't support a * "long long" type for 64-bit integers * MT_MACHINE_BITS must be either 32 or 64, reflecting the natural * size of the processor registers. If undefined, it * will default to 32. * * The first two variables above are defined in an inverted sense * because I expect that most compilers will support inline and * long-long. By inverting the sense, this common case will require * no special compiler flags. * * IMPORTANT NOTE: this software assumes that the inherent width of a * "long" is 32 bits. If you are running on a machine that uses * 64-bit longs, some of the declarations and code will have to be * modified. * * This software is based on LGPL-ed code by Takuji Nishimura. It has * also been heavily influenced by code written by Shawn Cokus, and * somewhat influenced by code written by Richard J. Wagner. It is * therefore also distributed under the LGPL: * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public License * as published by the Free Software Foundation; either version 2 of * the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Library General Public License for more details. You should have * received a copy of the GNU Library General Public License along * with this library; if not, write to the Free Foundation, Inc., 59 * Temple Place, Suite 330, Boston, MA 02111-1307 USA * * $Log: mtwist.c,v $ * Revision 1.1 2004/03/23 09:30:04 geoff * New addition * * Revision 1.19 2003/09/11 05:55:19 geoff * Get rid of some minor compiler warnings. * * Revision 1.18 2003/09/11 05:50:53 geoff * Don't #define inline to nothing, since that breaks standard include * files. Instead, use MT_INLINE as a synonym. * * Revision 1.17 2002/10/31 22:07:10 geoff * Make WIN32 detection work with GCC as well as MS C * * Revision 1.16 2002/10/31 22:04:59 geoff * Fix a typo in the WIN32 option * * Revision 1.15 2002/10/31 06:01:43 geoff * Incorporate Joseph Brill's Windows-portability changes * * Revision 1.14 2002/10/30 07:39:53 geoff * Reintroduce the old seeding functions (so that old code will still * produce the same results), and give the new versions new names. * * Revision 1.13 2002/10/30 01:08:26 geoff * Switch to M&T's new initialization method * * Revision 1.12 2001/06/18 05:40:12 geoff * Prefix the compile options with MT_. * * Revision 1.11 2001/06/14 10:26:59 geoff * Invert the sense of the #define flags so that the default is the * normal case (if gcc is normal!). Also default MT_MACHINE_BITS to 32. * * Revision 1.10 2001/06/14 10:10:38 geoff * Move the RNG functions into the header file so they can be inlined. * Add saving/loading of state. Add a function that marks the PRNG as * initialized while also calculating critical constants. Run the * refresh routine whenever seed32 is called. Add functions to seed * based on /dev/random or the time. * * Revision 1.9 2001/06/11 10:00:04 geoff * Major changes to improve flexibility and performance, and to prepare * for inlining. This code is about as fast as it can get without * inlining the various PRNG functions. Add seed/goodseed/bestseed for * seeding from random start values. Add the refresh routine a la Cokus, * but optimize it by unrolling loops. Change getstate to return a * complete state pointer, since knowing the position in the state vector * is critical to restoring state. Add more macros to improve * readability. Rename certain macros in preparation for inlining. Get * rid of leftover optimizer-bug stuff. Stop using mtwist_guts.h; * instead use direct code (via macros) and the refresh function. * * Revision 1.8 2001/04/23 08:36:03 geoff * Move the #defined code into a header file to ease stepping with a debugger. * * Revision 1.7 2001/04/23 08:00:13 geoff * Add code to work around optimizer bug * * Revision 1.6 2001/04/14 01:33:32 geoff * Clarify the license * * Revision 1.5 2001/04/09 08:45:00 geoff * Rename default_state to mt_default_state, and make it global so that * the random-distribution code can use it. * * Revision 1.4 2001/04/07 23:24:11 geoff * My guess in the commentary for the last delta was right: it's faster * on a x86 to convert the two halves of the PRN to double, multiplying * them by the appropriate value to scale them, and then add them as * doubles. I suspect the reason is that there is no instruction to * convert a 64-bit value directly to a double, so the work of building * the long long (which isn't easy anyway, without assembly access) is * worse than wasted. So add support for MT_MACHINE_BITS, and only go * the via-long-long route on a true 64-bit machine. * * Revision 1.3 2001/04/07 23:09:38 geoff * Get rid of MT_INLINE. Convert all of the code to use preprocessor * macros for the guts of the PRNG code. Take advantage of the * conversion to get rid of unnecessary calls initialization tests. Also * clean up the generation of long-double pseudorandom numbers on * machines that have the long long type (by converting first to a long * long, then to a double, saving one floating-point operation). The * latter change might be a mistake on 32-bit machines. The code is now * much faster as a result of macro-izing. * * Revision 1.2 2001/04/07 22:21:41 geoff * Make the long-double code a hair faster by always having a 64-bit * conversion constant. Add commentary to the PRNG loop. * * Revision 1.1 2001/04/07 09:43:41 geoff * Initial revision * */ #ifdef _WIN32 #undef WIN32 #define WIN32 #endif /* _WIN32 */ #include #include #ifdef WIN32 #include #else /* WIN32 */ #include #endif /* WIN32 */ /* * Before we include the Mersenne Twist header file, we must do a bit * of magic setup. The code for actual random-number generation * resides in that file rather than here. We need to arrange for the * code to be compiled into this .o file, either because inlines * aren't supported or because somebody might want to take a pointer * to a function. We do so with a couple of careful #defines. */ #undef MT_NO_INLINE /* Ask for code to be compiled */ #define MT_INLINE /* Disable the inline keyword */ #define MT_EXTERN /* Generate real code for functions */ #include "mtwist.h" /* * Table of contents: */ void mts_mark_initialized(mt_state* state); /* Mark a PRNG state as initialized */ void mts_seed32(mt_state* state, unsigned long seed); /* Set random seed for any generator */ void mts_seed32new(mt_state* state, unsigned long seed); /* Set random seed for any generator */ void mts_seedfull(mt_state* state, unsigned long seeds[MT_STATE_SIZE]); /* Set complicated seed for any gen. */ void mts_seed(mt_state* state); /* Choose seed from random input */ void mts_goodseed(mt_state* state); /* Choose seed from more random */ /* ..input than mts_seed */ static void mts_devseed(mt_state* state, char* seed_dev); /* Choose seed from a device */ void mts_bestseed(mt_state* state); /* Choose seed from extremely random */ /* ..input (can be *very* slow) */ void mts_refresh(mt_state* state); /* Generate 624 more random values */ int mts_savestate(FILE* statefile, mt_state* state); /* Save state to a file (ASCII) */ int mts_loadstate(FILE* statefile, mt_state* state); /* Load state from a file (ASCII) */ void mt_seed32(unsigned long seed); /* Set random seed for default gen. */ void mt_seed32new(unsigned long seed); /* Set random seed for default gen. */ void mt_seedfull(unsigned long seeds[MT_STATE_SIZE]); /* Set complicated seed for default */ void mt_seed(void); /* Choose seed from random input */ void mt_goodseed(void); /* Choose seed from more random */ /* ..input than mts_seed */ void mt_bestseed(void); /* Choose seed from extremely random */ /* ..input (can be *very* slow) */ extern mt_state* mt_getstate(void); /* Get current state of default */ /* ..generator */ int mt_savestate(FILE* statefile); /* Save state to a file (ASCII) */ int mt_loadstate(FILE* statefile); /* Load state from a file (ASCII) */ /* * The following values are fundamental parameters of the algorithm. * With the exception of the two masks, all of them were found * experimentally using methods described in Matsumoto and Nishimura's * paper. They are exceedingly magic; don't change them. */ /* MT_STATE_SIZE is defined in the header file. */ #define RECURRENCE_OFFSET 397 /* Offset into state space for the */ /* ..recurrence relation. The */ /* ..recurrence mashes together two */ /* ..values that are separated by */ /* ..this offset in the state */ /* ..space. */ #define MATRIX_A 0x9908b0df /* Constant vector A for the */ /* ..recurrence relation. The */ /* ..mashed-together value is */ /* ..multiplied by this vector to */ /* ..get a new value that will be */ /* ..stored into the state space. */ /* * Width of a long. Don't change this even if your longs are 64 bits. */ #define BIT_WIDTH 32 /* Work with 32-bit words */ /* * Masks for extracting the bits to be mashed together. The widths of these * masks are also fundamental parameters of the algorithm, determined * experimentally -- but of course the masks themselves are simply bit * selectors. */ #define UPPER_MASK 0x80000000 /* Most significant w-r bits */ #define LOWER_MASK 0x7fffffff /* Least significant r bits */ /* * Macro to simplify code in the generation loop. This function * combines the top bit of x with the bottom 31 bits of y. */ #define COMBINE_BITS(x, y) \ (((x) & UPPER_MASK) | ((y) & LOWER_MASK)) /* * Another generation-simplification macro. This one does the magic * scrambling function. */ #define MATRIX_MULTIPLY(original, new) \ ((original) ^ ((new) >> 1) \ ^ matrix_decider[(new) & 0x1]) /* * Parameters of Knuth's PRNG (Line 25, Table 1, p. 102 of "The Art of * Computer Programming, Vol. 2, 2nd ed, 1981). */ #define KNUTH_MULTIPLIER_OLD \ 69069 /* * Parameters of Knuth's PRNG (p. 106 of "The Art of Computer * Programming, Vol. 2, 3rd ed). */ #define KNUTH_MULTIPLIER_NEW \ 1812433253ul #define KNUTH_SHIFT 30 // Even on a 64-bit machine! /* * Default 32-bit random seed if mts_seed32 wasn't called */ #define DEFAULT_SEED32_OLD \ 4357 #define DEFAULT_SEED32_NEW \ 5489ul /* * Where to get random numbers */ #define DEVRANDOM "/dev/random" #define DEVURANDOM "/dev/urandom" /* * Many applications need only a single PRNG, so it's a nuisance to have to * specify a state. For those applications, we will provide a default * state, and functions to use it. */ mt_state mt_default_state; /* * To generate double-precision random numbers, we need to divide the result * of mts_lrand or mts_llrand by 2^32 or 2^64, respectively. The quickest * way to do that on most machines is to multiply by the inverses of those * numbers. However, I don't trust the compiler to correctly convert the * corresponding decimal constant. So we will compute the correct number at * run time as part of initialization, which will produce a nice exact * result. */ double mt_32_to_double; /* Multiplier to convert long to dbl */ double mt_64_to_double; /* Mult'r to cvt long long to dbl */ /* * In the recurrence relation, the new value is XORed with MATRIX_A only if * the lower bit is nonzero. Since most modern machines don't like to * branch, it's vastly faster to handle this decision by indexing into an * array. The chosen bit is used as an index into the following vector, * which produces either zero or MATRIX_A and thus the desired effect. */ static unsigned long matrix_decider[2] = {0x0, MATRIX_A}; /* * Mark a PRNG's state as having been initialized. This is the only * way to set that field nonzero; that way we can be sure that the * constants are set properly before the PRNG is used. * * As a side effect, set up some constants that the PRNG assumes are * valid. These are calculated at initialization time rather than * being written as decimal constants because I frankly don't trust * the compiler's ASCII conversion routines. */ void mts_mark_initialized( mt_state* state) /* State vector to mark initialized */ { int i; /* Power of 2 being calculated */ /* * Figure out the proper multiplier for long-to-double conversion. We * don't worry too much about efficiency, since the assumption is that * initialization is vastly rarer than generation of random numbers. */ mt_32_to_double = 1.0; for (i = 0; i < BIT_WIDTH; i++) mt_32_to_double /= 2.0; mt_64_to_double = mt_32_to_double; for (i = 0; i < BIT_WIDTH; i++) mt_64_to_double /= 2.0; state->initialized = 1; } /* * Initialize a Mersenne Twist PRNG from a 32-bit seed. * * According to Matsumoto and Nishimura's paper, the seed array needs to be * filled with nonzero values. (My own interpretation is that there needs * to be at least one nonzero value). They suggest using Knuth's PRNG from * Line 25, Table 1, p.102, "The Art of Computer Programming," Vol. 2 (2nd * ed.), 1981. I find that rather odd, since that particular PRNG is * sensitive to having an initial seed of zero (there are many other PRNGs * out there that have an additive component, so that a seed of zero does * not generate a repeating-zero sequence). However, one thing I learned * from reading Knuth is that you shouldn't second-guess mathematicians * about PRNGs. Also, by following M & N's approach, we will be compatible * with other implementations. So I'm going to stick with their version, * with the single addition that a zero seed will be changed to their * default seed. */ void mts_seed32( mt_state* state, /* State vector to initialize */ unsigned long seed) /* 32-bit seed to start from */ { int i; /* Loop index */ if (seed == 0) seed = DEFAULT_SEED32_OLD; /* * Fill the state vector using Knuth's PRNG. Be sure to mask down * to 32 bits in case we're running on a machine with 64-bit * longs. */ state->statevec[MT_STATE_SIZE - 1] = seed & 0xffffffff; for (i = MT_STATE_SIZE - 2; i >= 0; i--) state->statevec[i] = (KNUTH_MULTIPLIER_OLD * state->statevec[i + 1]) & 0xffffffff; state->stateptr = MT_STATE_SIZE; mts_mark_initialized(state); /* * Matsumoto and Nishimura's implementation refreshes the PRNG * immediately after running the Knuth algorithm. This is * probably a good thing, since Knuth's PRNG doesn't generate very * good numbers. */ mts_refresh(state); } /* * Initialize a Mersenne Twist PRNG from a 32-bit seed, using * Matsumoto and Nishimura's newer reference implementation (Jan. 9, * 2002). */ void mts_seed32new( mt_state* state, /* State vector to initialize */ unsigned long seed) /* 32-bit seed to start from */ { int i; /* Loop index */ unsigned long nextval; /* Next value being calculated */ /* * Fill the state vector using Knuth's PRNG. Be sure to mask down * to 32 bits in case we're running on a machine with 64-bit * longs. */ state->statevec[MT_STATE_SIZE - 1] = seed & 0xffffffffUL; for (i = MT_STATE_SIZE - 2; i >= 0; i--) { nextval = state->statevec[i + 1] >> KNUTH_SHIFT; nextval ^= state->statevec[i + 1]; nextval *= KNUTH_MULTIPLIER_NEW; nextval += (MT_STATE_SIZE - 1) - i; state->statevec[i] = nextval & 0xffffffffUL; } state->stateptr = MT_STATE_SIZE; mts_mark_initialized(state); /* * Matsumoto and Nishimura's implementation refreshes the PRNG * immediately after running the Knuth algorithm. This is * probably a good thing, since Knuth's PRNG doesn't generate very * good numbers. */ mts_refresh(state); } /* * Initialize a Mersenne Twist RNG from a 624-long seed. * * The 32-bit seeding routine given by Matsumoto and Nishimura has the * drawback that there are only 2^32 different PRNG sequences that can be * generated by calling that function. This function solves that problem by * allowing a full 624*32-bit state to be given. (Note that 31 bits of the * given state are ignored; see the paper for details.) * * Since an all-zero state would cause the PRNG to cycle, we detect * that case and abort the program (silently, since there is no * portable way to produce a message in both C and C++ environments). * An alternative would be to artificially force the state to some * known nonzero value. However, I feel that if the user is providing * a full state, it's a bug to provide all zeros and we we shouldn't * conceal the bug by generating apparently correct output. */ void mts_seedfull( mt_state* state, /* State vector to initialize */ unsigned long seeds[MT_STATE_SIZE]) /* Seed array to start from */ { int had_nz = 0; /* NZ if at least one NZ seen */ int i; /* Loop index */ for (i = 0; i < MT_STATE_SIZE; i++) { if (seeds[i] != 0) had_nz = 1; state->statevec[MT_STATE_SIZE - i - 1] = seeds[i]; } if (!had_nz) { /* * It would be nice to abort with a message. Unfortunately, fprintf * isn't compatible with all implementations of C++. In the * interest of C++ compatibility, therefore, we will simply abort * silently. It will unfortunately be up to a programmer to run * under a debugger (or examine the core dump) to discover the cause * of the abort. */ abort(); } state->stateptr = MT_STATE_SIZE; mts_mark_initialized(state); } /* * Choose a seed based on some moderately random input. Prefers * /dev/urandom as a source of random numbers, but uses the lower bits * of the current time if /dev/urandom is not available. In any case, * only provides 32 bits of entropy. */ void mts_seed( mt_state* state) /* State vector to seed */ { mts_devseed(state, DEVURANDOM); } /* * Choose a seed based on some fairly random input. Prefers * /dev/random as a source of random numbers, but uses the lower bits * of the current time if /dev/random is not available. In any case, * only provides 32 bits of entropy. */ void mts_goodseed( mt_state* state) /* State vector to seed */ { mts_devseed(state, DEVRANDOM); } /* * Choose a seed based on a random-number device given by the caller. * If that device can't be opened, use the lower 32 bits from the * current time. */ static void mts_devseed( mt_state* state, /* State vector to seed */ char* seed_dev) /* Device to seed from */ { int bytesread; /* Byte count read from device */ int nextbyte; /* Index of next byte to read */ FILE* ranfile; /* Access to device */ union { char ranbuffer[sizeof (unsigned long)]; /* Space for reading random int */ unsigned long randomvalue; /* Random value for initialization */ } randomunion; /* Union for reading random int */ #ifdef WIN32 struct _timeb tb; /* Time of day (Windows mode) */ #else /* WIN32 */ struct timeval tv; /* Time of day */ struct timezone tz; /* Dummy for gettimeofday */ #endif /* WIN32 */ ranfile = fopen(seed_dev, "rb"); if (ranfile != NULL) { for (nextbyte = 0; nextbyte < (int)sizeof randomunion.ranbuffer; nextbyte += bytesread) { bytesread = fread(&randomunion.ranbuffer[nextbyte], 1, sizeof randomunion.ranbuffer - nextbyte, ranfile); if (bytesread == 0) break; } fclose(ranfile); if (nextbyte == sizeof randomunion.ranbuffer) { mts_seed32new(state, randomunion.randomvalue); return; } } /* * The device isn't available. Use the time. We will * assume that the time of day is accurate to microsecond * resolution, which is true on most modern machines. */ #ifdef WIN32 (void) _ftime (&tb); #else /* WIN32 */ (void) gettimeofday (&tv, &tz); #endif /* WIN32 */ /* * We just let the excess part of the seconds field overflow */ #ifdef WIN32 randomunion.randomvalue = tb.time * 1000 + tb.millitm; #else /* WIN32 */ randomunion.randomvalue = tv.tv_sec * 1000000 + tv.tv_usec; #endif /* WIN32 */ mts_seed32new(state, randomunion.randomvalue); } /* * Choose a seed based on the best random input available. Prefers * /dev/random as a source of random numbers, and reads the entire * 624-long state from that device. Because of this approach, the * function can take a long time (in real time) to complete, since * /dev/random may have to wait quite a while before it can provide * that much randomness. If /dev/random is unavailable, falls back to * calling mts_goodseed. */ void mts_bestseed( mt_state* state) /* State vector to seed */ { int bytesread; /* Byte count read from device */ int nextbyte; /* Index of next byte to read */ FILE* ranfile; /* Access to device */ ranfile = fopen("/dev/random", "rb"); if (ranfile == NULL) { mts_goodseed(state); return; } for (nextbyte = 0; nextbyte < (int)sizeof state->statevec; nextbyte += bytesread) { bytesread = fread((char *)&state->statevec + nextbyte, 1, sizeof state->statevec - nextbyte, ranfile); if (bytesread == 0) { /* * Something went wrong. Fall back to time-based seeding. */ fclose(ranfile); mts_goodseed(state); return; } } } /* * Generate 624 more random values. This function is called when the * state vector has been exhausted. It generates another batch of * pseudo-random values. The performance of this function is critical * to the performance of the Mersenne Twist PRNG, so it has been * highly optimized. */ void mts_refresh( register mt_state* state) /* State for the PRNG */ { register int i; /* Index into the state */ register unsigned long* state_ptr; /* Next place to get from state */ register unsigned long value1; /* Scratch val picked up from state */ register unsigned long value2; /* Scratch val picked up from state */ /* * Start by making sure a random seed has been set. If not, set * one. */ if (!state->initialized) { mts_seed32(state, DEFAULT_SEED32_OLD); return; /* Seed32 calls us recursively */ } /* * Now generate the new pseudorandom values by applying the * recurrence relation. We use two loops and a final * 2-statement sequence so that we can handle the wraparound * explicitly, rather than having to use the relatively slow * modulus operator. * * In essence, the recurrence relation concatenates bits * chosen from the current random value (last time around) * with the immediately preceding one. Then it * matrix-multiplies the concatenated bits with a value * RECURRENCE_OFFSET away and a constant matrix. The matrix * multiplication reduces to a shift and two XORs. * * Some comments on the optimizations are in order: * * Strictly speaking, none of the optimizations should be * necessary. All could conceivably be done by a really good * compiler. However, the compilers available to me aren't quite * smart enough, so hand optimization needs to be done. * * Shawn Cokus was the first to achieve a major speedup. In the * original code, the first value given to COMBINE_BITS (in my * characterization) was re-fetched from the state array, rather * than being carried in a scratch variable. Cokus noticed that * the first argument to COMBINE_BITS could be saved in a register * in the previous loop iteration, getting rid of the need for an * expensive memory reference. * * Cokus also switched to using pointers to access the state * array and broke the original loop into two so that he could * avoid using the expensive modulus operator. Cokus used three * pointers; Richard J. Wagner noticed that the offsets between * the three were constant, so that they could be collapsed into a * single pointer and constant-offset accesses. This is clearly * faster on x86 architectures, and is the same cost on RISC * machines. A secondary benefit is that Cokus' version was * register-starved on the x86, while Wagner's version was not. * * I made several smaller improvements to these observations. * First, I reversed the contents of the state vector. In the * current version of the code, this change doesn't directly * affect the performance of the refresh loop, but it has the nice * side benefit that an all-zero state structure represents an * uninitialized generator. It also slightly speeds up the * random-number routines, since they can compare the state * pointer against zero instead of against a constant (this makes * the biggest difference on RISC machines). * * Second, I returned to Matsumoto and Nishimura's original * technique of using a lookup table to decide whether to xor the * constant vector A (MATRIX_A in this code) with the newly * computed value. Cokus and Wagner had used the ?: operator, * which requires a test and branch. Modern machines don't like * branches, so the table lookup is faster. * * Third, in the Cokus and Wagner versions the loop ends with a * statement similar to "value1 = value2", which is necessary to * carry the fetched value into the next loop iteration. I * recognized that if the loop were unrolled so that it generates * two values per iteration, a bit of variable renaming would get * rid of that assignment. A nice side effect is that the * overhead of loop control becomes only half as large. * * It is possible to improve the code's performance somewhat * further. In particular, since the second loop's loop count * factors into 2*2*3*3*11, it could be unrolled yet further. * That's easy to do, too: just change the "/ 2" into a division * by whatever factor you choose, and then use cut-and-paste to * duplicate the code in the body. To remove a few more cycles, * fix the code to decrement state_ptr by the unrolling factor, and * adjust the various offsets appropriately. However, the payoff * will be small. At the moment, the x86 version of the loop is * 25 instructions, of which 3 are involved in loop control * (including the decrementing of state_ptr). Further unrolling by * a factor of 2 would thus produce only about a 6% speedup. * * The logical extension of the unrolling * approach would be to remove the loops and create 624 * appropriate copies of the body. However, I think that doing * the latter is a bit excessive! * * I suspect that a superior optimization would be to simplify the * mathematical operations involved in the recurrence relation. * However, I have no idea whether such a simplification is * feasible. */ state_ptr = &state->statevec[MT_STATE_SIZE - 1]; value1 = *state_ptr; for (i = (MT_STATE_SIZE - RECURRENCE_OFFSET) / 2; --i >= 0; ) { state_ptr -= 2; value2 = state_ptr[1]; value1 = COMBINE_BITS(value1, value2); state_ptr[2] = MATRIX_MULTIPLY(state_ptr[-RECURRENCE_OFFSET + 2], value1); value1 = state_ptr[0]; value2 = COMBINE_BITS(value2, value1); state_ptr[1] = MATRIX_MULTIPLY(state_ptr[-RECURRENCE_OFFSET + 1], value2); } value2 = *--state_ptr; value1 = COMBINE_BITS(value1, value2); state_ptr[1] = MATRIX_MULTIPLY(state_ptr[-RECURRENCE_OFFSET + 1], value1); for (i = (RECURRENCE_OFFSET - 1) / 2; --i >= 0; ) { state_ptr -= 2; value1 = state_ptr[1]; value2 = COMBINE_BITS(value2, value1); state_ptr[2] = MATRIX_MULTIPLY(state_ptr[MT_STATE_SIZE - RECURRENCE_OFFSET + 2], value2); value2 = state_ptr[0]; value1 = COMBINE_BITS(value1, value2); state_ptr[1] = MATRIX_MULTIPLY(state_ptr[MT_STATE_SIZE - RECURRENCE_OFFSET + 1], value1); } /* * The final entry in the table requires the "previous" value * to be gotten from the other end of the state vector, so it * must be handled specially. */ value1 = COMBINE_BITS(value2, state->statevec[MT_STATE_SIZE - 1]); *state_ptr = MATRIX_MULTIPLY(state_ptr[MT_STATE_SIZE - RECURRENCE_OFFSET], value1); /* * Now that refresh is complete, reset the state pointer to allow more * pseudorandom values to be fetched from the state array. */ state->stateptr = MT_STATE_SIZE; } /* * Save state to a file. The save format is compatible with Richard * J. Wagner's format, although the details are different. Returns NZ * if the save succeeded. Produces one very long line containing 625 * numbers. */ int mts_savestate( FILE* statefile, /* File to save to */ mt_state* state) /* State to be saved */ { int i; /* Next word to save */ if (!state->initialized) mts_seed32(state, DEFAULT_SEED32_OLD); for (i = MT_STATE_SIZE; --i >= 0; ) { if (fprintf(statefile, "%lu ", state->statevec[i]) < 0) return 0; } if (fprintf(statefile, "%d\n", state->stateptr) < 0) return 0; return 1; } /* * Load state from a file. Returns NZ if the load succeeded. */ int mts_loadstate( FILE* statefile, /* File to load from */ mt_state* state) /* State to be loaded */ { int i; /* Next word to load */ /* * Set the state to "uninitialized" in case the load fails. */ state->initialized = state->stateptr = 0; for (i = MT_STATE_SIZE; --i >= 0; ) { if (fscanf(statefile, "%lu", &state->statevec[i]) != 1) return 0; } if (fscanf(statefile, "%d", &state->stateptr) != 1) return 0; /* * The only validity checking we can do is to insist that the * state pointer be valid. */ if (state->stateptr < 0 || state->stateptr > MT_STATE_SIZE) { state->stateptr = 0; return 0; } mts_mark_initialized(state); return 1; } /* * Initialize the default Mersenne Twist PRNG from a 32-bit seed. * * See mts_seed32 for full commentary. */ void mt_seed32( unsigned long seed) /* 32-bit seed to start from */ { mts_seed32(&mt_default_state, seed); } /* * Initialize the default Mersenne Twist PRNG from a 32-bit seed. * * See mts_seed32new for full commentary. */ void mt_seed32new( unsigned long seed) /* 32-bit seed to start from */ { mts_seed32new(&mt_default_state, seed); } /* * Initialize a Mersenne Twist RNG from a 624-long seed. * * See mts_seedfull for full commentary. */ void mt_seedfull( unsigned long seeds[MT_STATE_SIZE]) { mts_seedfull(&mt_default_state, seeds); } /* * Initialize the PRNG from random input. See mts_seed. */ void mt_seed() { mts_seed(&mt_default_state); } /* * Initialize the PRNG from random input. See mts_goodseed. */ void mt_goodseed() { mts_goodseed(&mt_default_state); } /* * Initialize the PRNG from random input. See mts_bestseed. */ void mt_bestseed() { mts_bestseed(&mt_default_state); } /* * Return a pointer to the current state of the PRNG. The purpose of * this function is to allow the state to be saved for later * restoration. The state should not be modified; instead, it should * be reused later as a parameter to one of the mts_xxx functions. */ extern mt_state* mt_getstate() { return &mt_default_state; } /* * Save state to a file. The save format is compatible with Richard * J. Wagner's format, although the details are different. */ int mt_savestate( FILE* statefile) /* File to save to */ { return mts_savestate(statefile, &mt_default_state); } /* * Load state from a file. */ int mt_loadstate( FILE* statefile) /* File to load from */ { return mts_loadstate(statefile, &mt_default_state); }