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 [2]:
def blum_blum_shub(x, m):
    return (x * x) % m

These parameters are from Wikipedia.

In [3]:
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 [4]:
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 [5]:
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 [6]:
class BlumBlumShub:
    def __init__(self, seed, mod):
        self.x = seed
        self.mod = mod

In order to calculate the next state, we are still going to use the same formula. But because we’re going to be working with large numbers, it’s not very efficient to square the number and then take the modulo as two separate steps.

There are efficient algorithms for computing modular exponentiation, which we can use instead. These algorithms are usually called powmod or modpow in programming languages. Python helpfully provides a built-in function for this called pow.

If you call pow with three arguments, it will compute the modular exponentiation of the first two arguments, using the third argument as the modulus.

In [7]:
class BlumBlumShub(BlumBlumShub):
    def next_state(self):
        self.x = pow(self.x, 2, self.mod)
        return self.x

Let’s see if we get the same output.

In [8]:
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. This is pretty simple.

  • To generate a bit, we take one bit from the state.
  • To generate a byte, we generate 8 bits.
  • To generate a buffer, we generate a bunch of bytes.
In [9]:
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 bytes(buf)
In [10]:
bbs = BlumBlumShub(seed, M)

bbs.buffer(64).hex()
Out [10]:
'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 [11]:
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):
    while True:
        n = random_int(rng, bits) | 1 # Random odd number
        if is_prime(n): return 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 the section above, we tested numbers for their primality with an is_prime function. This function uses the Miller-Rabin primality test to check if a number is (probably) prime.

If you want to know more about this, I have a notebook on the Miller-Rabin primality test.

In [12]:
urandom = lambda: os.urandom(1)[0]

get_suitable_prime(urandom, 128)
Out [12]:
438875015154062593704416815007831812283

Picking a seed

Continuing from the Constraints section, we need to pick a seed that is co-prime to p and q. This means that the greatest common divisor of the seed and p and q should be 1.

In [13]:
def pick_seed(p, q, rng, bits):
    while True:
        n = random_int(rng, bits)

        if n == 0 or n == 1:
            continue

        if math.gcd(n, p) == 1 and math.gcd(n, q) == 1:
            return n
In [14]:
p = get_suitable_prime(urandom, 128)
q = get_suitable_prime(urandom, 128)
pick_seed(p, q, urandom, 128)
Out [14]:
51852778346964364772912943704014255005

Putting it all together.

In [15]:
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 [16]:
get_parameters(urandom, 32)
Out [16]:
Parameters(p=3263052707, q=5847777359, m=19081605741218260813, seed=1033830034)

Keyed selection

Key derivation

In [17]:
class KeyedRNG:
    def __init__(self, key):
        self.key = key
        self.i = 0

    def __call__(self):
        buf = self.key + self.i.to_bytes(3, 'big')
        h = hashlib.sha256(buf).digest()
        self.i += 1
        return h[0]
In [18]:
rng = KeyedRNG(b"secret key")

bytes([rng() for _ in range(32)]).hex()
Out [18]:
'3e0b2f19af8ddb7b93ce65e1fd18e1027e662088fd7a5c6beb3861e9f42891a5'
In [19]:
get_parameters(KeyedRNG(b"hello"), 16)
get_parameters(KeyedRNG(b"world"), 16)
get_parameters(KeyedRNG(b"hello"), 16)
Out [19]:
Parameters(p=66107, q=58907, m=3894165049, seed=16044)
Out [19]:
Parameters(p=14159, q=25583, m=362229697, seed=28657)
Out [19]:
Parameters(p=66107, q=58907, m=3894165049, seed=16044)

Usage as a cipher

Encryption

In [20]:
def encrypt(key, data):
    rng = KeyedRNG(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 [21]:
plaintext = b"Hello, world! This is Blum Blum Shub."
ciphertext = encrypt(b"test key", plaintext)

ciphertext.hex()
Out [21]:
'082ae382e7478ada5bfb11985a54a454cda158100c9e69208ed5a61f6aed1fb79038630a6d'

Decryption

Decryption is the same as encryption.

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

Let’s try to decrypt with the wrong password.

In [24]:
decrypt(b"Test Key", ciphertext)
Out [24]:
b'7O\xbc"\xa8\x12\x96K("\x8eq\t\xcc2\x0e\xc3_\xd3\xd7G\xff\x11\x9f\xad\x17LD\x94l\x1aK\xb8\xa5\x08k\x04'

References and useful links

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