June 10, 2022 - Reading time: 5 minutes
Blum Blum Shub is a PRNG algorithm published in 1986.
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² 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.
3 (mod 4). This means that p % 4 and q % 4 both need to be 3.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 [1]:
def blum_blum_shub(x, m):
return (x * x) % m
These parameters are from Wikipedia:
In [2]:
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 [3]:
x = seed
for _ in range(21):
x = blum_blum_shub(x, M)
print(x, end=" ")
Out [1]:
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 [4]:
x = seed
for i in range(21):
x = blum_blum_shub(x, M)
bit = x & 1
print(bit, end="")
Out [2]:
110010100000110110011
Let’s encapsulate all of this in a simple Python class.
In [5]:
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 [6]:
bbs = BlumBlumShub(seed, M)
for _ in range(21):
print(bbs.next_state(), end=" ")
Out [3]:
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 [7]:
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 [8]:
bbs = BlumBlumShub(seed, M)
bbs.buffer(64).hex()
Out [4]:
'ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0d9ca0'
As we can see, the stream starts repeating very quickly. To mitigate this, we need better values for p, q and seed.
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 [9]:
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 [10]:
get_suitable_prime(urandom, 128)
Out [5]:
398695105695143609311998375660109791359
In [11]:
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 [12]:
p = get_suitable_prime(urandom, 128)
q = get_suitable_prime(urandom, 128)
pick_seed(p, q, urandom, 128)
Out [6]:
270385993931251981005558974782161044015
Putting it all together.
In [13]:
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 [14]:
get_parameters(urandom, 32)
Out [7]:
Parameters(p=6273284939, q=1221558803, m=7663186440962768017, seed=2819474525)
In [15]:
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 [16]:
rng = keyed_rng(b"secret key")
bytes([rng() for _ in range(32)]).hex()
Out [8]:
'3e0b2f19af8ddb7b93ce65e1fd18e1027e662088fd7a5c6beb3861e9f42891a5'
In [17]:
get_parameters(keyed_rng(b"hello"), 16)
get_parameters(keyed_rng(b"world"), 16)
get_parameters(keyed_rng(b"hello"), 16)
Out [9]:
Parameters(p=64319, q=39983, m=2571666577, seed=26252)
Parameters(p=22343, q=48563, m=1085043109, seed=3393)
Parameters(p=64319, q=39983, m=2571666577, seed=26252)
In [18]:
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 [19]:
plaintext = b"Hello, world! This is Blum Blum Shub."
ciphertext = encrypt(b"test key", plaintext)
ciphertext.hex()
Out [10]:
'ed17127e9217d04733e1af4de2cd27a6fdf2b5eb4c1ba529302743e31d4f064ad4f4c329df'
Decryption is the same as encryption.
In [20]:
decrypt = encrypt
In [21]:
decrypt(b"test key", ciphertext)
Out [11]:
b'Hello, world! This is Blum Blum Shub.'
Let’s try to decrypt with the wrong password.
In [22]:
decrypt(b"Test Key", ciphertext)
Out [12]:
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'
In [23]:
Out [13]: