## Random number generator that can be run in the head

Humans are very bad at coming up with random numbers. I wanted to learn how to generate “fairly random” numbers quickly. I didn’t need anything fancy, just a way to make up random numbers in half a minute. After searching online, I found it

old Usenet post

written by George Marsaglia:

Choose a two-digit number, say 23. This will be your seed value.

Create a new two-digit number: the number of tens plus six times the number of ones.

Example sequence: 23 -> (2 + 6 * 3) = 20 -> (2 + 6 * 0) = 02 -> 12 -> 13 -> 19 -> 55 -> 35 -> …

Its period will be the order of the multiplier (6) of the group of residues prime to the module, 10 (in the case of 59).

The “random digits” will be the number of units of two-digit numbers, i.e. 3,0,2,2,3,9,5, i.e. the terms of the mod 10 sequence.

Marsaglia is most famous for its

a set of tests of diehard generators

random numbers (RNG), so he understands this (here and below by RNG I mean a pseudo-random number generator (PRNG)). It made me wonder why it worked and how he picked six.

We will be writing in Raku, the language for gremlins. In case you’re a gremlin too, I’ll explain all the weird features below the spoilers.

Contents

## ▍ Introduction

The sequence is periodic, that is, if we iteratively apply it, sooner or later we will get the same element. Let’s start with the function (“subroutine”) that creates the entire sequence:

```
my sub next-rng(Int $start, Int $unitmult = 6, --> List) {
my @out;
my $next = $start;
repeat while $next !(elem) @out {
@out.append($next);
$next = sum($next.polymod(10) Z* $unitmult,1);
};
return @out;
}
```

**Explanation**

Raku is an extremely strange language, but I will try to explain as clearly as possible.

`@`

and`$`

are sigils for “positional” (resembling lists) and “scalar” variables. If you define a positional variable without giving it anything, it defaults to an empty array.`(elem)`

checks for membership, and`!`

can be used to negate any infix operator.`(elem)`

is the ASCII version, Raku can also accept ∈.- polymod divides a number into a remainder and a quotient, e.g.
`346.polymod(10) = (6 34)`

. This method can take multiple parameters, so you can perform operations like`num.polymod(60, 60, 24)`

to get hours-minutes-seconds. - Z is a “zip” metaoperator that applies an infix operator elementwise between two lists.
`(4, 6) Z* (6, 1)`

=`(4*6, 6*1)`

.

Once we have the sequence, we can output it using commands

`put`

or

`say`

whose behavior

is a little different

. I will not go into details.

`put next-rng(1);`

```
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
```

Remember that we need random numbers – this

*the last*

numbers That is, the RNG works like this: 1 -> 6 -> 6 -> 9 -> 7 -> 7 -> …

## ▍ We study the properties

If the RNG is even, each digit must occur in sequence the same number of times. We can check this by converting the last digits to a multiset, or “bag”.

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

**Explanation**

`<<op>>`

is a hyper operator that “maps” the inner operator on both lists,*also recursively entering a list of lists*. For example, the result`((1, 2), 3) <<+>> 10`

will be`((11, 12), 13)`

! Hyperoperators have many other strange properties that make them useful and confusing.- Objects of the Bag class count the number of elements in something.
`bag((4, 5, 4)){4} = 2`

. Oddly enough, they can only contain scalar variables, but not arrays, lists, or anything like that.

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

It looks like the distribution is fairly even, although the chance of getting a 0 or 9 is a little lower.

I got my next idea from the diehard tests. Quote from Wikipedia:

Overlapping PermutationsSequences of five random consecutive numbers are analyzed. The 120 possible permutations must occur with statistically equivalent probability.

There are only 54 five-digit sequences in the dataset, so I’ll apply the test to two-digit “transitions” instead. I will do this by displaying a 10 x 10 grid, where the (i, j)th index (from 0) is the number of transitions from

`i`

of the last digit to

`j`

. For example, a sequence contains a transition

`28 -> 50`

and no other species transitions

`X8 -> Y0`

therefore, there must be a value in cell (8, 0).

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

**Explanation**

`| @f, $x`

concatenates`$x`

directly in`@f`

. Without`|`

it would be a list of two elements.`*`

in``* % 10``

is whatever, a strange little operator that does a lot of things in different contexts, but usually in this case it turns an expression into a closure.*Usually*. It is the same as writing`map({$_ % 10})`

^{1}.- rotor
`(2 => -1)`

gets two elements, then goes back one element, then gets two more, and so on. As a result`[1, 2, 3, 4].rotor(2 => -1)`

will be`[(1, 2), (2, 3), (3, 4)]`

. You can also do it`rotor(2)`

To get`[(1, 2), (3, 4)]`

or`rotor(1 => 1)`

To get`[1, 3]`

. Rotor is a very cool thing. `^10`

– that’s all`0..9`

. Finally something simple!- X is the metaoperator of the vector product. That is, if
`$x = 2`

then`$x X ^4`

will give`((2 0), (2 1), (2 2), (2 3))`

. And yes, this operator can get a lot more expensive. `grep(foo)`

returns a list of all elements, performing a smart comparison`foo`

and`.elems`

– This is the number of items in the list. That is,`@pairs.grep($_).elems`

is the number of matching list elements`$_`

. It took me a while to figure it out*too*a lot of time

- Actually,
`-> $x {$x % 10}`

but it’s close enough

```
> 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
```

The table shows that some transitions are impossible. If I generated a 0, then immediately after it I will not be able to get a 6. Obviously, this is not a perfect RNG, but I did not expect something particularly good.

## ▍ Why 6?

What if I multiply the last number by 4 instead of 6?

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

Well, I don’t know at all, I like the RNG that never gives out 3. Specific sequences are called

*orbits*

And their length

*periods*

. Let’s look at all the possible orbits that can be obtained using 4 as a multiplier:

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

**Explanation**

`&`

is a sigil for “invoked object” or~~functions~~routines. Method`.assuming`

performs a partial application of the function, and the transfer`*`

causes it to be partially applied*second*parameter.^{1}`map`

returns the sequence of lists we pass`unique`

.`as => &set`

converts each sequence from`map`

in multitude and compares to uniqueness*their*, rather than the original lists. But in the end, the elements are used before the transformation.If this is not clear, then you can give a simpler example:

`[-1, 1].unique(as => &abs)`

returns`[-1]`

and`[1, -1].unique(as => &abs)`

returns`[1]`

.

- Daniel Sockwell (codesections) kindly agreed to read the first draft of this post and let me know about
`assume`

. Thanks to him!

```
> say orbits-for-mod(4, 38).map(*.gist).join("\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]
```

**Explanation**

I quote

Daniel Sockwell

:

`.map(*.gist).join("\n")`

used here only for good derivation.`cycles-for-mod`

returns`Seq`

with`Array`

; transformation of each`Array`

by using`.gist`

converts it to a string surrounded by square brackets, a`.join("\n")`

places a newline after each line.

If you choose 13 as the initial value, then the random numbers will be the sequence 3, 3, 3, 3, 3, 3.

For obvious reasons, 4 should never be a factor. In fact, for a multiplier to produce a “good” RNG, it must have one orbit.

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

**Explanation**

`.&`

implements a higher-order function as a method.`grep(*.&f)`

– This is equivalent`grep({f($_)})`

.`&orbits-for-mod`

returns a list.`==`

returns both input values of a number, and casting a list of numbers returns the number of elements. That is, we check whether the returned list contains only one element (in other words, only one orbit). (If you don’t want to do the comparison without pointing, use either === or eqv.)

This method is quite slow and only checks orbits, *are starting* from the number up to 20. Therefore, it will miss the orbit `26 -> 26`

for `x=4`

. I will fix both of these issues later.

Therefore, “good” options

`n`

– These are 6, 11 and 18.

If the result is a three-digit number, then the first two digits should be treated as one number. In case `n=11`

numeric `162`

leads to the meaning `16 + 22`

and not to `6 + 22`

(or to `6 + 1 + 22`

).

## ▍ Why does it work?

This part of the explanation was unclear to me:

Its period will be the order of the multiplier (6) of the group of residues prime to the module, 10 (in the case of 59).

After talking to friends and reading Wikipedia, I started to understand something. We calculate the RNG in our head

multiplication with transfer

(multiply with carry, MWC) with constants

`a=x`

and

`c=10`

. This choice has a convenient property: if

`MWC(x) = y`

then

`10y mod (10n-1) = x`

!

```
MWC: 01 -> 06 -> 36 -> ... -> 41 -> 10 -> 01
10y%59: 01 -> 10 -> 41 -> ... -> 36 -> 06 -> 01
```

And it is very wonderful! It’s easier for me to think about mathematically

`10y mod 59`

, than about “multiplying the last digit by six and adding the first digit.” For example, it is clear why the RNG generates 0 and 9 slightly less often than other numbers: whatever multiplier we choose, the sequence generated will go from 1 to 10n-2, “excluding” 10n-1 (which ends in 9) and 10n (which ends in 0).

“Multiplication and Division with Remainder” is also called Lemer’s RNG.

## ▍ Looking for better RNGs

And how do other numbers work? We know that a good multiplier only creates one orbit, and above I have shown the code to calculate this. Unfortunately, this algorithm has a complexity of O(n²) in the worst case.

*[Если множитель n имеет одну орбиту, тогда мы выполним next-rng примерно для 10n-2 чисел, и функция выполнит итерации 10n-2 раз (потому что ей нужно обойти каждое число в орбите). Если бы я сделал так, чтобы код пропускал числа, которые я уже видел в орбите, то время исполнения сократилось бы до O(n).]*

If we perceive the MWC algorithm as “Lemer’s reverse”, we can get a better way: if

`n`

is a good factor, then the period of the orbit starting with 1 should be

`10n-2`

.

Lemaire’s method also gives us a faster way to calculate the orbit:

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

**Explanation**

- Having indicated
`\x`

instead`$x`

as a parameter we can use in the body`x`

instead`$x`

. `...`

– This is a sequence operator. He can perform*many*different actions, but what is important for us is what if we write`10, &f … 1`

it will start with 10 and then apply`&f`

until it finally generates`1`

.- In
`* *10%[etc]`

: the first`*`

– is Whatever, and the second`*`

– Normal multiplication. This is then turned into a function`-> $a {$a * 10 % (10*x - 1)}`

.

This code actually creates the orbit in reverse, but we’re only interested in the period, so it doesn’t matter.

Then we check the period using the same trick with

`==`

that before

```
> say (1..100).grep({oneorbit($_) == 10*$_-2});
(2 3 6 11 15 18 23 27 38 39 42 50 51 62 66 71)
```

Now I understand why Marsaglia chose 6: most programmers know the multiplication table by 6 and understand that it will never return to a three-bit number, so the addition step is very simple. An orbit consists of only 58 numbers, and we won’t get sequences of the same numbers, but if you need to grab a few random numbers quickly, this is perfectly acceptable.

If you want more randomness, then I have a couple of candidates. 50 has a period of 498 and is incredibly easy to calculate. If the last digit is even, no carryover is required: **23**8 -> 4**23**!

Sequence for 50 though *looks* as random as other sequences. At a certain point, it generates 9 even numbers, followed by 8 odd numbers. Do not use it to simulate coin flips.

Another interesting number is 18. It has a decent period (178) and has all the possible number transitions:

```
> 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
```

Its disadvantage is that you will have to learn the multiplication table by 18. It is not very difficult: I mastered it in about ten minutes of practice. But I still am

*not perfect*

managing to run the entire MWC stage, but can consistently get the next random number every five seconds. I am quite satisfied with it.

The Raku code I used for research can be viewed here. It is configured as a CLI, but can be used in many different ways; see details in the file.

Discounts, raffle results and news about the RUVDS satellite — in our Telegram channel 🚀