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

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 } = 1 ) { @!x1 = @!x2 = $seed, 0, 0; }

   method next-int {
       my $x1i = (a1[0] * @!x1[0] + a1[1] * @!x1[1] + a1[2] * @!x1[2]) % m1;
       my $x2i = (a2[0] * @!x2[0] + a2[1] * @!x2[1] + a2[2] * @!x2[2]) % m2;
       @!x1 = $x1i, |@!x1[0, 1];
       @!x2 = $x2i, |@!x2[0, 1];
       ($x1i - $x2i) % 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 and increment

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
1012402089
1268414425
3353586348