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
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
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:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
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
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 ;
! 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
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
FreeBASIC
Const As Longint a1_1 = 0
Const As Longint a1_2 = 1403580
Const As Longint a1_3 = -810728
Const As Longint m1 = 4294967087 ' 2^32 - 209
Const As Longint a2_1 = 527612
Const As Longint a2_2 = 0
Const As Longint a2_3 = -1370589
Const As Longint m2 = 4294944443 ' 2^32 - 22853
Const d As Double = m1 + 1
Dim Shared As Longint x1(2), x2(2)
Sub seed(seed_state As Longint)
If seed_state <= 0 Or seed_state >= m1 Then
Print "Error: Invalid seed state"
End
End If
x1(0) = seed_state : x1(1) = 0 : x1(2) = 0
x2(0) = seed_state : x2(1) = 0 : x2(2) = 0
End Sub
Function next_int() As Longint
Dim As Longint x1i = (a1_1 * x1(0) + a1_2 * x1(1) + a1_3 * x1(2)) Mod m1
If x1i < 0 Then x1i += m1
Dim As Longint x2i = (a2_1 * x2(0) + a2_2 * x2(1) + a2_3 * x2(2)) Mod m2
If x2i < 0 Then x2i += m2
x1(2) = x1(1) : x1(1) = x1(0) : x1(0) = x1i
x2(2) = x2(1) : x2(1) = x2(0) : x2(0) = x2i
Dim As Longint z = (x1i - x2i) Mod m1
If z < 0 Then z += m1
Return z + 1
End Function
Function next_float() As Double
Return next_int() / d
End Function
' Main program
Dim As Integer i, rdx
seed(1234567)
For i = 1 To 5
Print next_int()
Next i
Print
seed(987654321)
Dim As Longint r(4)
For i = 1 To 100000
rdx = Int(next_float() * 5)
r(rdx) += 1
Next i
For i = 0 To 4
Print i & ": " & r(i);
If i < 4 Then Print ", ";
Next i
Sleep
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931
Go
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
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
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
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), 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
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
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
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
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
// 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