Pseudo-random numbers/Combined recursive generator MRG32k3a: Difference between revisions
uBasic/4tH version added |
m Removed some superfluous semi-colons |
||
Line 1,184: | Line 1,184: | ||
' keep last three values of the first generator |
' keep last three values of the first generator |
||
@(8) = @(7) |
@(8) = @(7) |
||
@(7) = @(6) |
@(7) = @(6) |
||
@(6) = a@ |
@(6) = a@ |
||
' keep last three values of the second generator |
' keep last three values of the second generator |
||
@(11) = @(10) |
@(11) = @(10) |
||
@(10) = @(9) |
@(10) = @(9) |
||
@(9) = b@ |
@(9) = b@ |
||
Return (c@ + 1)</lang> |
Return (c@ + 1)</lang> |
Revision as of 09:44, 25 June 2021
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
<lang 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)</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 [0 = 20002, 1 = 20060, 2 = 19948, 3 = 20059, 4 = 19931]
Ada
<lang 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; </lang>
<lang Ada> 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; </lang>
<lang Ada>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; </lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
C
<lang 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;
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
C++
<lang cpp>#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;
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Factor
<lang 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 .</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 H{ { 0 20002 } { 1 20060 } { 2 19948 } { 3 20059 } { 4 19931 } }
Go
<lang 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]) }
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 The counts for 100,000 repetitions are: 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
Java
<lang 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]); } }
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Julia
<lang 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))
</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 1: 20002 2: 20060 3: 19948 4: 20059 5: 19931
Kotlin
<lang scala>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}") }
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Nim
<lang 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(", ")</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931
Perl
<lang 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 ($class,$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 { next_int(@_) / (m1 + 1) }
}
say 'Seed: 1234567, first 5 values:'; my $rng = MRG32k3a->new( seed => 1234567 ); say MRG32k3a->next_int($rng) for 1..5;
my %h; say "\nSeed: 987654321, values histogram:"; $rng = MRG32k3a->new( seed => 987654321 ); $h{int 5 * MRG32k3a->next_float($rng)}++ for 1..100_000; say "$_ $h{$_}" for sort keys %h;</lang>
- 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
<lang Phix>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
r[floor(next_float()*5)+1] += 1
end for ?r</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 {20002,20060,19948,20059,19931}
Python
<lang 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)</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 {0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931}
Raku
All constants are encapsulated within the class.
<lang perl6>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;</lang>
- 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
<lang 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"
}</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 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. <lang>@(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)</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0 OK, 0:398
Wren
<lang ecmascript>// 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])")</lang>
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 The counts for 100,000 repetitions are: 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931