Mersenne Twister Source Code Library

Article

Organization

Home » Linoleum Source Code » Mersenne Twister » Source Code Library

Scope

This is the source code of the Mersenne Twister library, implemented in the Linoleum programming language.

Author

Herbert Glarner

Source Code

(Random Number Generator "RNG.txt" - Copyright 2006-2007 Herbert Glarner)

(Algorithm: MT19937 algorithm copyright by Takuji Nishimura and Makoto 
Matsumoto, 2002)

(*****************************************************************************)

    (=========================================================================)
    (    Copyright by Herbert Glarner [2006, 2007]
        Mail: herbert.glarner@bluewin.ch
         Free usage granted for all purposes, provided above credits are given. 
         No liability for any damage.    )
    (=========================================================================)

    (=========================================================================)
    (Routines summary:
        Name of routine                    Short description)
    (-------------------------------------------------------------------------)
    (    PUBLIC ROUTINES:
            Rng Init                    Array initialization
            Rng Long                    Returns a Random Integer 0...2^32-1

        INTERNAL ROUTINES:
            Rng Init Genrand            Initialization
            
        Of course one can use also Internal Routines from the outside world,
        but read the assumptions in the comment header of those routines if
        you do.    )
    (=========================================================================)

(*****************************************************************************)
"libraries"

(*****************************************************************************)
"directors"
        program name = { RNG_library };
        unit = 32;
         
(*****************************************************************************)
"constants"                                 (Convention: "RNG ...", UPPER CASE)

    (=========================================================================)
    (Error numbers for variable "rng error")
    (-------------------------------------------------------------------------)
    (Errors codes)
        RNG ERR NONE                = 0;    (No error)
    (=========================================================================)

    (=========================================================================)
    (Period parameters)
    (-------------------------------------------------------------------------)
        RNG N                       = 624;
        RNG M                       = 397;
        RNG MATRIX A                = 9908b0dfh;        (Constant vector A)
        RNG UPPER MASK              = 80000000h;
        RNG LOWER MASK              = 7fffffffh;
    (=========================================================================)
        
(*****************************************************************************)
"variables"                                 (Convention: "rng ...", lower case)

    (=========================================================================)
    (Static variable for 'Rng Long')
    (-------------------------------------------------------------------------)
        rng mag                     = 0;
                                      RNG MATRIX A;
    (=========================================================================)

    (=========================================================================)
    (Initialization)
    (-------------------------------------------------------------------------)
        rng mti                     = RNG N plus 1;     (init with 625)
    (=========================================================================)

    (=========================================================================)
    (Takes one of the RAM ERR constants)
    (-------------------------------------------------------------------------)
        rng error                   = 0;
    (=========================================================================)

(*****************************************************************************)
"workspace"                                 (Convention: "rng ...", lower case)

    (=========================================================================)
    (Array for the state vector)
    (-------------------------------------------------------------------------)
        rng mt                      = RNG N;            (624 elements)
    (=========================================================================)

(*****************************************************************************)
"programme"                                 (Convention: "Rng ...", Mixed case)

    (=========================================================================)
    (This is a library only: It can not be called as a stand-alone 
    application.)
    (-------------------------------------------------------------------------)
        end;
    (=========================================================================)

(*****************************************************************************)
"Rng Long"

    (=========================================================================)
    (Generates a random integer 0...ffffffffh)
    (-------------------------------------------------------------------------)
    (Input:        -
    Output:        E: Random integer 0...ffffffffh
    
    All registers unchanged.
    )
    (-------------------------------------------------------------------------)
        A -->; B -->; C -->; D -->; 
    (-------------------------------------------------------------------------)
        (unsigned long y;)
        (y is register E)
    (-------------------------------------------------------------------------)
        (static unsigned long mag01[2]={0x0UL, MATRIX_A};)
        (addressed with [rng mag] and [rng mag plus 1])
        (initialized when declared)
    (-------------------------------------------------------------------------)
        (if {mti >= N} { /* generate N words at one time */)
        ? [rng mti] < RNG N -> Rng Long No New Words;
    (-------------------------------------------------------------------------)
        (int kk;)
        C = 0;
    (-------------------------------------------------------------------------)
        ( if {mti == N+1}   /* if init_genrand{} has not been called, */)
        ? [rng mti] != RNG N plus 1 -> Rng Long Init Ok;
        
        (init_genrand{5489UL}; /* a default initial seed is used */)
        E = 5489;
        => Rng Init Genrand;
        
    "Rng Long Init Ok"
    (-------------------------------------------------------------------------)
        (First loop)
        
        (for {kk=0;kk<N-M;kk++} )
        ? C >= RNG N minus RNG M -> Rng Long First Loop Ok;
    "Rng Long First Loop"
    (-------------------------------------------------------------------------)
        (y = {mt[kk]&UPPER_MASK}|{mt[kk+1]&LOWER_MASK};)
        E = [C plus rng mt];

        (remains: y = {E&UPPER_MASK}|{mt[kk+1]&LOWER_MASK};)
        E & RNG UPPER MASK;

        (remains: y = {E}|{mt[kk+1]&LOWER_MASK};)
        B = [C plus rng mt plus 1];

        (remains: y = {E}|{B&LOWER_MASK};)
        B & RNG LOWER MASK;

        (remains: E = E|B;)
        E | B;
    (-------------------------------------------------------------------------)
        (mt[kk] = mt[kk+M] ^ {y >> 1} ^ mag01[y & 0x1UL];)
        B = E;
        B & 1h;
        B = [B plus rng mag];               (mag01[y & 0x1UL])

        (remains: mt[kk] = mt[kk+M] ^ {E >> 1} ^ B;)
        E > 1;

        (remains: mt[kk] = mt[kk+M] ^ E ^ B;)
        (sequence irrelevant: xor is associative)
        E # B;

        (remains: mt[kk] = mt[kk+M] ^ E;)
        B = [C plus RNG M plus rng mt];

        (remains: mt[kk] = B ^ E;)
        E # B;

        (remains: mt[kk] = E;)
        [C plus rng mt] = E;
    (-------------------------------------------------------------------------)
        (Next element)
        C +;                                                (kk++)
        ? C < RNG N minus RNG M -> Rng Long First Loop;     (while kk<N-M)

    "Rng Long First Loop Ok"
    (-------------------------------------------------------------------------)
        (Second Loop)
        
        (for {;kk<N-1;kk++} )
        ? C >= RNG N minus 1 -> Rng Long Second Loop Ok;
    "Rng Long Second Loop"
    (-------------------------------------------------------------------------)
        (y = {mt[kk]&UPPER_MASK}|{mt[kk+1]&LOWER_MASK};)
        E = [C plus rng mt];

        (remains: y = {E&UPPER_MASK}|{mt[kk+1]&LOWER_MASK};)
        E & RNG UPPER MASK;

        (remains: y = {E}|{mt[kk+1]&LOWER_MASK};)
        B = [C plus rng mt plus 1];

        (remains: y = {E}|{B&LOWER_MASK};)
        B & RNG LOWER MASK;

        (remains: y = E|B;)
        E | B;
    (-------------------------------------------------------------------------)
        (mt[kk] = mt[kk+{M-N}] ^ {y >> 1} ^ mag01[y & 0x1UL];)
        B = E;
        B & 1h;
        B = [B plus rng mag];               (mag01[y & 0x1UL])

        (remains: mt[kk] = mt[kk+{M-N}] ^ {E >> 1} ^ B;)
        E > 1;

        (remains: mt[kk] = mt[kk+{M-N}] ^ E ^ B;)
        (sequence irrelevant: xor is associative)
        E # B;

        (remains: mt[kk] = mt[kk+{M-N}] ^ E;)
        B = [C plus RNG M minus RNG N plus rng mt];

        (remains: mt[kk] = B ^ E;)
        E # B;

        (remains: mt[kk] = E;)
        [C plus rng mt] = E;
    (-------------------------------------------------------------------------)
        C +;
        ? C < RNG N minus 1 -> Rng Long Second Loop;    (while kk<N-1)

    "Rng Long Second Loop Ok"
    (-------------------------------------------------------------------------)
        (y = {mt[N-1]&UPPER_MASK}|{mt[0]&LOWER_MASK};)
        E = [rng mt];

        (remains: y = {mt[N-1]&UPPER_MASK}|{E&LOWER_MASK};)
        E & RNG LOWER MASK;

        (remains: y = {mt[N-1]&UPPER_MASK}|{E};)
        B = [RNG N minus 1 plus rng mt];

        (remains: y = {B&UPPER_MASK}|{E};)
        B & RNG UPPER MASK;

        (remains: y = B|E;)
        E | B;
    (-------------------------------------------------------------------------)
        (mt[N-1] = mt[M-1] ^ {y >> 1} ^ mag01[y & 0x1UL];)
        B = E;
        B & 1h;

        (remains: mt[N-1] = mt[M-1] ^ {E >> 1} ^ mag01[B];)
        B = [B plus rng mag];

        (remains: mt[N-1] = mt[M-1] ^ {E >> 1} ^ B;)
        E > 1;

        (remains: mt[N-1] = mt[M-1] ^ E ^ B;)
        E # B;

        (remains: mt[N-1] = mt[M-1] ^ E;)
        B = [RNG M minus 1 plus rng mt];

        (remains: mt[N-1] = B ^ E;)
        E # B;

        (remains: mt[N-1] = E;)
        [rng mt plus RNG N minus 1] = E;

    (-------------------------------------------------------------------------)
        (mti = 0;)
        [rng mti] = 0;
    (-------------------------------------------------------------------------)
    "Rng Long No New Words"
        (y = mt[mti++];)
        E = [rng mti];
        E = [E plus rng mt];
        [rng mti] +;
    (-------------------------------------------------------------------------)
        (Tempering)
        
        (y ^= {y >> 11};)
        B = E;
        B > 11;
        E # B;
    (-------------------------------------------------------------------------)
        (y ^= {y << 7} & 0x9d2c5680UL;)
        B = E;
        B < 7;
        B & 9d2c5680h;
        E # B;
    (-------------------------------------------------------------------------)
        (y ^= {y << 15} & 0xefc60000UL;)
        B = E;
        B < 15;
        B & efc60000h;
        E # B;
    (-------------------------------------------------------------------------)
        (y ^= {y >> 18};)
        B = E;
        B > 18;
        E # B;
    (-------------------------------------------------------------------------)
        <-- D; <-- C; <-- B; <-- A;
        [rng error] = RNG ERR NONE;
        end;
    (=========================================================================)

(*****************************************************************************)
"Rng Init"

    (=========================================================================)
    (Initializize by array with array length)
    (-------------------------------------------------------------------------)
    (Input:        E:    addr of init key[]  {32 bits each}
                   D:  length of key        {num of elements}
    Output:        -
    
    All registers unchanged.
    )
    (-------------------------------------------------------------------------)
        A -->; B -->; C -->; D -->; E -->; 
    (-------------------------------------------------------------------------)
        (init_genrand{19650218UL};)
        (This number is just the birthday of one of the authors, chosen
        without a reason ;)
        E = 19650218;
        => Rng Init Genrand;
    (-------------------------------------------------------------------------)
        <-- E;                              (restore E)
        E -->; 
    (-------------------------------------------------------------------------)
        (i=1; j=0;)
        A = 1; B = 0;                       (A: i;  B: j)
    (-------------------------------------------------------------------------)
        (k = {N>key_length ? N : key_length};)
        (equal to: k = {key_length<=N ? N : key_length};)
        C = RNG N;
        ? D <= RNG N -> Rng Init Key Length;
        C = D;
    "Rng Init Key Length"
    (-------------------------------------------------------------------------)
        (First loop)
        
        (for {; k; k--} )                   (k is in Reg. C)
        ? C = 0 -> Rng Init First Loop Ok;
    "Rng Init First Loop"
    (-------------------------------------------------------------------------)
        (mt[i] = {mt[i] ^ {{mt[i-1] ^ {mt[i-1] >> 30}} * 1664525UL}}
          + init_key[j] + j; /* non linear */ )

        C -->; A -->; B -->;                (need registers, save k, i and j)

        B = A; B -; B = [B plus rng mt];    (mt[i-1])

        (Remains: mt[i] = {mt[i] ^ {{B ^ {B >> 30}} * 1664525UL}} 
          + init_key[j] + j;)
        C = B;                                
        C > 30;                             (B >> 30)

        (Remains: mt[i] = {mt[i] ^ {{B ^ C} * 1664525UL}} + init_key[j] + j;)
        B # C;                              (B ^ C)

        (Remains: mt[i] = {mt[i] ^ {B * 1664525UL}} + init_key[j] + j;)
        B * 1664525;

        (Remains: mt[i] = {mt[i] ^ B} + init_key[j] + j;)
        C = [A plus rng mt];                (mt[i])

        (Remains: mt[i] = {C ^ B} + init_key[j] + j;)
        C # B;                              (C ^ B)

        (Remains: mt[i] = C + init_key[j] + j;)
        <-- B;                              (restore j)
        B -->;                              (Still need later)

        (Remains: mt[i] = C + init_key[B] + B;)
        C + B;

        (Remains: mt[i] = C + init_key[B];)
        (Address of vector init_key is in E, B is index)
        B + E;                              (Address! of element)

        (Remains: mt[i] = C + [B];)
        C + [B];

        (mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */)
        C & ffffffffh;

        (Remains: mt[i] = C;)
        <-- B;                              (restore j)
        <-- A;                              (restore i)
        
        (Remains: mt[A] = C;)
        [A plus rng mt] = C;                (mt[A] = C;)

        (Restoring of k in Reg. C delayed)        
    (-------------------------------------------------------------------------)
        (i++; j++;)
        A +; B +;
    (-------------------------------------------------------------------------)
        ( if (i>=N) { mt[0] = mt[N-1]; i=1; } )
        ( equal to: if (i<N) {} else { mt[0] = mt[N-1]; i=1; } )
        
        ? A < RNG N -> Rng Init First Loop No I Corr;
        C = RNG N; C -;                     (N-1)
        C = [C plus rng mt];                (mt[N-1])
        [rng mt] = C;                       (mt[0] = mt[N-1];)
        A = 1;                              (i=1;)
    "Rng Init First Loop No I Corr"
        (Restoring of k in Reg. C delayed)        
    (-------------------------------------------------------------------------)
        (if {j>=key_length} j=0;)
        (equal to: if {j<key_length} {} j=0;)
        
        ? B < D -> Rng Init First Loop No J Corr;
        B = 0;                              (j=0;)
    "Rng Init First Loop No J Corr"
    (-------------------------------------------------------------------------)
        (Next element first loop)
        <-- C;                              (restore k)
        C ^ Rng Init First Loop;
    "Rng Init First Loop Ok"
    (-------------------------------------------------------------------------)
        (Second loop)

        (for {k=N-1; k; k--})
        C = RNG N; C -;
        ? C = 0 -> Rng Init Second Loop Ok;
    "Rng Init Second Loop"
    (-------------------------------------------------------------------------)
        C -->; B -->;                       (need registers, save k and j)

        (mt[i] = {mt[i] ^ {{mt[i-1] ^ {mt[i-1] >> 30}} * 1566083941UL}} - i;)
        B = A; B -; B = [B plus rng mt];    (mt[i-1])

        (Remains: mt[i] = {mt[i] ^ {{B ^ {B >> 30}} * 1566083941UL}} - i;)
        C = B;
        C > 30;                             (B >> 30)

        (Remains: mt[i] = {mt[i] ^ {{B ^ C} * 1566083941UL}} - i;)
        B # C;                              (B ^ C)

        (Remains: mt[i] = {mt[i] ^ {B * 1566083941UL}} - i;)
        B * 1566083941;                     (B * 1566083941UL)

        (Remains: mt[i] = {mt[i] ^ B} - i;)
        C = [A plus rng mt];                (mt[i])

        (Remains: mt[i] = {C ^ B} - i;)
        C # B;                              (C ^ B)

        (Remains: mt[i] = C - i;)
        C - A;

        (mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */)
        C & ffffffffh;

        (Remains: mt[i] = C;)
        [A plus rng mt] = C;                (mt[i] = C)

        <-- B;                              (restore j)
        (Restauration of k into C delayed)
    (-------------------------------------------------------------------------)
        (i++;)
        A +;
        (Restauration of k into C delayed)
    (-------------------------------------------------------------------------)
        (if {i>=N} { mt[0] = mt[N-1]; i=1; })
        (equals to: if {i<N} {} { mt[0] = mt[N-1]; i=1; })
        
        ? A < RNG N -> Rng Init Second Loop No I Corr;
        C = RNG N; C -;
        C = [C plus rng mt];                (mt[N-1];)
        [rng mt] = C;                       (mt[0] = mt[N-1];)
        A = 1;                              (i=1;)
    "Rng Init Second Loop No I Corr"
    (-------------------------------------------------------------------------)
        (Next element second loop)
        <-- C;                              (restore k)
        C ^ Rng Init Second Loop;
    "Rng Init Second Loop Ok"
    (-------------------------------------------------------------------------)
        (mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */)
        [rng mt] = 80000000h;
    (-------------------------------------------------------------------------)
    "Rng Init Exit"
        <-- E; <-- D; <-- C; <-- B; <-- A;
        [rng error] = RNG ERR NONE;
        end;
    (=========================================================================)
    
(*****************************************************************************)
"Rng Init Genrand"

    (=========================================================================)
    (Initializizes rng mt[N] with a seed)
    (-------------------------------------------------------------------------)
    (Input:        E:    Seed 32 bit
    Output:        -
    
    All registers unchanged.
    )
    (-------------------------------------------------------------------------)
        A -->; B -->; C -->; D -->; E -->; 
    (-------------------------------------------------------------------------)
        (mt[0]= s & 0xffffffffUL;)
        E & ffffffffh;
        [rng mt] = E;
    (-------------------------------------------------------------------------)
        (for {mti=1; mti<N; mti++} )
        [rng mti] = 1;                      (mti=1;)
    "Rng Init Genrand Loop"
    (-------------------------------------------------------------------------)
        (mt[mti] = {1812433253UL * {mt[mti-1] ^ {mt[mti-1] >> 30}} + mti};)
        B = [rng mti]; B -; B = [B plus rng mt];    (mt[mti-1])
        
        (Remains: mt[mti] = {1812433253UL * {B ^ {B >> 30}} + mti};)
        C = B; C > 30;                      (B >> 30)

        (Remains: mt[mti] = {1812433253UL * {B ^ C} + mti};)
        B # C;                              (B ^ C)

        (Remains: mt[mti] = {1812433253UL * {B} + mti};)
        (This constant is a multiplier for a linear congruential generator,
        here chosen without reason)
        B * 1812433253;                     (1812433253UL * {B} )

        (Remains: mt[mti] = {B + A};)
        B + [rng mti];                      (B + mti)

        (mt[mti] &= 0xffffffffUL; for >32 bit machines)
        B & ffffffffh;                      (mt[mti] &= 0xffffffffUL)

        (Remains: mt[mti] = {B};)
        A = [rng mti];
        [A plus rng mt] = B;                (mt[mti] = {B} )
    (-------------------------------------------------------------------------)
        (Next element)
        [rng mti] +;
        ? [rng mti] < RNG N -> Rng Init Genrand Loop;
    (-------------------------------------------------------------------------)
        <-- E; <-- D; <-- C; <-- B; <-- A;
        [rng error] = RNG ERR NONE;
        end;
    (=========================================================================)

(*****************************************************************************)

Linoleum Syntax Highlighting produced with LSH (© 2007 by Herbert Glarner)