Blum Blum Shub

10, Jul, 2015

Blum Blum Shub is a pseudo random number generator which has been proven to be secure when used appropriately. We’re not going to use it appropriately here.

What we are going to do is exploit the fact that you can compute any number in the sequence directly, which makes it suitable for procedural content generation. A hashing algorithm would probably be quicker, but BBS makes it easier to generate a long list of values without any repetition. Nobody wants two planets called Arse.

The formula for generating values sequentially is:

Xn+1 = Xn2 mod M

M is the product of two large primes (P * Q), where P mod 4 = 3, and Q mod 4 = 3. Here we will use a value of M which is representable in 32 bits, which is where cryptographers in the audience will start wincing. The reason for using a 32 bit value is that the square will fit in a 64 bit integer, and we don’t need to use a big number library.

The original paper recommends using special primes to produce a large cycle length, where P = 2 * P1 + 1, and P1 = 2 * P2 + 1. As I had a limited number of primes of an appropriate magnitude, I simply let the computer churn away at possible combinations. Within a couple of minutes it spat out 64007 and 66959, which give a cycle length of 535681477. Here, 64007 is a special prime but 66959 is not. Interestingly, I did not check for P and Q being equal, and a little while later it found that using 64007 for both gives an even longer cycle length (1024160005). This doesn’t fit the BBS criteria for P and Q to be distinct, but hey, it works!

The seed value (X0) isn’t too important, but an odd prime that isn’t P or Q usually works well. Computing a value directly uses Euler’s theorem:

Xi = (X02i mod λ(M)) mod M

Where λ(M) is the lowest common multiple of (P-1) and (Q-1). Here’s some code:

#include <stdio.h>
#include <stdint.h>  // For uint32_t / uint64_t
#include <algorithm> // For std::swap

uint32_t GreatestCommonDenominator(uint32_t x, uint32_t y)
{
	if (x < y)
		std::swap(x, y);

	while (y > 0)
	{
		x = x % y;
		std::swap(x, y);
	}
	return x;
}

uint32_t APowBModC(uint32_t a, uint32_t b, uint32_t c)
{
	if (b == 0)
		return 1;

	if (b == 1)
		return a % c;

	uint32_t b1 = b >> 1;
	uint32_t b2 = b - b1;
	uint64_t a_pow_b1_mod_c = APowBModC(a, b1, c);
	uint64_t a_pow_b2_mod_c;
	if (b1 == b2)
	{
		a_pow_b2_mod_c = a_pow_b1_mod_c;
	}
	else
	{
		a_pow_b2_mod_c = APowBModC(a, b2, c);
	}

	// a^b1 * a^b2 == a^(b1+b2)
	uint64_t result = (a_pow_b1_mod_c * a_pow_b2_mod_c) % c;
	return static_cast<uint32_t>(result);
}

class BlumBlumShub
{
	uint32_t P, Q, M, Seed, Value;

public:
	BlumBlumShub(uint32_t p, uint32_t q, uint32_t s)
	{
		P = p;
		Q = q;
		Seed = s;
		M = p*q;
		Value = s;
	}

	uint32_t NextInt()
	{
		uint64_t value64 = static_cast<uint64_t>(Value);
		Value = (value64*value64) % M;
		return Value;
	}

	uint32_t GetInt(uint32_t index)
	{
		uint32_t gcd = GreatestCommonDenominator(P - 1, Q - 1);
		uint32_t lcm = (P - 1)*(Q - 1) / gcd;
		uint32_t temp = APowBModC(2, index, lcm);
		return APowBModC(Seed, temp, M);
	}

	uint32_t RandMax() const
	{
		return M - 1;
	}
};

int main()
{
	BlumBlumShub bbs(64007, 66959, 48023);

	for (uint32_t i = 1; i <= 100; ++i)
	{
		printf("%u t %un", bbs.NextInt(), bbs.GetInt(i));
	}

	uint32_t Bitmap[32] = { 0 };
	for (uint32_t i = 0; i<100; ++i)
	{
		uint32_t x = bbs.NextInt() % 32;
		uint32_t y = bbs.NextInt() % 32;
		Bitmap[x] |= 1 << y;
	}

	for (uint32_t y = 0; y<32; ++y)
	{
		for (uint32_t x = 0; x<32; ++x)
		{
			printf("%c", Bitmap[x] & 1 << y ? '#' : ' ');
		}
		printf("n");
	}

	return 0;
}

Which will output the first 100 numbers calculated sequentially and directly, then a 32 x 32 grid with 100 random positions plotted.

# #            #   #   #      ##
                            #   
     #     #                    
               #     #          
                         #      
   #    # #    # #              
# #                        #    
                         #      
               # #  #           
#   ##   #      #               
  #            #              # 
            #              # #  
            #      #  # #    #  
      #               #     #   
     #     #   ##  # #          
                                
 #       #                      
                                
   # #                    #     
     #                       #  
 #     ##        #        #     
               #        ## #    
   #                            
                 #        # # # 
# #                  #          
 #                 #            
        #     #  ##       #     
         #   #  #    #          
       #                    #   
   ##         #                 
                  #   #         
##                        #     

References
A Simple Unpredictable Pseudo-Random Number Generator (Blum Blum Shub, 1986)