Marsaglia’s Mental RNG

January 29, 2019

Since this is a mental exercise, there is really no program to write. Here is the sequence of seeds, starting from 23:

> (define rs
    (do ((rs (list 23)
          (cons (+ (quotient (car rs) 10)
                   (* (modulo (car rs) 10) 6))
                rs)))
        ((member (car rs) (cdr rs))
          (reverse (cdr rs)))))
> rs
(23 20 2 12 13 19 55 35 33 21 8 48 52 17 43 22 14 25 32 15
 31 9 54 29 56 41 10 1 6 36 39 57 47 46 40 4 24 26 38 51 11 7
 42 16 37 45 34 27 44 28 50 5 30 3 18 49 58 53)

Here we compute the random digits by mapping the seed sequence mod 10, and compute a potentially-winning sequence for your next high-stakes rock-paper-scissors match:

> (define (mod n) (lambda (x) (modulo x n)))
> (map (mod 10) rs)
(3 0 2 2 3 9 5 5 3 1 8 8 2 7 3 2 4 5 2 5 1 9 4 9 6 1 0 1 6 6
 9 7 7 6 0 4 4 6 8 1 1 7 2 6 7 5 4 7 4 8 0 5 0 3 8 9 8 3)
> (map (lambda (n)
         (case n ((0) 'rock) ((1) 'paper) ((2) 'scissors)))
       (map (mod 3) rs))
(scissors scissors scissors rock paper paper paper scissors
 rock rock scissors rock paper scissors paper paper scissors
 paper scissors rock paper rock rock scissors scissors
 scissors paper paper rock rock rock rock scissors paper
 paper paper rock scissors scissors rock scissors paper rock
 paper paper rock paper rock scissors paper scissors scissors
 rock rock rock paper paper scissors)

As expected, there are nearly-equal counts of all digits:

> (uniq-c = (sort < (map (mod 10) rs)))
((0 . 5) (1 . 6) (2 . 6) (3 . 6) (4 . 6) (5 . 6) (6 . 6)
  (7 . 6) (8 . 6) (9 . 5))

But looking at adjacent pairs of digits show some bias:

> (define ps
    (let loop ((rs (append (map (mod 10) rs)
                           (list ((mod 10) (car rs)))))
               (ps (list)))
      (if (null? (cdr rs)) (reverse ps)
        (loop (cdr rs) (cons (cons (car rs) (cadr rs)) ps)))))
> ps
((3 . 0) (0 . 2) (2 . 2) (2 . 3) (3 . 9) (9 . 5) (5 . 5)
 (5 . 3) (3 . 1) (1 . 8) (8 . 8) (8 . 2) (2 . 7) (7 . 3)
 (3 . 2) (2 . 4) (4 . 5) (5 . 2) (2 . 5) (5 . 1) (1 . 9)
 (9 . 4) (4 . 9) (9 . 6) (6 . 1) (1 . 0) (0 . 1) (1 . 6)
 (6 . 6) (6 . 9) (9 . 7) (7 . 7) (7 . 6) (6 . 0) (0 . 4)
 (4 . 4) (4 . 6) (6 . 8) (8 . 1) (1 . 1) (1 . 7) (7 . 2)
 (2 . 6) (6 . 7) (7 . 5) (5 . 4) (4 . 7) (7 . 4) (4 . 8)
 (8 . 0) (0 . 5) (5 . 0) (0 . 3) (3 . 8) (8 . 9) (9 . 8)
 (8 . 3) (3 . 3))
> (sort (lambda (a b) (or (< (car a) (car b))
                          (and (= (car a) (car b))
                               (< (cdr a) (cdr b)))))
        ps)
((0 . 1) (0 . 2) (0 . 3) (0 . 4) (0 . 5) (1 . 0) (1 . 1)
 (1 . 6) (1 . 7) (1 . 8) (1 . 9) (2 . 2) (2 . 3) (2 . 4)
 (2 . 5) (2 . 6) (2 . 7) (3 . 0) (3 . 1) (3 . 2) (3 . 3)
 (3 . 8) (3 . 9) (4 . 4) (4 . 5) (4 . 6) (4 . 7) (4 . 8)
 (4 . 9) (5 . 0) (5 . 1) (5 . 2) (5 . 3) (5 . 4) (5 . 5)
 (6 . 0) (6 . 1) (6 . 6) (6 . 7) (6 . 8) (6 . 9) (7 . 2)
 (7 . 3) (7 . 4) (7 . 5) (7 . 6) (7 . 7) (8 . 0) (8 . 1)
 (8 . 2) (8 . 3) (8 . 8) (8 . 9) (9 . 4) (9 . 5) (9 . 6)
 (9 . 7) (9 . 8))

Rearranging that so it makes a little bit more sense, we see:

0: 1 2 3 4 5
1: 6 7 8 9 0 1
2: 2 3 4 5 6 7
3: 8 9 0 1 2 3
4: 4 5 6 7 8 9
5: 0 1 2 3 4 5
6: 6 7 8 9 0 1
7: 2 3 4 5 6 7
8: 8 9 0 1 2 3
9: 4 5 6 7 8

The bias is evident in the ascending sequences of digits, some of which wrap around, although the ascending sequences aren’t ordered in the random sequence. Of course, there must be bias, since the calculations are so small there isn’t time for the generator to overcome the bias. But it’s still a worthy way to generate random digits in your head. If you want to know more, Marsaglia’s method is an instance of his multiply-with-carry random number generator.

You can run the program at https://ideone.com/Y9pNuZ.

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: