Mersenne Twister implentation

[KeithB]

Based on the discussion in the random number thread, here is an implementation of the Mersenne Twister.

To use, call Twister which will set up the state in the list L9(). I used a list instead of a matrix since a list can store integers. After the initial Setup you can call:

MTSEED(number) to reset the seed.
RANDMT() to get a uniform number from 0 to 1
RANDMTI() to get the raw hex numbers
RANDMTN() to get normally distributed numbers with a mean of 0 and SD of 1.

Observations:
The output matches the reference implementations exactly.
Bit twiddling is not the Prime's strong suit. (Dogs walking on their hind legs come to mind)
It is much slower than the built in functions. My random pi program is slowed down by a factor of about 20, 21.5 sec for 10000 iterations vs 1.2 sec for 10000 iterations using the built-in rand.

Code:


//
// Twister
//An Implementation of the Mersenne
//Twister Random Number Algorithm by 
//Makoto Matsumoto and Takuji Nishimura
//Normalizer from Abramowitz and Stegun
//
// Basic flow cribbed from Stephan Brumme
//
// Keith Barkley, Feb 2017
//
//DECLARES
MTSEED();
TWIST();
RANDMTI();
RANDMT();
RANDMTN();

EXPORT Twister()
BEGIN
//THE STATE SIZE IS 624
// 625 FOR INDEX
// 626 FOR NORMAL "PHASE"
//INIT
LOCAL i;
FOR i FROM 1 TO 624 DO
L9(i):=#0d;
END;
L9(625):=625;
L9(626):=0;

MTSEED(1);
END;

EXPORT MTSEED(S)
BEGIN
LOCAL i;
L9(1):=R→B(S);
FOR i FROM 2 TO 624 DO
L9(i):=(#1812433253d *BITXOR(L9(i-1),BITSR(L9(i-1),30)))+R→B(i-1);
END;
TWIST();
L9(625):=1;
END;

TWIST()
BEGIN
LOCAL b,b1,b2,i;
//TWIST
FOR i FROM 1 TO 227 DO
b:=BITOR(BITAND(L9(i),#80000000h),BITAND(L9(i+1),#7FFFFFFFh));
b1:=BITSR(R→B(b),1);
b2:=BITAND(R→B(b),#1h) * #9908B0DFh;
L9(i):=BITXOR(L9(i+397),R→B(b1),R→B(b2));
END;

FOR i FROM 228 TO 623 DO
b:=BITOR(BITAND(L9(i),#80000000h),BITAND(L9(i+1),#7FFFFFFFh));
b1:=BITSR(R→B(b),1);
b2:=BITAND(R→B(b),#1h) * #9908B0DFh;
L9(i):=BITXOR(L9(i-227),R→B(b1),R→B(b2));
END;

b:=BITOR(BITAND(L9(624),#80000000h),BITAND(L9(1),#7FFFFFFFh));
b1:=BITSR(R→B(b),1);
b2:=BITAND(R→B(b),#1h) * #9908B0DFh;
L9(624):=BITXOR(L9(397),R→B(b1),R→B(b2));

END;

EXPORT RANDMTI()
BEGIN
LOCAL x;

IF (L9(625)=625) THEN 
TWIST();
L9(625):=1;
END;
x:=L9(L9(625));

x:= BITXOR(R→B(x),BITSR(R→B(x),11));
x:= BITXOR(R→B(x),BITAND(BITSL(R→B(x),7),#9D2C5680h));
x:= BITXOR(R→B(x),BITAND(BITSL(R→B(x),15),#EFC60000h));
x:= BITXOR(R→B(x),BITSR(R→B(x),18));

L9(625):=L9(625)+1;
RETURN x;
END;

EXPORT RANDMT()
BEGIN
RETURN B→R(RANDMTI())/B→R(#FFFFFFFFh);
END;
EXPORT RANDMTN()
BEGIN
LOCAL n,u1,u2;
IF (L9(626)≠0) THEN
n:=L9(626);
L9(626):=0;
RETURN n;
END;

u1:=RANDMT();
u2:=RANDMT();
L9(626):=(−2*LN(u1))^0.5*COS(2*π*u2);
RETURN ((−2*LN(u1))^0.5*SIN(2*π*u2));
END;