Pseudo-random numbers/Combined recursive generator MRG32k3a

Pseudo-random numbers/Combined recursive generator MRG32k3a is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

<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

Translation of: Python

<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

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

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

Works with: Rakudo version 2020.07

Translation of: Python

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

}

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;

1. 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;

1. 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

Wren

Translation of: Python

<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