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) = x

^{2}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.

def blum_blum_shub(x, m): return (x * x) % m

These parameters are from Wikipedia.

P = 11 Q = 23 M = P * Q seed = 3

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

x = seed for _ in range(21): x = blum_blum_shub(x, M) print(x, end=" ")

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.

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

110010100000110110011

# Python implementation

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

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.

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.

bbs = BlumBlumShub(seed, M) for _ in range(21): print(bbs.next_state(), end=" ")

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.

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)

bbs = BlumBlumShub(seed, M) bbs.buffer(64).hex()

'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.

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.

urandom = lambda: os.urandom(1)[0] get_suitable_prime(urandom, 128)

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.

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

p = get_suitable_prime(urandom, 128) q = get_suitable_prime(urandom, 128) pick_seed(p, q, urandom, 128)

51852778346964364772912943704014255005

Putting it all together.

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)

get_parameters(urandom, 32)

Parameters(p=3263052707, q=5847777359, m=19081605741218260813, seed=1033830034)

## Keyed selection

### Key derivation

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]

rng = KeyedRNG(b"secret key") bytes([rng() for _ in range(32)]).hex()

'3e0b2f19af8ddb7b93ce65e1fd18e1027e662088fd7a5c6beb3861e9f42891a5'

get_parameters(KeyedRNG(b"hello"), 16) get_parameters(KeyedRNG(b"world"), 16) get_parameters(KeyedRNG(b"hello"), 16)

Parameters(p=66107, q=58907, m=3894165049, seed=16044)

Parameters(p=14159, q=25583, m=362229697, seed=28657)

Parameters(p=66107, q=58907, m=3894165049, seed=16044)

# Usage as a cipher

## Encryption

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)

plaintext = b"Hello, world! This is Blum Blum Shub." ciphertext = encrypt(b"test key", plaintext) ciphertext.hex()

'082ae382e7478ada5bfb11985a54a454cda158100c9e69208ed5a61f6aed1fb79038630a6d'

## Decryption

Decryption is the same as encryption.

decrypt = encrypt

decrypt(b"test key", ciphertext)

b'Hello, world! This is Blum Blum Shub.'

Let’s try to decrypt with the wrong password.

decrypt(b"Test Key", ciphertext)

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