#ifndef lint #ifdef __GNUC__ #define ATTRIBUTE(attrs) __attribute__(attrs) #else #define ATTRIBUTE(attrs) #endif static char Rcs_Id[] ATTRIBUTE((used)) = "$Id: mtwist.c,v 1.28 2014-01-23 21:11:42-08 geoff Exp $"; #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.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html * * These functions were written by Geoff Kuenning, Claremont, CA. * * IMPORTANT NOTE: this implementation assumes a modern compiler. In * particular, it assumes that the "inline" keyword is available, and * that the "inttypes.h" header file is present. * * IMPORTANT NOTE: this software requires access to a 32-bit type. * The Mersenne Twist algorithms are not guaranteed to produce correct * results with a 64-bit type. * * 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.28 2014-01-23 21:11:42-08 geoff * Fix an obsolete gettimeofday call * * Revision 1.27 2013-06-12 23:22:03-07 geoff * Validity-check state pointer when saving. * * Revision 1.26 2013-01-01 01:18:52-08 geoff * Fix a lot of comiler warnings. Try to fall back to /dev/random if * /dev/urandom doesn't exist. * * Revision 1.25 2012-12-30 16:24:49-08 geoff * Use gcc attributes to suppress warnings on Rcs_Id. Also get rid of * most of the table of contents, to suppress a different gcc warning. * Fix mts_seed and mts_goodseed to return the 32-bit seeds they chose. * Fix all three /dev/Xrandom seeding functions to use only the entropy * they need, rather than letting stdio eat far more than necessary. * (Thanks to Markus Armbruster for all three ideas.) * * Revision 1.24 2011-02-18 00:49:12-08 geoff * Fix a (harmless) typecasting error in mts_bestseed. Thanks to Waclaw * Kusnierczyk for pointing it out. * * Revision 1.23 2010-12-11 00:28:18+13 geoff * Add support for GENERATE_CODE_IN_HEADER. Fix the URL for the original * Web page. * * Revision 1.22 2010-06-24 20:53:59+12 geoff * Switch to using types and formats from inttypes.h. Get rid of all * compilation options. * * Revision 1.21 2010-06-24 00:29:38-07 geoff * Correctly save and restore the state vector even if longs are wider * than 32 bits. * * Revision 1.20 2007-10-26 00:21:06-07 geoff * Use the new mt_u32bit_t type. * * 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 #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. */ #define MT_INLINE /* Disable the inline keyword */ #define MT_EXTERN /* Generate real code for functions */ #undef MT_GENERATE_CODE_IN_HEADER #define MT_GENERATE_CODE_IN_HEADER 1 /* Generate code when #including */ #include "mtwist.h" /* * Table of (non-global) contents: */ static uint32_t mts_devseed(mt_state* state, char* seed_dev); /* Choose seed from a device */ /* * 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 an unsigned int. Don't change this even if your ints 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 uint32_t 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 */ uint32_t 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 * ints. */ 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 */ uint32_t seed) /* 32-bit seed to start from */ { int i; /* Loop index */ uint32_t 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 * ints. */ 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-int 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 */ uint32_t 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. */ uint32_t mts_seed( mt_state* state) /* State vector to seed */ { return 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. */ uint32_t mts_goodseed( mt_state* state) /* State vector to seed */ { return 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 uint32_t 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 (uint32_t)]; /* Space for reading random int */ uint32_t 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 */ #endif /* WIN32 */ ranfile = fopen(seed_dev, "rb"); /* * Some machines have /dev/random but not /dev/urandom. On those * machines, /dev/random is nonblocking, so we'll try it before we * fall back to using the time. */ if (ranfile == NULL) ranfile = fopen(DEVRANDOM, "rb"); if (ranfile != NULL) { setbuf(ranfile, NULL); for (nextbyte = 0; nextbyte < (int)sizeof randomunion.ranbuffer; nextbyte += bytesread) { bytesread = fread(&randomunion.ranbuffer[nextbyte], (size_t)1, sizeof randomunion.ranbuffer - nextbyte, ranfile); if (bytesread == 0) break; } fclose(ranfile); if (nextbyte == sizeof randomunion.ranbuffer) { mts_seed32new(state, randomunion.randomvalue); return randomunion.randomvalue; } } /* * 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, NULL); #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); return 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-int 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[0] + nextbyte, (size_t)1, sizeof state->statevec - nextbyte, ranfile); if (bytesread == 0) { /* * Something went wrong. Fall back to time-based seeding. */ fclose(ranfile); mts_goodseed(state); return; } } fclose(ranfile); } /* * 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 uint32_t* state_ptr; /* Next place to get from state */ register uint32_t value1; /* Scratch val picked up from state */ register uint32_t 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); /* * Ensure the state pointer is valid. */ if (state->stateptr < 0 || state->stateptr > MT_STATE_SIZE) { fprintf(stderr, "Mtwist internal: Trying to write invalid state pointer %d\n", state->stateptr); mts_refresh(state); } for (i = MT_STATE_SIZE; --i >= 0; ) { if (fprintf(statefile, "%" PRIu32 " ", 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, "%" SCNu32, &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( uint32_t 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( uint32_t seed) /* 32-bit seed to start from */ { mts_seed32new(&mt_default_state, seed); } /* * Initialize a Mersenne Twist RNG from a 624-int seed. * * See mts_seedfull for full commentary. */ void mt_seedfull( uint32_t seeds[MT_STATE_SIZE]) { mts_seedfull(&mt_default_state, seeds); } /* * Initialize the PRNG from random input. See mts_seed. */ uint32_t mt_seed(void) { return mts_seed(&mt_default_state); } /* * Initialize the PRNG from random input. See mts_goodseed. */ uint32_t mt_goodseed(void) { return mts_goodseed(&mt_default_state); } /* * Initialize the PRNG from random input. See mts_bestseed. */ void mt_bestseed(void) { 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. */ mt_state* mt_getstate(void) { 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); }