## Marsaglia’s Mental RNG

### January 29, 2019

We saw an “in-your-head” random number generator in a previous exercise, but I found it hard to operate. Today’s exercise, due to George Marsaglia is a random number generator that even a dummy like me can work in his head:

Choose a two-digit number as a seed. Then form a new two-digit number by adding the ten’s digit and six times the unit digit. For instance, if you start from 23, the sequence is 23 → 20 → 02 → 12 → 13 → 19 → 55 → 35 …. Then the sequence of random digits is the unit digit of each two-digit number. This operation produces the numbers 1 through 58 as seeds; the ten digits each occur six times, except that nine and zero occur five times (because 59 and 60 are not part of the sequence).

Your task is to implement Marsaglia’s mental RNG, and experiment with it to determine its effectiveness. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 5 Responses to “Marsaglia’s Mental RNG”

1. Ernie said

I found the last digits were not distributed randomly. For a sequence of 10000 numbers the digits 1…8 appeared 1035 +- 2 times while 0 and 9 appeared ~862 times

2. Alex B said

Here’s a Python version as a generator:

```def marsaglia_generator(seed):
"""Generate a pseudorandom sequence based on Marsaglia’s Mental RNG"""

if seed < 10 or seed > 99:
raise ValueError('Seed value must be two digits.')
if seed == 59:
raise ValueError('Seed value 59 will cause an infinite loop.')

def next_seed(seed):
return seed // 10 + seed % 10 * 6
while seed > 58:
seed = next_seed(seed)
new_seed = seed
while True:
yield new_seed % 10
new_seed = next_seed(new_seed)
# chain is exhausted when original number comes round again
if new_seed == seed:
break
```

I’m not sure it’s a particularly good PRNG, as each seed simply produces the same sequence in rotation, but it certainly is simple to use.
Note: DO NOT use 59 as the seed :)

3. matthew said

Well, it’s a pretty good RNG considering it’s only got 6 bits of state! If the shortage of 0 and 9 is bothersome, you can always skip values greater than 50:

```def rng(n):
def f():
nonlocal n
while True:
n = n//10+n%10*6
if n <= 50: return n
return f

print([f() for f in [rng(1)] for _ in range(50)])
```
```[6, 36, 39, 47, 46, 40, 4, 24, 26, 38,
11, 7, 42, 16, 37, 45, 34, 27, 44, 28,
50, 5, 30, 3, 18, 49, 23, 20, 2, 12,
13, 19, 35, 33, 21, 8, 48, 17, 43, 22,
14, 25, 32, 15, 31, 9, 29, 41, 10, 1]
```
4. matthew said

It’s interesting to generalize a bit and work in base b (this is a lag-1 Multiply with Carry RNG, as Praxis says). The inverse operation to bx + y -> ay + x is just n -> bn mod (ab-1) ie. the sequence, starting and ending at 1 is just the reverse of the sequence 1,b,b^2,b^3,… taking each element mod m = ab-1 so the period is the period of b, mod m. For a given b, this is longest when m is a prime and b is a primitive root mod m, and we do particularly well if we can take a = b-1 since the sequence then covers nearly all of the range [1..b*b]. For example, for b = 127, a = 126, we get a sequence of length 16000 including all but 129 of the 16129 two digit base-127 numbers. Here’s a program that calculates such maximal sequences:

```def f(b,a):
s = 1
while True:
yield(s)
s1 = s//b + s%b*a
assert(s == s1*b%(a*b-1))
s = s1
if s == 1: break

def test():
for b in range(3,1000):
a = b-1
m = a*b-1
s = list(f(b,a))
if len(s) == m-1:
print(b,a,m-1,b*b)
test()
```

Like 127, 999 is interesting (this is a CMWC generator, I believe):

```3 2 4 9
7 6 40 49
11 10 108 121
27 26 700 729
31 30 928 961
39 38 1480 1521
51 50 2548 2601
67 66 4420 4489
87 86 7480 7569
115 114 13108 13225
127 126 16000 16129
131 130 17028 17161
135 134 18088 18225
...
967 966 934120 935089
975 974 949648 950625
995 994 989028 990025
999 998 997000 998001
```

I don’t think that sequence is in the OEIS.

5. Daniel said

Here’s a solution in C.

```#include <stdio.h>
#include <stdlib.h>

int main(int argc, char* argv[]) {
if (argc != 3) {
fprintf(stderr, "Usage: %s <seed> <count>\n", argv[0]);
return EXIT_FAILURE;
}
int x = atoi(argv[1]);
int count = atoi(argv[2]);
while (count-- > 0) {
printf("%d\n", x % 10);
x = (x / 10) + (x % 10) * 6;
}
return EXIT_SUCCESS;
}
```

Example usage:

\$ ./a.out 32 10
2
5
1
9
4
9
6
1
0
1