Scott David Daniels wrote:
No, I don't really need MT. The others would be fine. I'd love further details.
The one I've been working with is due to Pierre L'Ecuyer [1] and is known as MRG32k3a. It's a combined multiple recursive linear congruential generator with 6 words of state. The formulas are r1[i] = (a12 * r1[i-2] + a13 * r1[i-3]) % m1 r2[i] = (a21 * r2[i-1] + a23 * r2[i-3]) % m2 r[i] = (r1[i] - r2[i]) * m1 where m1 = 2**32 - 209 m2 = 2**32 - 22835 a12 = 1403580 a13 = -810728 a21 = 527612 a23 = -1370589 If you consider the state to be made up of two 3-word state vectors, then there are two 3x3 matrices which map a given state onto the next state. So to jump ahead n steps in the sequence, you raise these matrices to the power of n. I've attached some code implementing this generator together with the jumping-ahead. (Sorry it's in C++, I hadn't discovered Python when I wrote it.) [1] Pierre L'Ecuyer, Good Parameters and Implementations for Combined Multiple Recursive Random Number Generators, Operations Research v47 no1 Jan-Feb 1999 http://www.iro.umontreal.ca/~lecuyer/myftp/papers/combmrg2.ps -- Greg /* * cmr_random_generator.C * ====================== * * Combined Multiple Recursive random number generator. * * This is an implementation of Pierre L'Ecuyer's * MRG32k3a generator, described in: * * Pierre L'Ecuyer, Good Parameters and Implementations for * Combined Multiple Recursive Random Number Generators, * Operations Research v47 no1 Jan-Feb 1999 */ #include "cmr_random_generator.H" static const double norm = 2.328306549295728e-10, m1 = 4294967087.0, m2 = 4294944443.0, a12 = 1403580.0, a13 = -810728.0, a21 = 527612.0, a23 = -1370589.0; static double a[2][3][3] = { { {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {a13, a12, 0.0} }, { {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {a23, 0.0, a21} } }; static double m[2] = { m1, m2 }; static double init_s[2][3] = { {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0} }; inline static double mod(double x, double m) { long k = (long)(x / m); x -= k * m; if (x < 0.0) x += m; return x; } /* * Initialisation */ CMRRandomGenerator::CMRRandomGenerator() { for (int i = 0; i <= 1; i++) for (int j = 0; j <= 2; j++) s[i][j] = init_s[i][j]; } /* * Advance CMRG one step and return next number */ double CMRRandomGenerator::Next() { double p1 = mod(a12 * s[0][1] + a13 * s[0][0], m1); s[0][0] = s[0][1]; s[0][1] = s[0][2]; s[0][2] = p1; double p2 = mod(a21 * s[1][2] + a23 * s[1][0], m2); s[1][0] = s[1][1]; s[1][1] = s[1][2]; s[1][2] = p2; double p = p1 - p2; if (p < 0.0) p += m1; return (p + 1) * norm; } typedef unsigned long long Int64; typedef Int64 CMRG_Vector[3]; typedef Int64 CMRG_Matrix[3][3]; static Int64 ftoi(double x, double m) { if (x >= 0.0) return Int64(x); else return Int64((long double)x + (long double)m); } static double itof(Int64 i, Int64 m) { return i; } static void v_ftoi(double u[], CMRG_Vector v, double m) { for (int i = 0; i <= 2; i++) v[i] = ftoi(u[i], m); } static void v_itof(CMRG_Vector u, double v[], Int64 m) { for (int i = 0; i <= 2; i++) v[i] = itof(u[i], m); } static void v_copy(CMRG_Vector u, CMRG_Vector v) { for (int i = 0; i <= 2; i++) v[i] = u[i]; } static void m_ftoi(double a[][3], CMRG_Matrix b, double m) { for (int i = 0; i <= 2; i++) for (int j = 0; j <= 2; j++) b[i][j] = ftoi(a[i][j], m); } static void m_copy(CMRG_Matrix a, CMRG_Matrix b) { for (int i = 0; i <= 2; i++) for (int j = 0; j <= 2; j++) b[i][j] = a[i][j]; } static void mv_mul(CMRG_Matrix a, CMRG_Vector u, CMRG_Vector v, Int64 m) { CMRG_Vector w; int i, j; for (i = 0; i <= 2; i++) { w[i] = 0; for (j = 0; j <= 2; j++) w[i] = (a[i][j] * u[j] + w[i]) % m; } v_copy(w, v); } static void mm_mul(CMRG_Matrix a, CMRG_Matrix b, CMRG_Matrix c, Int64 m) { CMRG_Matrix d; int i, j, k; for (i = 0; i <= 2; i++) { for (j = 0; j <= 2; j++) { d[i][j] = 0; for (k = 0; k <= 2; k++) d[i][j] = (a[i][k] * b[k][j] + d[i][j]) % m; } } m_copy(d, c); } /* * Advance the CMRG by n*2^e steps */ void CMRRandomGenerator::Advance(unsigned long n, unsigned int e) { CMRG_Matrix B[2]; CMRG_Vector S[2]; Int64 M[2]; int i; for (i = 0; i <= 1; i++) { m_ftoi(a[i], B[i], m[i]); v_ftoi(s[i], S[i], m[i]); M[i] = Int64(m[i]); } while (e--) { for (i = 0; i <= 1; i++) mm_mul(B[i], B[i], B[i], M[i]); } while (n) { if (n & 1) for (i = 0; i <= 1; i++) mv_mul(B[i], S[i], S[i], M[i]); n >>= 1; if (n) for (i = 0; i <= 1; i++) mm_mul(B[i], B[i], B[i], M[i]); } for (i = 0; i <= 1; i++) v_itof(S[i], s[i], M[i]); } /* * cmr_random_generator.H * ====================== * * Combined Multiple Recursive random number generator. */ #ifndef cmr_random_generator_H #define cmr_random_generator_H class CMRRandomGenerator { public: CMRRandomGenerator(); double Next(); void Advance(unsigned long mant, unsigned int exp); protected: double s[2][3]; }; #endif