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
Greg Ewing