# Pseudo-random numbers/Combined recursive generator MRG32k3a

Pseudo-random numbers/Combined recursive generator MRG32k3a
You are encouraged to solve this task according to the task description, using any language you may know.
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

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 */

```
• 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
```

## 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]
```

```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;
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;
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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## 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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## 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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## 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 ;

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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## 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
```

```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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## jq

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

foreach range(0; 5) as \$i (seed(1234567); nextInt ) | .nextInt;

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)");

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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

## 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```

## Pari/GP

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),
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

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:
```1459213977
2827710106
4245671317
3877608661
2595287583

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

### 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
```