This article has multiple issues. Please help improve it or discuss these issues on the talk page. (Learn how and when to remove these messages) This article includes a list of general references, but it lacks sufficient corresponding inline citations. Please help to improve this article by introducing more precise citations. (September 2013) (Learn how and when to remove this message) This article relies excessively on references to primary sources. Please improve this article by adding secondary or tertiary sources. Find sources: "Blum Blum Shub" – news · newspapers · books · scholar · JSTOR (September 2013) (Learn how and when to remove this message) (Learn how and when to remove this message)

Blum Blum Shub (B.B.S.) is a pseudorandom number generator proposed in 1986 by Lenore Blum, Manuel Blum and Michael Shub^[1] that is derived from Michael O. Rabin's one-way function.

Blum Blum Shub takes the form

${\displaystyle x_{n+1}=x_{n}^{2}{\bmod {M}}}$,

where M = pq is the product of two large primes p and q. At each step of the algorithm, some output is derived from x_n+1; the output is commonly either the bit parity of x_n+1 or one or more of the least significant bits of x_n+1.

The seed x₀ should be an integer that is co-prime to M (i.e. p and q are not factors of x₀) and not 1 or 0.

The two primes, p and q, should both be congruent to 3 (mod 4) (this guarantees that each quadratic residue has one square root which is also a quadratic residue), and should be safe primes with a small gcd((p-3)/2, (q-3)/2) (this makes the cycle length large).

An interesting characteristic of the Blum Blum Shub generator is the possibility to calculate any x_i value directly (via Euler's theorem):

${\displaystyle x_{i}=\left(x_{0}^{2^{i}{\bmod {\lambda }}(M)}\right){\bmod {M}}}$,

where ${\displaystyle \lambda }$ is the Carmichael function. (Here we have ${\displaystyle \lambda (M)=\lambda (p\cdot q)=\operatorname {lcm} (p-1,q-1)}$).

There is a proof reducing its security to the computational difficulty of factoring.^[1] When the primes are chosen appropriately, and O(log log M) lower-order bits of each x_n are output, then in the limit as M grows large, distinguishing the output bits from random should be at least as difficult as solving the quadratic residuosity problem modulo M.

The performance of the BBS random-number generator depends on the size of the modulus M and the number of bits per iteration j. While lowering M or increasing j makes the algorithm faster, doing so also reduces the security. A 2005 paper gives concrete, as opposed to asymptotic, security proof of BBS, for a given M and j. The result can also be used to guide choices of the two numbers by balancing expected security against computational cost.^[2]

Let ${\displaystyle p=11}$, ${\displaystyle q=23}$ and ${\displaystyle s=3}$ (where ${\displaystyle s}$ is the seed). We can expect to get a large cycle length for those small numbers, because ${\displaystyle {\rm {gcd}}((p-3)/2,(q-3)/2)=2}$. The generator starts to evaluate ${\displaystyle x_{0}}$ by using ${\displaystyle x_{-1}=s}$ and creates the sequence ${\displaystyle x_{0}}$, ${\displaystyle x_{1}}$, ${\displaystyle x_{2}}$, ${\displaystyle \ldots }$ ${\displaystyle x_{5}}$ = 9, 81, 236, 36, 31, 202. The following table shows the output (in bits) for the different bit selection methods used to determine the output.

The following is a Python implementation that does check for primality.

import sympy
def blum_blum_shub(p1, p2, seed, iterations):
  assert p1 % 4 == 3
  assert p2 % 4 == 3
  assert sympy.isprime(p1//2)
  assert sympy.isprime(p2//2)
  n = p1 * p2
  numbers = []
  for _ in range(iterations):
    seed = (seed ** 2) % n
    if seed in numbers:
      print(f"The RNG has fallen into a loop at {len(numbers)} steps")
      return numbers
    numbers.append(seed)
  return numbers

print(blum_blum_shub(11, 23, 3, 100))

The following Common Lisp implementation provides a simple demonstration of the generator, in particular regarding the three bit selection methods. It is important to note that the requirements imposed upon the parameters p, q and s (seed) are not checked.

(defun get-number-of-1-bits (bits)
  "Returns the number of 1-valued bits in the integer-encoded BITS."
  (declare (type (integer 0 *) bits))
  (the (integer 0 *) (logcount bits)))

(defun get-even-parity-bit (bits)
  "Returns the even parity bit of the integer-encoded BITS."
  (declare (type (integer 0 *) bits))
  (the bit (mod (get-number-of-1-bits bits) 2)))

(defun get-least-significant-bit (bits)
  "Returns the least significant bit of the integer-encoded BITS."
  (declare (type (integer 0 *) bits))
  (the bit (ldb (byte 1 0) bits)))

(defun make-blum-blum-shub (&key (p 11) (q 23) (s 3))
  "Returns a function of no arguments which represents a simple
   Blum-Blum-Shub pseudorandom number generator, configured to use the
   generator parameters P, Q, and S (seed), and returning three values:
   (1) the number x[n+1],
   (2) the even parity bit of the number,
   (3) the least significant bit of the number.
   ---
   Please note that the parameters P, Q, and S are not checked in
   accordance to the conditions described in the article."
  (declare (type (integer 0 *) p q s))
  (let ((M    (* p q))       ;; M  = p * q
        (x[n] s))            ;; x0 = seed
    (declare (type (integer 0 *) M x[n]))
    #'(lambda ()
        ;; x[n+1] = x[n]^2 mod M
        (let ((x[n+1] (mod (* x[n] x[n]) M)))
          (declare (type (integer 0 *) x[n+1]))
          ;; Compute the random bit(s) based on x[n+1].
          (let ((even-parity-bit       (get-even-parity-bit       x[n+1]))
                (least-significant-bit (get-least-significant-bit x[n+1])))
            (declare (type bit even-parity-bit))
            (declare (type bit least-significant-bit))
            ;; Update the state such that x[n+1] becomes the new x[n].
            (setf x[n] x[n+1])
            (values x[n+1]
                    even-parity-bit
                    least-significant-bit))))))

;; Print the exemplary outputs.
(let ((bbs (make-blum-blum-shub :p 11 :q 23 :s 3)))
  (declare (type (function () (values (integer 0 *) bit bit)) bbs))
  (format T "~&Keys: E = even parity, L = least significant")
  (format T "~2%")
  (format T "~&x[n+1] | E | L")
  (format T "~&--------------")
  (loop repeat 6 do
    (multiple-value-bind (x[n+1] even-parity-bit least-significant-bit)
        (funcall bbs)
      (declare (type (integer 0 *) x[n+1]))
      (declare (type bit           even-parity-bit))
      (declare (type bit           least-significant-bit))
      (format T "~&~6d | ~d | ~d"
                x[n+1] even-parity-bit least-significant-bit))))

• GMPBBS, a C-language implementation by Mark Rossmiller
• BlumBlumShub, a Java-language implementation by Mark Rossmiller
• An implementation in Java
• Randomness tests