An RNG that runs in your mind

People are notoriously unhealthy at developing with random numbers. I needed to have the ability to rapidly generate “random sufficient” numbers. I’m not searching for something that nice, I simply need to have the ability to provide you with the random digits in half a minute. Some wanting round introduced me to an old usenet post by George Marsaglia:

Select a 2-digit quantity, say 23, your “seed”.

Type a brand new 2-digit quantity:
the ten’s digit plus 6 occasions the items digit.

The instance sequence is
23 –> 20 –> 02 –> 12 –> 13 –> 19 –> 55 –> 35 –> …

and its interval is the order of the multiplier, 6, within the group of
residues comparatively prime to the modulus, 10. (59 on this case).

The “random digits” are the items digits of the 2-digit numbers,
ie, 3,0,2,2,3,9,5,… the sequence mod 10.

Marsaglia is most well-known for the diehard suite of RNG exams, so he is aware of his stuff. I’m curious why this works and why he selected 6.

We’re going to make use of Raku, the language for gremlins. I’ll be explaining all of the bizarre options I exploit in dropdowns in case you’re a little bit of a gremlin, too.

The sequence is periodic, that means that if we iteratively apply it we’ll ultimately get the identical factor. Let’s begin with a operate (“subroutine”) that produces the entire sequence:

my sub next-rng(Int $begin, Int $unitmult = 6, --> Listing) {
    my @out;
    my $subsequent = $begin;
    repeat whereas $subsequent !(elem) @out {
        @out.append($subsequent);
        $subsequent = sum($subsequent.polymod(10) Z* $unitmult,1);
    };
    return @out;
}

As soon as we’ve got a sequence we are able to print it with the put or say instructions, which have subtly different behaviors I’m not going to get into.

01 06 36 39 57 47 46 40 04 24 26 38 51 11 07 42 16 
37 45 34 27 44 28 50 05 30 03 18 49 58 53 23 20 02 
12 13 19 55 35 33 21 08 48 52 17 43 22 14 25 32 15 
31 09 54 29 56 41 10 01

Keep in mind, the random numbers are the final digit. So the RNG goes 1 -> 6 -> 6 -> 9 -> 7 -> 7 -> …

Investigating Properties

If the RNG is uniform then every digit ought to seem within the sequence the identical variety of occasions. We will examine this by casting the final digits to a multiset, or “bag”.

say bag(next-rng(1) <<%>> 10);

Bag(0(5) 1(6) 2(6) 3(6) 4(6) 5(6) 6(6) 7(6) 8(6) 9(5))

That appears to be a uniform-enough distribution, although I’m a bit much less prone to get a 0 or a 9.

My subsequent concept comes from the diehard exams. From the wiki:

Overlapping permutations: Analyze sequences of 5 consecutive random numbers. The 120 attainable orderings ought to happen with statistically equal likelihood.

There are solely 54 5-number sequences within the dataset, so I’ll as a substitute apply this to 2-number “transitions”. I’ll do that by outputting a 10-by-10 grid the place the (i, j)th index (from 0) is the variety of transitions from last-digit i to j. For instance, the sequence consists of the transition 28 -> 50, and no different transitions of type X8 -> Y0, so cell (8, 0) must be a 1.

sub successions-grid(@orbit) {
  my @pairs = (|@orbit , @orbit[0]).map(* % 10).rotor(2 => -1);
  for ^10 -> $x {put ($x X ^10).map({@pairs.grep($_).elems})}
}

> successions-grid(next-rng(1, 6))

0 1 1 1 1 1 0 0 0 0
1 1 0 0 0 0 1 1 1 1
0 0 1 1 1 1 1 1 0 0
1 1 1 1 0 0 0 0 1 1
0 0 0 0 1 1 1 1 1 1
1 1 1 1 1 1 0 0 0 0
1 1 0 0 0 0 1 1 1 1
0 0 1 1 1 1 1 1 0 0
1 1 1 1 0 0 0 0 1 1
0 0 0 0 1 1 1 1 1 0

We will see from this desk that some transitions are inconceivable. If I generate a 0, I can’t get a 6 proper after. Clearly not an awesome RNG, however my expectations have been fairly low anyway.

Why 6?

What if as a substitute of multiplying the final digit by 6, I multiply by 4?

> say next-rng(1, 4);
01 04 16 25 22 10 

I dunno, I kinda like an RNG that by no means provides me 3. The distinct sequences are known as orbits and their lengths are known as durations. Let’s see all of the attainable orbits we are able to get by utilizing 4 because the multiplier:

sub orbits-for-mod(int $mult, $prime = 20) {
  my &f = &next-rng.assuming(*, $mult);
  (1..$prime).map(&f).distinctive(as => &set)
}

> say orbits-for-mod(4, 38).map(*.gist).be part of("n");
[1 4 16 25 22 10]
[2 8 32 11 5 20]
[3 12 9 36 27 30]
[6 24 18 33 15 21]
[7 28 34 19 37 31]
[13]
[14 17 29 38 35 23]
[26]

In case you picked 13 as your beginning worth, your random digits could be 3, 3, 3, 3, 3, 3.

For apparent causes, 4 ought to by no means be our multiplier. Actually for a multiplier to offer a “good” RNG, it must have precisely one orbit.

> say (1..30).grep(*.&orbits-for-mod == 1)
(3 6 11 15 18 23 27)

Rationalization

See Also

What occurred with the Net Monetization API? – Chris Coyier

• .& applies a top-level routine as a way. grep(*.&f) is the equal of grep({f($_)}).
• &orbits-for-mod returns a listing. == coerces each inputs to numbers, and coercing a listing to a quantity returns the variety of parts. So we’re testing if the returned listing has one factor, ie there’s precisely one orbit. (In case you don’t need to examine with out coercion, use both === or eqv.)

This fashion of doing issues is fairly gradual and likewise solely appears for orbits that begin with a quantity as much as 20. So it will miss the 26 -> 26 orbit for x=4. We’ll repair each of those points later.

So some “good” selections for n are 6, 11, and 18.

Be aware that if you find yourself with a 3 digit quantity, you deal with the primary two digits as a single quantity. For n=11, 162 results in 16 + 22, not 6 + 22 (or 6 + 1 + 22).

Why does this work?

Right here’s part of the reason that actually confused me:

and its interval is the order of the multiplier, 6, within the group of
residues comparatively prime to the modulus, 10. (59 on this case).

After speaking with some pals and a number of studying Wiki articles, it began making extra sense. I’m mentally computing a “multiply with carry” RNG with constants a=x and c=10. This alternative has a cool property: if MWC(x) = y, then 10y mod (10n-1) = x!

MWC:    01 -> 06 -> 36 -> ... -> 41 -> 10 -> 01
10ypercent59: 01 -> 10 -> 41 -> ... -> 36 -> 06 -> 01

That’s fairly neat! It’s simpler for me to mathematically purpose about 10y mod 59 than “multiply the final digit by six and add the primary digit”. For instance, it’s clear why the RNG generates 0 and 9 barely much less usually than the opposite digits: regardless of which multiplier we decide, the generated sequence will go from 1 to 10n-2, “leaving out” 10n-1 (which ends with 9) and 10n (ends with 0).

“Multiply and modulo” is also called the Lehmer RNG.

Discovering higher RNGs

So what different numbers work? We already know {that a} good multiplier will produce just one orbit, and I confirmed some code above for calculating that. Sadly, it’s an O(n²) worst-case algorithm. Fascinated with the MWC algorithm as “Lehmer in reverse” provides us a greater methodology: if n is an efficient multiplier, then the interval of the orbit ranging from 1 must be 10n-2.

The Lehmer strategy additionally provides us a quicker means of computing the orbit:

sub oneorbit(x) {
  10, * *10% (10*x - 1) … 1
}

Rationalization

Then we examine the interval utilizing the identical “== coerces lists to lengths” trick as earlier than.

> say (1..100).grep({oneorbit($_) == 10*$_-2});

(2 3 6 11 15 18 23 27 38 39 42 50 51 62 66 71)

I can see why Marsaglia selected 6: most programmers know their 6 times-table and it by no means returns a 3-digit quantity, so the addition step is actual simple. The orbit has solely 58 numbers and also you received’t get some digit sequences, but when you could pull out a number of random digits rapidly it’s completely nice.

If you would like extra randomness, I see a few candidates. 50 has a interval of 498 and is extremely simple to compute. If the ultimate digit is even then you definately don’t have to do any carries: 238 -> 423!

That stated, the 50-sequence doesn’t appear as random as different sequences. There’s some extent the place it generates 9 even numbers adopted by 8 odd ones. Don’t use it to simulate coin flips.

The final attention-grabbing quantity is eighteen. It has a decent interval of 178 and has each attainable digit transition:

> successions-grid(next-rng(1, 18))

1 2 2 2 2 2 2 2 1 1
2 2 2 2 2 2 1 1 2 2
2 2 2 2 1 1 2 2 2 2
2 2 1 1 2 2 2 2 2 2
1 1 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 1 1
2 2 2 2 2 2 1 1 2 2
2 2 2 2 1 1 2 2 2 2
2 2 1 1 2 2 2 2 2 2
1 1 2 2 2 2 2 2 2 1

The draw back is that you must be taught the 18 times-table. This isn’t too unhealthy: I internalized it with perhaps 10 minutes of observe. I’m nonetheless not nice at doing the entire MWC step however I can constantly produce one other random digit each 5 seconds or so. That’s adequate for me.

You may see the Raku code I used to analysis this here. It’s arrange as a CLI so you need to use it in a number of alternative ways; see the file for more information.

Due to Codesections for suggestions and Quinn Wilton and Jeremy Kun for serving to me perceive the maths.