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

Content added Content deleted

Inline

Latest revision as of 09:46, 11 June 2024

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:

jq

Adapted from Wren

Works with jq, the C implementation of jq

Works with gojq, the Go implementation of jq

Works with jaq, the Rust implementation of jq provided the MRG32k3a module (minus its declaration) is inlined.

The jq module MRG32k3a is available at MRG32k3a.jq.

include "MRG32k3a" {search: "."};  # see above

def task1:
  foreach range(0; 5) as $i (seed(1234567); nextInt ) | .nextInt;

def task2($n):
  seed(987654321)
  | reduce range(0; $n) as $i (.counts = [range(0;5)|0];
      nextFloat
      | .counts[ (.nextFloat * 5) | floor ] += 1 )
 | "\nThe counts for \($n) repetitions are:",
   (range(0;5) as $i | "  \($i) : \(.counts[$i] // 0)");

task1,
task2(100000)

Output:

1459213977
2827710106
4245671317
3877608661
2595287583

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

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