Pseudo-random numbers/Combined recursive generator MRG32k3a
From Rosetta Code
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.
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.
Factor[edit]
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 ] [email protected] - 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 } }
Go[edit]
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
Julia[edit]
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
Perl[edit]
use strict;
use warnings;
use feature 'say';
package MRG32k3a {
our(@x1,@x2);
our $m1 = 2**32 - 209;
our $m2 = 2**32 - 22853;
our @a1 = < 0 1403580 -810728>;
our @a2 = <527612 0 -1370589>;
sub new { @x1 = @x2 = ($_[2], 0, 0) }
sub next_int {
unshift @x1, ($a1[0] * $x1[0] + $a1[1] * $x1[1] + $a1[2] * $x1[2]) % $m1; pop @x1;
unshift @x2, ($a2[0] * $x2[0] + $a2[1] * $x2[1] + $a2[2] * $x2[2]) % $m2; pop @x2;
($x1[0] - $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;
- 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[edit]
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
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 {20002,20060,19948,20059,19931}
Python[edit]
# 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[edit]
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
Wren[edit]
// 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