Pseudo-random numbers/Combined recursive generator MRG32k3a: Difference between revisions

Content added Content deleted

Inline

Revision as of 11:29, 28 August 2022

MRG32k3a Combined recursive generator (pseudo-code)

   /* Constants */
   /* First generator */
   a1 = [0, 1403580, -810728]
   m1 = 2**32 - 209
   /* Second Generator */
   a2 = [527612, 0, -1370589]
   m2 = 2**32 - 22853
    
   d = m1 + 1
   
   class MRG32k3a
       x1 = [0, 0, 0]  /* list of three last values of gen #1 */
       x2 = [0, 0, 0]  /* list of three last values of gen #2 */
       
       method seed(u64 seed_state)
           assert seed_state in range >0 and < d 
           x1 = [seed_state, 0, 0]
           x2 = [seed_state, 0, 0]
       end method
           
       method next_int()
           x1i = (a1[0]*x1[0] + a1[1]*x1[1] + a1[2]*x1[2]) mod m1
           x2i = (a2[0]*x2[0] + a2[1]*x2[1] + a2[2]*x2[2]) mod m2
           x1 = [x1i, x1[0], x1[1]]    /* Keep last three */
           x2 = [x2i, x2[0], x2[1]]    /* Keep last three */
           z = (x1i - x2i) % m1
           answer = (z + 1)
           
           return answer
       end method
       
       method next_float():
           return float next_int() / d
       end method
       
   end class

MRG32k3a Use:

   random_gen = instance MRG32k3a
   random_gen.seed(1234567)
   print(random_gen.next_int())   /* 1459213977 */
   print(random_gen.next_int())   /* 2827710106 */
   print(random_gen.next_int())   /* 4245671317 */
   print(random_gen.next_int())   /* 3877608661 */
   print(random_gen.next_int())   /* 2595287583 */

Task

Generate a class/set of functions that generates pseudo-random

numbers as shown above.

Show that the first five integers generated with the seed `1234567`

are as shown above

Show that for an initial seed of '987654321' the counts of 100_000

repetitions of

   floor(random_gen.next_float() * 5)

Is as follows:

   0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931

Show your output here, on this page.

11l

Translation of: Python

V a1 = [Int64(0), 1403580, -810728]
V m1 = Int64(2) ^ 32 - 209
V a2 = [Int64(527612), 0, -1370589]
V m2 = Int64(2) ^ 32 - 22853
V d = m1 + 1

T MRG32k3a
   [Int64] x1, x2

   F (seed_state = 123)
      .seed(seed_state)

   F seed(Int64 seed_state)
      assert(seed_state C Int64(0) <.< :d, ‘Out of Range 0 x < #.’.format(:d))
      .x1 = [Int64(seed_state), 0, 0]
      .x2 = [Int64(seed_state), 0, 0]

   F next_int()
      ‘return random int in range 0..d’
      V x1i = (sum(zip(:a1, .x1).map((aa, xx) -> aa * xx)) % :m1 + :m1) % :m1
      V x2i = (sum(zip(:a2, .x2).map((aa, xx) -> aa * xx)) % :m2 + :m2) % :m2
      .x1 = [x1i] [+] .x1[0.<2]
      .x2 = [x2i] [+] .x2[0.<2]
      V z = ((x1i - x2i) % :m1 + :m1) % :m1
      R z + 1

   F next_float()
      ‘return random float between 0 and 1’
      R Float(.next_int()) / :d

V random_gen = MRG32k3a()
random_gen.seed(1234567)
L 5
   print(random_gen.next_int())

random_gen.seed(987654321)
V hist = Dict(0.<5, i -> (i, 0))
L 100'000
   hist[Int(random_gen.next_float() * 5)]++
print(hist)

Output:

1459213977
2827710106
4245671317
3877608661
2595287583
[0 = 20002, 1 = 20060, 2 = 19948, 3 = 20059, 4 = 19931]

Ada

package MRG32KA is
   type I64 is range -2**63..2**63 - 1;
   m1 : constant I64 := 2**32 - 209;
   m2 : constant I64 := 2**32 - 22853;

   subtype state_value is I64 range 1..m1;
   
   procedure Seed (seed_state : state_value);
   function Next_Int return I64;
   function Next_Float return Long_Float;
end MRG32KA;

package body MRG32KA is
   
   type Data_Array is array (0..2) of I64;
   
   d : constant I64 := m1 + 1;
   ----------------
   -- Generators --
   ----------------
   
   a1 : Data_Array := (0, 1403580, -810728);
   a2 : Data_Array := (527612, 0, -1370589);
   
   x1 : Data_Array := (0, 0, 0);
   x2 : Data_Array := (0, 0, 0);
   ----------
   -- Seed --
   ----------

   procedure Seed (seed_state : state_value) is
   begin
      x1 := (seed_state, 0, 0);
      x2 := (seed_state, 0, 0);
   end Seed;

   --------------
   -- Next_Int --
   --------------

   function Next_Int return I64 is
      x1i : i64;
      x2i : I64;
      z   : I64;
      answer : I64;
   begin
      x1i := (a1(0) * x1(0) + a1(1) * x1(1) + a1(2) * x1(2)) mod m1;
      x2i := (a2(0) * x2(0) + a2(1) * x2(1) + a2(2) * x2(2)) mod m2;
      x1  := (x1i, x1(0), x1(1));
      x2  := (x2i, x2(0), x2(1));
      z := (x1i - x2i) mod m1;
      answer := z + 1;
      return answer;
   end Next_Int;

   ----------------
   -- Next_Float --
   ----------------

   function Next_Float return Long_Float is
   begin
      return Long_float(Next_Int) / Long_Float(d);
   end Next_Float;

end MRG32KA;

with Ada.Text_IO; use Ada.Text_IO;
with mrg32ka; use mrg32ka;

procedure Main is
   counts : array(0..4) of Natural := (Others => 0);
   J : Natural;
begin
   
   seed(1234567);
   for I in 1..5 loop
      Put_Line(I64'Image(Next_Int));
   end loop;
   New_Line;
   seed(987654321);
   
   for I in 1..100_000 loop
      J := Natural(Long_Float'Floor(Next_Float * 5.0));
      Counts(J) := Counts(J) + 1;
   end loop;
   
   for I in Counts'Range loop
      Put(I'Image & " :" & Counts(I)'Image);
   end loop;
   
end Main;

Output:

 1459213977
 2827710106
 4245671317
 3877608661
 2595287583

 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931

C

#include <math.h>
#include <stdio.h>
#include <stdint.h>

int64_t mod(int64_t x, int64_t y) {
    int64_t m = x % y;
    if (m < 0) {
        if (y < 0) {
            return m - y;
        } else {
            return m + y;
        }
    }
    return m;
}

// Constants
// First generator
const static int64_t a1[3] = { 0, 1403580, -810728 };
const static int64_t m1 = (1LL << 32) - 209;
// Second generator
const static int64_t a2[3] = { 527612, 0, -1370589 };
const static int64_t m2 = (1LL << 32) - 22853;

const static int64_t d = (1LL << 32) - 209 + 1; // m1 + 1

// the last three values of the first generator
static int64_t x1[3];
// the last three values of the second generator
static int64_t x2[3];

void seed(int64_t seed_state) {
    x1[0] = seed_state;
    x1[1] = 0;
    x1[2] = 0;

    x2[0] = seed_state;
    x2[1] = 0;
    x2[2] = 0;
}

int64_t next_int() {
    int64_t x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
    int64_t x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
    int64_t z = mod(x1i - x2i, m1);

    // keep last three values of the first generator
    x1[2] = x1[1];
    x1[1] = x1[0];
    x1[0] = x1i;

    // keep last three values of the second generator
    x2[2] = x2[1];
    x2[1] = x2[0];
    x2[0] = x2i;

    return z + 1;
}

double next_float() {
    return (double)next_int() / d;
}

int main() {
    int counts[5] = { 0, 0, 0, 0, 0 };
    int i;

    seed(1234567);
    printf("%lld\n", next_int());
    printf("%lld\n", next_int());
    printf("%lld\n", next_int());
    printf("%lld\n", next_int());
    printf("%lld\n", next_int());
    printf("\n");

    seed(987654321);
    for (i = 0; i < 100000; i++) {
        int64_t value = floor(next_float() * 5);
        counts[value]++;
    }
    for (i = 0; i < 5; i++) {
        printf("%d: %d\n", i, counts[i]);
    }

    return 0;
}

Output:

C++

Translation of: C

#include <array>
#include <iostream>

int64_t mod(int64_t x, int64_t y) {
    int64_t m = x % y;
    if (m < 0) {
        if (y < 0) {
            return m - y;
        } else {
            return m + y;
        }
    }
    return m;
}

class RNG {
private:
    // First generator
    const std::array<int64_t, 3> a1{ 0, 1403580, -810728 };
    const int64_t m1 = (1LL << 32) - 209;
    std::array<int64_t, 3> x1;
    // Second generator
    const std::array<int64_t, 3> a2{ 527612, 0, -1370589 };
    const int64_t m2 = (1LL << 32) - 22853;
    std::array<int64_t, 3> x2;
    // other
    const int64_t d = (1LL << 32) - 209 + 1; // m1 + 1

public:
    void seed(int64_t state) {
        x1 = { state, 0, 0 };
        x2 = { state, 0, 0 };
    }

    int64_t next_int() {
        int64_t x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
        int64_t x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
        int64_t z = mod(x1i - x2i, m1);

        // keep last three values of the first generator
        x1 = { x1i, x1[0], x1[1] };
        // keep last three values of the second generator
        x2 = { x2i, x2[0], x2[1] };

        return z + 1;
    }

    double next_float() {
        return static_cast<double>(next_int()) / d;
    }
};

int main() {
    RNG rng;

    rng.seed(1234567);
    std::cout << rng.next_int() << '\n';
    std::cout << rng.next_int() << '\n';
    std::cout << rng.next_int() << '\n';
    std::cout << rng.next_int() << '\n';
    std::cout << rng.next_int() << '\n';
    std::cout << '\n';

    std::array<int, 5> counts{ 0, 0, 0, 0, 0 };
    rng.seed(987654321);
    for (size_t i = 0; i < 100000; i++) 		{
        auto value = floor(rng.next_float() * 5.0);
        counts[value]++;
    }
    for (size_t i = 0; i < counts.size(); i++) 		{
        std::cout << i << ": " << counts[i] << '\n';
    }

    return 0;
}

Output:

D

Translation of: C++

import std.math;
import std.stdio;

long mod(long x, long y) {
    long m = x % y;
    if (m < 0) {
        if (y < 0) {
            return m - y;
        } else {
            return m + y;
        }
    }
    return m;
}

class RNG {
private:
    // First generator
    immutable(long []) a1 = [0, 1403580, -810728];
    immutable long m1 = (1L << 32) - 209;
    long[3] x1;
    // Second generator
    immutable(long []) a2 = [527612, 0, -1370589];
    immutable long m2 = (1L << 32) - 22853;
    long[3] x2;
    // other
    immutable long d = m1 + 1;

public:
    void seed(long state) {
        x1 = [state, 0, 0];
        x2 = [state, 0, 0];
    }

    long next_int() {
        long x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
        long x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
        long z = mod(x1i - x2i, m1);

        // keep the last three values of the first generator
        x1 = [x1i, x1[0], x1[1]];
        // keep the last three values of the second generator
        x2 = [x2i, x2[0], x2[1]];

        return z + 1;
    }

    double next_float() {
        return cast(double) next_int() / d;
    }
}

void main() {
    auto rng = new RNG();

    rng.seed(1234567);
    writeln(rng.next_int);
    writeln(rng.next_int);
    writeln(rng.next_int);
    writeln(rng.next_int);
    writeln(rng.next_int);
    writeln;

    int[5] counts;
    rng.seed(987654321);
    foreach (i; 0 .. 100_000) {
        auto value = cast(int) floor(rng.next_float * 5.0);
        counts[value]++;
    }
    foreach (i,v; counts) {
        writeln(i, ": ", v);
    }
}

Output:

Factor

USING: arrays kernel math math.order math.statistics
math.vectors prettyprint sequences ;

CONSTANT: m1 4294967087
CONSTANT: m2 4294944443

: seed ( n -- seq1 seq2 )
    dup 1 m1 between? t assert= 0 0 3array dup ;

: new-state ( seq1 seq2 n -- new-seq )
    [ dup ] [ vdot ] [ rem prefix but-last ] tri* ;

: next-state ( a b -- a' b' )
    [ { 0 1403580 -810728 } m1 new-state ]
    [ { 527612 0 -1370589 } m2 new-state ] bi* ;

: next-int ( a b -- a' b' n )
    next-state 2dup [ first ] bi@ - m1 rem 1 + ;

: next-float ( a b -- a' b' x ) next-int m1 1 + /f ;

! Task
1234567 seed 5 [ next-int . ] times 2drop

987654321 seed 100,000 [ next-float 5 * >integer ] replicate
2nip histogram .

Output:

1459213977
2827710106
4245671317
3877608661
2595287583
H{ { 0 20002 } { 1 20060 } { 2 19948 } { 3 20059 } { 4 19931 } }

Forth

Translation of: uBasic/4tH

Works with: 4tH v3.64

6 array (seed)                         \ holds the seed
6 array (gens)                         \ holds the generators
                                       \ set up constants
       0 (gens) 0 th !                 \ 1st generator
 1403580 (gens) 1 th !
 -810728 (gens) 2 th !
  527612 (gens) 3 th !                 \ 2nd generator
       0 (gens) 4 th !
-1370589 (gens) 5 th !

1 32 lshift   209 - value (m)          \ 1st generator constant
1 32 lshift 22853 - value (n)          \ 2nd generator constant
                                       ( n1 n2 -- n3)
: (mod) tuck mod tuck 0< if abs + ;then drop ;
: (generate) do (seed) i th @ (gens) i th @ * + loop swap (mod) ;
: (reseed) ?do (seed) i th ! loop ;    ( n1 n2 n3 limit index --)
: randomize 6 0 do dup i 3 mod if >zero then (seed) i th ! loop drop ;
                                       ( n --)
: random                               ( -- n)
  (m) 0 3 0 (generate) (n) 0 6 3 (generate) over over
  (seed) 4 th @ (seed) 3 th @ rot 6 3 (reseed)
  (seed) 1 th @ (seed) 0 th @ rot 3 0 (reseed) - (m) (mod) 1+
;

include lib/fp1.4th                    \ simple floating point support
include lib/zenfloor.4th               \ for FLOOR

5 array (count)                        \ setup an array of 5 elements

: test
  1234567 randomize
  random . cr random . cr random . cr
  random . cr random . cr cr           \ perform the first test

  987654321 randomize (m) 1+ s>f       \ set up denominator

  100000 0 ?do                         \ do this 100,000 times
    random s>f fover f/ 5 s>f f* floor f>s cells (count) + 1 swap +!
  loop fdrop
                                       \ show the results
  5 0 ?do i . ." : " (count) i th ? cr loop
;

test

Output:

Go

Translation of: Python

package main

import (
    "fmt"
    "log"
    "math"
)

var a1 = []int64{0, 1403580, -810728}
var a2 = []int64{527612, 0, -1370589}

const m1 = int64((1 << 32) - 209)
const m2 = int64((1 << 32) - 22853)
const d = m1 + 1

// Python style modulus
func mod(x, y int64) int64 {
    m := x % y
    if m < 0 {
        if y < 0 {
            return m - y
        } else {
            return m + y
        }
    }
    return m
}

type MRG32k3a struct{ x1, x2 [3]int64 }

func MRG32k3aNew() *MRG32k3a { return &MRG32k3a{} }

func (mrg *MRG32k3a) seed(seedState int64) {
    if seedState <= 0 || seedState >= d {
        log.Fatalf("Argument must be in the range [0, %d].\n", d)
    }
    mrg.x1 = [3]int64{seedState, 0, 0}
    mrg.x2 = [3]int64{seedState, 0, 0}
}

func (mrg *MRG32k3a) nextInt() int64 {
    x1i := mod(a1[0]*mrg.x1[0]+a1[1]*mrg.x1[1]+a1[2]*mrg.x1[2], m1)
    x2i := mod(a2[0]*mrg.x2[0]+a2[1]*mrg.x2[1]+a2[2]*mrg.x2[2], m2)
    mrg.x1 = [3]int64{x1i, mrg.x1[0], mrg.x1[1]} /* keep last three */
    mrg.x2 = [3]int64{x2i, mrg.x2[0], mrg.x2[1]} /* keep last three */
    return mod(x1i-x2i, m1) + 1
}

func (mrg *MRG32k3a) nextFloat() float64 { return float64(mrg.nextInt()) / float64(d) }

func main() {
    randomGen := MRG32k3aNew()
    randomGen.seed(1234567)
    for i := 0; i < 5; i++ {
        fmt.Println(randomGen.nextInt())
    }

    var counts [5]int
    randomGen.seed(987654321)
    for i := 0; i < 1e5; i++ {
        j := int(math.Floor(randomGen.nextFloat() * 5))
        counts[j]++
    }
    fmt.Println("\nThe counts for 100,000 repetitions are:")
    for i := 0; i < 5; i++ {
        fmt.Printf("  %d : %d\n", i, counts[i])
    }
}

Output:

1459213977
2827710106
4245671317
3877608661
2595287583

The counts for 100,000 repetitions are:
  0 : 20002
  1 : 20060
  2 : 19948
  3 : 20059
  4 : 19931

Haskell

import Data.List

randoms :: Int -> [Int]
randoms seed = unfoldr go ([seed,0,0],[seed,0,0])
  where
    go (x1,x2) =
      let x1i = sum (zipWith (*) x1 a1) `mod` m1
          x2i = sum (zipWith (*) x2 a2) `mod` m2
      in Just $ ((x1i - x2i) `mod` m1, (x1i:init x1, x2i:init x2))
    
    a1 = [0, 1403580, -810728]
    m1 = 2^32 - 209
    a2 = [527612, 0, -1370589]
    m2 = 2^32 - 22853

randomsFloat = map ((/ (2^32 - 208)) . fromIntegral) . randoms

*Main> take 5 $ randoms 1234567
[1459213976,2827710105,4245671316,3877608660,2595287582]

*Main> let hist = map length . group . sort
*Main> hist . take 100000 $ (floor . (*5)) <$> randomsFloat 987654321
[20002,20060,19948,20059,19931]

As a RandomGen instanse

import System.Random

newtype MRG32k3a = MRG32k3a ([Int],[Int])

mkMRG32k3a s = MRG32k3a ([s,0,0],[s,0,0])

instance RandomGen MRG32k3a where
  next (MRG32k3a (x1,x2)) =
    let x1i = sum (zipWith (*) x1 a1) `mod` m1
        x2i = sum (zipWith (*) x2 a2) `mod` m2
    in ((x1i - x2i) `mod` m1, MRG32k3a (x1i:init x1, x2i:init x2))
    where
      a1 = [0, 1403580, -810728]
      m1 = 2^32 - 209
      a2 = [527612, 0, -1370589]
      m2 = 2^32 - 22853

  split _ = error "MRG32k3a is not splittable"

In this case the sequence or numbers differs from direct unfolding, due to internal uniform shuffling.

*Main> take 5 $ randoms (mkMRG32k3a 1234567)
[2827710105,3877608660,3642754129,1259674122,3002249941]

*Main> let hist = map length . group . sort
*Main> hist . take 100000 $ (floor . (*5)) <$> (randoms (mkMRG32k3a 987654321) :: [Float])
[20015,19789,20024,20133,20039]

Java

Translation of: C++

public class App {
    private static long mod(long x, long y) {
        long m = x % y;
        if (m < 0) {
            if (y < 0) {
                return m - y;
            } else {
                return m + y;
            }
        }
        return m;
    }

    public static class RNG {
        // first generator
        private final long[] a1 = {0, 1403580, -810728};
        private static final long m1 = (1L << 32) - 209;
        private long[] x1;
        // second generator
        private final long[] a2 = {527612, 0, -1370589};
        private static final long m2 = (1L << 32) - 22853;
        private long[] x2;
        // other
        private static final long d = m1 + 1;

        public void seed(long state) {
            x1 = new long[]{state, 0, 0};
            x2 = new long[]{state, 0, 0};
        }

        public long nextInt() {
            long x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1);
            long x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2);
            long z = mod(x1i - x2i, m1);

            // keep the last three values of the first generator
            x1 = new long[]{x1i, x1[0], x1[1]};
            // keep the last three values of the second generator
            x2 = new long[]{x2i, x2[0], x2[1]};

            return z + 1;
        }

        public double nextFloat() {
            return 1.0 * nextInt() / d;
        }
    }

    public static void main(String[] args) {
        RNG rng = new RNG();

        rng.seed(1234567);
        System.out.println(rng.nextInt());
        System.out.println(rng.nextInt());
        System.out.println(rng.nextInt());
        System.out.println(rng.nextInt());
        System.out.println(rng.nextInt());
        System.out.println();

        int[] counts = {0, 0, 0, 0, 0};
        rng.seed(987654321);
        for (int i = 0; i < 100_000; i++) {
            int value = (int) Math.floor(rng.nextFloat() * 5.0);
            counts[value]++;
        }
        for (int i = 0; i < counts.length; i++) {
            System.out.printf("%d: %d%n", i, counts[i]);
        }
    }
}

Output:

Julia

const a1 = [0, 1403580, -810728]
const m1 = 2^32 - 209
const a2 = [527612, 0, -1370589]
const m2 = 2^32 - 22853
const d = m1 + 1

mutable struct MRG32k3a
    x1::Tuple{Int64, Int64, Int64}
    x2::Tuple{Int64, Int64, Int64}
    MRG32k3a() = new((0, 0, 0), (0, 0, 0))
    MRG32k3a(seed_state) = new((seed_state, 0, 0), (seed_state, 0, 0))
end
seed(sd) = begin @assert(0 < sd < d); MRG32k3a(sd) end

function next_int(x::MRG32k3a)
    x1i = mod1(a1[1] * x.x1[1] + a1[2] * x.x1[2] + a1[3] * x.x1[3], m1)
    x2i = mod1(a2[1] * x.x2[1] + a2[2] * x.x2[2] + a2[3] * x.x2[3], m2)
    x.x1 = (x1i, x.x1[1], x.x1[2])
    x.x2 = (x2i, x.x2[1], x.x2[2])
    return mod1(x1i - x2i, m1) + 1
end

next_float(x::MRG32k3a) = next_int(x) / d

const g1 = seed(1234567)
for _ in 1:5
    println(next_int(g1))
end
const g2 = seed(987654321)
hist = fill(0, 5)
for _ in 1:100_000
    hist[Int(floor(next_float(g2) * 5)) + 1] += 1
end
foreach(p -> print(p[1], ": ", p[2], "  "), enumerate(hist))

Output:

1459213977
2827710106
4245671317
3877608661
2595287583
1: 20002  2: 20060  3: 19948  4: 20059  5: 19931

Kotlin

Translation of: C++

import kotlin.math.floor

fun mod(x: Long, y: Long): Long {
    val m = x % y
    return if (m < 0) {
        if (y < 0) {
            m - y
        } else {
            m + y
        }
    } else m
}

class RNG {
    // first generator
    private val a1 = arrayOf(0L, 1403580L, -810728L)
    private val m1 = (1L shl 32) - 209
    private var x1 = arrayOf(0L, 0L, 0L)

    // second generator
    private val a2 = arrayOf(527612L, 0L, -1370589L)
    private val m2 = (1L shl 32) - 22853
    private var x2 = arrayOf(0L, 0L, 0L)

    private val d = m1 + 1

    fun seed(state: Long) {
        x1 = arrayOf(state, 0, 0)
        x2 = arrayOf(state, 0, 0)
    }

    fun nextInt(): Long {
        val x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1)
        val x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2)
        val z = mod(x1i - x2i, m1)

        // keep last three values of the first generator
        x1 = arrayOf(x1i, x1[0], x1[1])
        // keep last three values of the second generator
        x2 = arrayOf(x2i, x2[0], x2[1])

        return z + 1
    }

    fun nextFloat(): Double {
        return nextInt().toDouble() / d
    }
}

fun main() {
    val rng = RNG()

    rng.seed(1234567)
    println(rng.nextInt())
    println(rng.nextInt())
    println(rng.nextInt())
    println(rng.nextInt())
    println(rng.nextInt())
    println()

    val counts = IntArray(5)
    rng.seed(987654321)
    for (i in 0 until 100_000) {
        val v = floor((rng.nextFloat() * 5.0)).toInt()
        counts[v]++
    }
    for (iv in counts.withIndex()) {
        println("${iv.index}: ${iv.value}")
    }
}

Output:

Nim

import algorithm, math, sequtils, strutils, tables

const
  # First generator.
  a1 = [int64 0, 1403580, -810728]
  m1: int64 = 2^32 - 209
  # Second generator.
  a2 = [int64 527612, 0, -1370589]
  m2: int64 = 2^32 - 22853

  d = m1 + 1

type MRG32k3a = object
  x1: array[3, int64]   # List of three last values of gen #1.
  x2: array[3, int64]   # List of three last values of gen #2.


func seed(gen: var MRG32k3a; seedState: int64) =
  assert seedState in 1..<d
  gen.x1 = [seedState, 0, 0]
  gen.x2 = [seedState, 0, 0]

func nextInt(gen: var MRG32k3a): int64 =
  let x1i = floormod(a1[0] * gen.x1[0] + a1[1] * gen.x1[1] + a1[2] * gen.x1[2], m1)
  let x2i = floormod(a2[0] * gen.x2[0] + a2[1] * gen.x2[1] + a2[2] * gen.x2[2], m2)
  # In version 1.4, the following two lines doesn't work.
  # gen.x1 = [x1i, gen.x1[0], gen.x1[1]]    # Keep last three.
  # gen.x2 = [x2i, gen.x2[0], gen.x2[1]]    # Keep last three.
  gen.x1[2] = gen.x1[1]; gen.x1[1] = gen.x1[0]; gen.x1[0] = x1i
  gen.x2[2] = gen.x2[1]; gen.x2[1] = gen.x2[0]; gen.x2[0] = x2i
  result = floormod(x1i - x2i, m1) + 1

func nextFloat(gen: var MRG32k3a): float =
  gen.nextInt().float / d.float

when isMainModule:
  var gen: MRG32k3a

  gen.seed(1234567)
  for _ in 1..5:
    echo gen.nextInt()

  echo ""
  gen.seed(987654321)
  var counts: CountTable[int]
  for _ in 1..100_000:
    counts.inc int(gen.nextFloat() * 5)
  echo sorted(toSeq(counts.pairs)).mapIt($it[0] & ": " & $it[1]).join(", ")

Output:

1459213977
2827710106
4245671317
3877608661
2595287583

0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931

Pretty straightforward translation from the directions. Used column/vector multiplication (essentially he dot product) instead of the more tedious form given in the definition of x1i and x2i; rationals (t_FRAC) used in place of floating-point since GP lacks floating-point.

a1 = [0, 1403580, -810728];
m1 = 2^32-209;
a2 = [527612, 0, -1370589];
m2 = 2^32-22853;
d = m1+1;
seed(s)=x1=x2=[s,0,0];
next_int()=
{
  my(x1i=a1*x1~%m1, x2i=a2*x2~%m2);
  x1 = [x1i, x1[1], x1[2]];
  x2 = [x2i, x2[1], x2[2]];
  (x1i-x2i)%m1 + 1;
}
next_float()=next_int()/d;

seed(1234567);
vector(5,i,next_int())
seed(987654321);
v=vector(5); for(i=1,1e5, v[next_float()*5\1+1]++); v

Output:

%1 = [1459213977, 2827710106, 4245671317, 3877608661, 2595287583]
%2 = [20002, 20060, 19948, 20059, 19931]

Perl

use strict;
use warnings;
use feature 'say';

package MRG32k3a {

    use constant {
        m1 => 2**32 - 209,
        m2 => 2**32 - 22853
    };

    use Const::Fast;
    const my @a1 => <     0 1403580  -810728>;
    const my @a2 => <527612       0 -1370589>;

    sub new {
        my ($class,undef,$seed) = @_;
        my @x1 = my @x2 = ($seed, 0, 0);
        bless {x1 => \@x1, x2 => \@x2}, $class;
    }

    sub next_int {
        my ($self) = @_;
        unshift @{$$self{x1}}, ($a1[0] * $$self{x1}[0] + $a1[1] * $$self{x1}[1] + $a1[2] * $$self{x1}[2]) % m1; pop @{$$self{x1}};
        unshift @{$$self{x2}}, ($a2[0] * $$self{x2}[0] + $a2[1] * $$self{x2}[1] + $a2[2] * $$self{x2}[2]) % m2; pop @{$$self{x2}};
        ($$self{x1}[0] - $$self{x2}[0]) % (m1 + 1)
    }

    sub next_float { $_[0]->next_int / (m1 + 1) }
}

say 'Seed: 1234567, first 5 values:';
my $rng = MRG32k3a->new( seed => 1234567 );
say $rng->next_int for 1..5;

my %h;
say "\nSeed: 987654321, values histogram:";
$rng = MRG32k3a->new( seed => 987654321 );
$h{int 5 * $rng->next_float}++ for 1..100_000;
say "$_ $h{$_}" for sort keys %h;

Output:

Seed: 1234567, first 5 values:
1459213977
2827710106
4245671317
3877608661
2595287583

Seed: 987654321, values histogram:
0 20002
1 20060
2 19948
3 20059
4 19931

Phix

with javascript_semantics
constant
    -- First generator
    a1 = {0, 1403580, -810728},
    m1 = power(2,32) - 209,
    -- Second Generator
    a2 = {527612, 0, -1370589},
    m2 = power(2,32) - 22853,
     d = m1 + 1
 
sequence x1 = {0, 0, 0},  /* list of three last values of gen #1 */
         x2 = {0, 0, 0}   /* list of three last values of gen #2 */
 
procedure seed(integer seed_state)
    assert(seed_state>0 and seed_state<d)
    x1 = {seed_state, 0, 0}
    x2 = {seed_state, 0, 0}
end procedure
 
function next_int()
    atom x1i = mod(a1[1]*x1[1] + a1[2]*x1[2] + a1[3]*x1[3],m1),
         x2i = mod(a2[1]*x2[1] + a2[2]*x2[2] + a2[3]*x2[3],m2)
    x1 = {x1i, x1[1], x1[2]}    /* Keep last three */
    x2 = {x2i, x2[1], x2[2]}    /* Keep last three */
    atom z = mod(x1i-x2i,m1),
    answer = (z + 1)
    return answer
end function
 
function next_float()
    return next_int() / d
end function
 
seed(1234567)
for i=1 to 5 do
    printf(1,"%d\n",next_int())
end for
seed(987654321)
sequence r = repeat(0,5)
for i=1 to 100_000 do
    integer rdx = floor(next_float()*5)+1
    r[rdx] += 1
end for
?r

Output:

1459213977
2827710106
4245671317
3877608661
2595287583
{20002,20060,19948,20059,19931}

Python

# Constants
a1 = [0, 1403580, -810728]
m1 = 2**32 - 209
#
a2 = [527612, 0, -1370589]
m2 = 2**32 - 22853
#
d = m1 + 1

class MRG32k3a():
    
    def __init__(self, seed_state=123):
        self.seed(seed_state)
    
    def seed(self, seed_state):
        assert 0 <seed_state < d, f"Out of Range 0 x < {d}"
        self.x1 = [seed_state, 0, 0]
        self.x2 = [seed_state, 0, 0]
        
    def next_int(self):
        "return random int in range 0..d"
        x1i = sum(aa * xx  for aa, xx in zip(a1, self.x1)) % m1
        x2i = sum(aa * xx  for aa, xx in zip(a2, self.x2)) % m2
        self.x1 = [x1i] + self.x1[:2]
        self.x2 = [x2i] + self.x2[:2]

        z = (x1i - x2i) % m1
        answer = (z + 1)
        
        return answer
    
    def  next_float(self):
        "return random float between 0 and 1"
        return self.next_int() / d
    

if __name__ == '__main__':
    random_gen = MRG32k3a()
    random_gen.seed(1234567)
    for i in range(5):
        print(random_gen.next_int())
        
    random_gen.seed(987654321)
    hist = {i:0 for i in range(5)}
    for i in range(100_000):
        hist[int(random_gen.next_float() *5)] += 1
    print(hist)

Output:

1459213977
2827710106
4245671317
3877608661
2595287583
{0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931}

Raku

Works with: Rakudo version 2020.07

Translation of: Python

All constants are encapsulated within the class.

class MRG32k3a {
    has @!x1;
    has @!x2;

    constant a1 = 0, 1403580, -810728;
    constant a2 = 527612, 0, -1370589;
    constant m1 = 2**32 - 209;
    constant m2 = 2**32 - 22853;

    submethod BUILD ( Int :$seed where 0 < * <= m1 = 1 ) { @!x1 = @!x2 = $seed, 0, 0 }

    method next-int {
        @!x1.unshift: (a1[0] * @!x1[0] + a1[1] * @!x1[1] + a1[2] * @!x1[2]) % m1; @!x1.pop;
        @!x2.unshift: (a2[0] * @!x2[0] + a2[1] * @!x2[1] + a2[2] * @!x2[2]) % m2; @!x2.pop;
        (@!x1[0] - @!x2[0]) % (m1 + 1)
    }

    method next-rat { self.next-int / (m1 + 1) }
}


# Test next-int with custom seed
say 'Seed: 1234567; first five Int values:';
my $rng = MRG32k3a.new :seed(1234567);
.say for $rng.next-int xx 5;


# Test next-rat (since these are rational numbers by default)
say "\nSeed: 987654321; first 1e5 Rat values histogram:";
$rng = MRG32k3a.new :seed(987654321);
say ( ($rng.next-rat * 5).floor xx 100_000 ).Bag;


# Test next-int with default seed
say "\nSeed: default; first five Int values:";
$rng = MRG32k3a.new;
.say for $rng.next-int xx 5;

Output:

Seed: 1234567; first five Int values:
1459213977
2827710106
4245671317
3877608661
2595287583

Seed: 987654321; first 1e5 Rat values histogram:
Bag(0(20002) 1(20060) 2(19948) 3(20059) 4(19931))

Seed: default; first five Int values:
4294439476
798392476
1012402088
1268414424
3353586348

Ruby

Translation of: C

def mod(x, y)
    m = x % y
    if m < 0 then
        if y < 0 then
            return m - y
        else
            return m + y
        end
    end
    return m
end

# Constants
# First generator
A1 = [0, 1403580, -810728]
A1.freeze
M1 = (1 << 32) - 209
# Second generator
A2 = [527612, 0, -1370589]
A2.freeze
M2 = (1 << 32) - 22853

D = M1 + 1

# the last three values of the first generator
$x1 = [0, 0, 0]
# the last three values of the second generator
$x2 = [0, 0, 0]

def seed(seed_state)
    $x1 = [seed_state, 0, 0]
    $x2 = [seed_state, 0, 0]
end

def next_int()
    x1i = mod((A1[0] * $x1[0] + A1[1] * $x1[1] + A1[2] * $x1[2]), M1)
    x2i = mod((A2[0] * $x2[0] + A2[1] * $x2[1] + A2[2] * $x2[2]), M2)
    z = mod(x1i - x2i, M1)

    $x1 = [x1i, $x1[0], $x1[1]]
    $x2 = [x2i, $x2[0], $x2[1]]

    return z + 1
end

def next_float()
    return 1.0 * next_int() / D
end

########################################

seed(1234567)
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print "\n"

counts = [0, 0, 0, 0, 0]
seed(987654321)
for i in 1 .. 100000
    value = (next_float() * 5.0).floor
    counts[value] = counts[value] + 1
end
counts.each_with_index { |v,i|
    print i, ": ", v, "\n"
}

Output:

Mimicking the Pseudo-code

# Constants
# First generator
A1 = [0, 1403580, -810728]
M1 = 2**32 - 209
# Second Generator
A2 = [527612, 0, -1370589]
M2 = 2**32 - 22853
    
D = M1 + 1
   
class MRG32k3a

  def seed(seed_state)
    raise ArgumentError unless seed_state.between?(0, D) 
    @x1 = [seed_state, 0, 0]
    @x2 = [seed_state, 0, 0]
  end 
           
  def next_int
    x1i = (A1[0]*@x1[0] + A1[1]*@x1[1] + A1[2]*@x1[2]).modulo M1
    x2i = (A2[0]*@x2[0] + A2[1]*@x2[1] + A2[2]*@x2[2]).modulo M2
    @x1 = [x1i, @x1[0], @x1[1]]   # Keep last three
    @x2 = [x2i, @x2[0], @x2[1]]   # Keep last three
    z   = (x1i - x2i) % M1
    return z + 1
  end
       
  def next_float
    next_int.to_f / D
  end 
       
end

random_gen = MRG32k3a.new
random_gen.seed(1234567)
5.times{ puts random_gen.next_int}

random_gen = MRG32k3a.new
random_gen.seed(987654321)
p 100_000.times.map{(random_gen.next_float() * 5).floor}.tally.sort.to_h

Sidef

Translation of: Perl

class MRG32k3a(seed) {

    define(
        m1 = (2**32 - 209)
        m2 = (2**32 - 22853)
    )

    define(
        a1 = %n<     0 1403580  -810728>
        a2 = %n<527612       0 -1370589>
    )

    has x1 = [seed, 0, 0]
    has x2 = x1.clone

    method next_int {
        x1.unshift(a1.map_kv {|k,v| v * x1[k] }.sum % m1); x1.pop
        x2.unshift(a2.map_kv {|k,v| v * x2[k] }.sum % m2); x2.pop
        (x1[0] - x2[0]) % (m1 + 1)
    }

    method next_float { self.next_int / (m1 + 1) -> float }
}

say "Seed: 1234567, first 5 values:"
var rng = MRG32k3a(seed: 1234567)
5.of { rng.next_int }.each { .say }

say "\nSeed: 987654321, values histogram:";
var rng = MRG32k3a(seed: 987654321)
var freq = 100_000.of { rng.next_float * 5 -> int }.freq
freq.sort.each_2d {|k,v| say "#{k} #{v}" }

Output:

Seed: 1234567, first 5 values:
1459213977
2827710106
4245671317
3877608661
2595287583

Seed: 987654321, values histogram:
0 20002
1 20060
2 19948
3 20059
4 19931

uBasic/4tH

Works with: v3.64

Translation of: C

Since uBasic/4tH has no floating point support, only the integer part of the task can be implemented.

@(0) = 0                               ' First generator
@(1) = 1403580
@(2) = -810728
m = SHL(1, 32) - 209

@(3) = 527612                          ' Second generator
@(4) = 0
@(5) = -1370589
n = SHL(1, 32) - 22853

d = SHL(1, 32) - 209 + 1               ' m + 1

Proc  _Seed(1234567)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print
End

_Mod Param(2)
  Local(1)
  c@ = a@ % b@
  If c@ < 0 Then
     If b@ < 0 Then
       Return (c@-b@)
     Else
       Return (c@+b@)
     Endif
  EndIf
Return (c@)

_Seed Param(1)                         ' seed the PRNG
  @(6) = a@
  @(7) = 0
  @(8) = 0

  @(9) = a@
  @(10) = 0
  @(11) = 0
Return

_NextInt                               ' get the next random integer value
  Local(3)

  a@ = FUNC(_Mod((@(0) * @(6) + @(1) *  @(7) + @(2) *  @(8)), m))
  b@ = FUNC(_Mod((@(3) * @(9) + @(4) * @(10) + @(5) * @(11)), n))
  c@ = FUNC(_Mod(a@ - b@, m))

  ' keep last three values of the first generator
  @(8) = @(7)
  @(7) = @(6)
  @(6) = a@

  ' keep last three values of the second generator
  @(11) = @(10)
  @(10) = @(9)
  @(9) = b@

Return (c@ + 1)

Output:

1459213977
2827710106
4245671317
3877608661
2595287583


0 OK, 0:398

Wren

Translation of: Python

// constants
var A1 = [0, 1403580, -810728]
var M1 = 2.pow(32) - 209
var A2 = [527612, 0, -1370589]
var M2 = 2.pow(32) - 22853
var D = M1 + 1

// Python style modulus
var Mod = Fn.new { |x, y|
    var m = x % y
    return (m < 0) ?  m + y.abs : m
}

class MRG32k3a {
    construct new() {
        _x1 = [0, 0, 0]
        _x2 = [0, 0, 0]
    }

    seed(seedState) {
        if (seedState <= 0 || seedState >= D) {
            Fiber.abort("Argument must be in the range [0, %(D)].")
        }
        _x1 = [seedState, 0, 0]
        _x2 = [seedState, 0, 0]
    }

    nextInt {
        var x1i = Mod.call(A1[0]*_x1[0] + A1[1]*_x1[1] + A1[2]*_x1[2], M1)
        var x2i = Mod.call(A2[0]*_x2[0] + A2[1]*_x2[1] + A2[2]*_x2[2], M2)
        _x1 = [x1i, _x1[0], _x1[1]]    /* keep last three */
        _x2 = [x2i, _x2[0], _x2[1]]    /* keep last three */
        return Mod.call(x1i - x2i, M1) + 1
    }

    nextFloat { nextInt / D }
}

var randomGen = MRG32k3a.new()
randomGen.seed(1234567)
for (i in 0..4) System.print(randomGen.nextInt)

var counts = List.filled(5, 0)
randomGen.seed(987654321)
for (i in 1..1e5) {
    var i = (randomGen.nextFloat * 5).floor
    counts[i] = counts[i] + 1
}
System.print("\nThe counts for 100,000 repetitions are:")
for (i in 0..4) System.print("  %(i) : %(counts[i])")

Output:

1459213977
2827710106
4245671317
3877608661
2595287583

The counts for 100,000 repetitions are:
  0 : 20002
  1 : 20060
  2 : 19948
  3 : 20059
  4 : 19931