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.

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

These parameters are from Wikipedia.

In :
```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 :
```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 :
```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 :
```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 :
```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 :
```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 :
```bbs = BlumBlumShub(seed, M)

bbs.buffer(64).hex()```
Out :
`'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 :
```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 :
`get_suitable_prime(urandom, 128)`
Out :
`398695105695143609311998375660109791359`

## Picking a seed

In :
```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 :
```p = get_suitable_prime(urandom, 128)
q = get_suitable_prime(urandom, 128)
pick_seed(p, q, urandom, 128)```
Out :
`270385993931251981005558974782161044015`

Putting it all together.

In :
```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 :
`get_parameters(urandom, 32)`
Out :
`Parameters(p=6273284939, q=1221558803, m=7663186440962768017, seed=2819474525)`

## Keyed selection

### Key derivation

In :
```def keyed_rng(key):
i = 0

def inner():
nonlocal i
buf = key + i.to_bytes(3, 'big')
h = sha256(buf)
i += 1
return h

return inner```
In :
```rng = keyed_rng(b"secret key")

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

# Usage as a cipher

## Encryption

In :
```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 :
```plaintext = b"Hello, world! This is Blum Blum Shub."
ciphertext = encrypt(b"test key", plaintext)

ciphertext.hex()```
Out :
`'ed17127e9217d04733e1af4de2cd27a6fdf2b5eb4c1ba529302743e31d4f064ad4f4c329df'`

## Decryption

Decryption is the same as encryption.

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

Let’s try to decrypt with the wrong password.

In :
`decrypt(b"Test Key", ciphertext)`
Out :
`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'`