Blum Blum Shub is a PRNG algorithm published in 1986.

Formula

The algorithm is very short and simple. Starting from the seed, the next state can be computed by passing the current state through the following formula.

f(x) = x2 mod M

In this formula, M is the product of p and q, two large primes.

The complexity in this algorithm is hidden in the parameters; the seed and the modulus M. In order to have a long cycle length and fulfill its security promises, Blum Blum Shub has a few constraints on its parameters.

In contrast, some more complex PRNG algorithms can work with pretty much any randomized seed.

Constraints

  • The seed should be co-prime to p and q. This means their greatest common divisor should be 1.
  • p and q need to be congruent to 3 (mod 4). This means that p % 4 and q % 4 both need to be 3.
  • p and q should be safe primes.

Exploration

Before implementing the algorithm properly, I am going to play around with example values and see how the function behaves. Afterwards, I’ll encapsulate our discoveries here into more useful methods.

Let’s start with a one-to-one translation of the formula from before.

In [4]:
def blum_blum_shub(x, m):
    return (x * x) % m

These parameters are from Wikipedia.

In [5]:
P = 11
Q = 23
M = P * Q

seed = 3

The numbers are very small, resulting in a very short cycle length of 20 elements.

In [6]:
x = seed

for _ in range(21):
    x = blum_blum_shub(x, M)
    print(x, end=" ")
Out:
9 81 236 36 31 202 71 234 108 26 170 58 75 59 192 179 163 4 16 3 9

Instead of using the full state, we will be taking one bit from each iteration.

In [7]:
x = seed

for i in range(21):
    x = blum_blum_shub(x, M)
    bit = x & 1
    print(bit, end="")
Out:
110010100000110110011

Python implementation

Let’s encapsulate all of this in a simple Python class.

In [8]:
class BlumBlumShub:
    def __init__(self, seed, mod):
        self.x = seed
        self.mod = mod
    
    def next_state(self):
        self.x = powmod(self.x, 2, self.mod)
        return self.x

Let’s see if we get the same output.

In [9]:
bbs = BlumBlumShub(seed, M)

for _ in range(21):
    print(bbs.next_state(), end=" ")
Out:
9 81 236 36 31 202 71 234 108 26 170 58 75 59 192 179 163 4 16 3 9

Looks the same. We can now add our helpers to generate bits and bytes from this number stream.

In [10]:
class BlumBlumShub(BlumBlumShub):
    def next_bit(self):
        return self.next_state() & 1
    
    def next_byte(self):
        byte = 0
        
        for _ in range(8):
            byte <<= 1
            byte |= self.next_bit()
        
        return byte
    
    def buffer(self, size):
        buf = bytearray(size)
        
        for i in range(size):
            buf[i] = self.next_byte()
        
        return buf
In [11]:
bbs = BlumBlumShub(seed, M)

bbs.buffer(64).hex()
Out [11]:
'ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0'

As we can see, the stream starts repeating very quickly. To mitigate this, we need better values for p, q and seed.

Parameter selection

An RNG is not very useful unless we can generate different streams of numbers. In order to do that, we need to generate the parameters M and seed.

For some PRNG algorithms, you can pick these as uniform random values. But Blum Blum Shub has extra requirements as discussed in the Constraints section.

In [12]:
def random_int(rng, bits):
    size = int(bits / 8)
    buf = [rng() for _ in range(size)]
    buf = bytes(buf)
    return int.from_bytes(buf, 'big')

def random_prime(rng, bits):
    n = random_int(rng, bits)
    n = next_prime(n)
    return int(n)

def get_safe_prime(rng, bits):
    while True:
        n = random_prime(rng, bits)
        n = 2 * n + 1
        if is_prime(n):
            return n

def get_suitable_prime(rng, bits):
    while True:
        n = get_safe_prime(rng, bits)
        if n % 4 == 3:
            return n
In [13]:
get_suitable_prime(urandom, 128)
Out [13]:
398695105695143609311998375660109791359

Picking a seed

In [14]:
def pick_seed(p, q, rng, bits):
    while True:
        n = random_int(rng, bits)
        
        if n == 0 or n == 1:
            continue
        
        if gcd(n, p) == 1 and gcd(n, q) == 1:
            return n
In [15]:
p = get_suitable_prime(urandom, 128)
q = get_suitable_prime(urandom, 128)
pick_seed(p, q, urandom, 128)
Out [15]:
270385993931251981005558974782161044015

Putting it all together.

In [16]:
Parameters = namedtuple("Parameters", "p q m seed")

def get_parameters(rng, bits):
    p = get_suitable_prime(rng, bits)
    q = get_suitable_prime(rng, bits)
    m = p * q
    
    seed = pick_seed(p, q, rng, bits)
    return Parameters(p, q, m, seed)
In [17]:
get_parameters(urandom, 32)
Out [17]:
Parameters(p=6273284939, q=1221558803, m=7663186440962768017, seed=2819474525)

Keyed selection

Key derivation

In [18]:
def keyed_rng(key):
    i = 0
    
    def inner():
        nonlocal i
        buf = key + i.to_bytes(3, 'big')
        h = sha256(buf)
        i += 1
        return h[0]
    
    return inner
In [19]:
rng = keyed_rng(b"secret key")

bytes([rng() for _ in range(32)]).hex()
Out [19]:
'3e0b2f19af8ddb7b93ce65e1fd18e1027e662088fd7a5c6beb3861e9f42891a5'
In [20]:
get_parameters(keyed_rng(b"hello"), 16)
get_parameters(keyed_rng(b"world"), 16)
get_parameters(keyed_rng(b"hello"), 16)
Out [20]:
Parameters(p=64319, q=39983, m=2571666577, seed=26252)
Out [20]:
Parameters(p=22343, q=48563, m=1085043109, seed=3393)
Out [20]:
Parameters(p=64319, q=39983, m=2571666577, seed=26252)

Usage as a cipher

Encryption

In [21]:
def encrypt(key, data):
    rng = keyed_rng(key)
    params = get_parameters(rng, 256)
    bbs = BlumBlumShub(params.seed, params.m)
    
    res = bytearray(len(data))
    
    for i, c in enumerate(data):
        res[i] = c ^ bbs.next_byte()
    
    return bytes(res)
In [22]:
plaintext = b"Hello, world! This is Blum Blum Shub."
ciphertext = encrypt(b"test key", plaintext)

ciphertext.hex()
Out [22]:
'ed17127e9217d04733e1af4de2cd27a6fdf2b5eb4c1ba529302743e31d4f064ad4f4c329df'

Decryption

Decryption is the same as encryption.

In [23]:
decrypt = encrypt
In [24]:
decrypt(b"test key", ciphertext)
Out [24]:
b'Hello, world! This is Blum Blum Shub.'

Let’s try to decrypt with the wrong password.

In [25]:
decrypt(b"Test Key", ciphertext)
Out [25]:
b'\xd0\xc9+\xdd\xefL\xe3$\xd3\xf6\xd8\xda\\\x85\x81\x07/\xd3\xfclX\xebo\x8c\xc8\xff\xd5\x0f\x0b\xbe\xc7\xa7\xbf\xac\xc1\xfa\xee'

References and useful links

  • Blum, Blum, and Shub, “A simple unpredictable pseudo-random number generator”, May 1986. PDF
  • Blum Blum Shub on Wikipedia