(*******************************************************) 
 (* 
      The BBS - 
        Generator in Mathematica                   *) 
 (*                                                    *) 
 (* 
      Author :  
        Klaus Pommerening                          *) 
 (*        Johannes - Gutenberg - 
        Universität              *) 
 (*        Mainz,  
      Germany                              *) 
 (* 
      27.  November 1997                                  *) 
 (* 
      last change :  
          28.  July 2003                         *)  (*******************************************************)

(*******************************************************) 
 (* 
      Generate smallest prime p ≥ 
        x such that (p - 1)/
            2 is    *) 
 (*   prime, 
      in particular p is congruent to 3 modulo 4.  *) 
 (*   (This is called a special BLUM prime.)            *) 
 (* WARNING :  
        Since PrimeQ uses the probabilistic RABIN  *) 
 (*   test,  it is not absolutely certain,  
      that the      *) 
 (*   result                                            *) 
 (*   is actually prime.                                *) 
 (*******************************************************) 
BlumPrime[
        x_Integer] := 
Module[{p = x}, 
        If[EvenQ[p], p++]; 
If[EvenQ[(p - 1)/2], 
          p += 2]; 
While[! PrimeQ[(p - 1)/
                  2] || ! PrimeQ[p], p += 4]; 
Return[
          p]
]

(*******************************************************) 
 (* 
      Generate a random special BLUM prime with n Bits.   *) 
 (* 
      REMARK :  
        The probabilitiy that a BLUM prime p occurs *) 
 (*   as output is proportional to the difference of p  *) 
 (*   and the largest Blum prime  < p .  That is,  
      the       *) 
 (*   distribution is not uniform.                       *) 
 (*******************************************************) BlumRandomPrime[n_Integer] := 
Module[{x, p}, 
        x = Random[Integer, {2^(n - 1), 2^n - 1}]; 
p = 
          BlumPrime[x]; 
If[
          p ≥ 2^n,  (*wrap around*)  
          x = 2^(n - 1); 
p = 
            BlumPrime[x]
]; 
Return[
          p]
]

(*******************************************************) 
 (* 
      Generate a BLUM number with n or n + 
        1 Bits.          *) 
 (*   (A BLUM number (after Leonore and Manul BLUM) is  *) 
 (*   a product of two primes congruent 3 mod 4.         *) 
 (* 
      Here we take two BLUM primes that differ 1 or 2 in  *) 
 (*   Bit size.)                                        *) 
 (*******************************************************) BlumNumber[
        n_Integer] := 
Module[{k, l, p, q, 
          m}, 
k = Quotient[n, 2]; 
l = 
          n - k + 1; 
p = 
          BlumRandomPrime[k]; 
q = 
          BlumRandomPrime[l]; 
m = 
          p*q; 
Return[
          m]
]

\(\(\( (*******************************************************) \)\(
\)\( (*\ 
      The\ BBS\ \(\(generator\)\(.\)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) \)\(
\)\( (*\ 
      Generate\ a\ list\ of\ n\ pseudorandom\ bits\ by\ n\ calls\ \ \ *) \)\(
\)\( (*\ \ "square x, extract lowest significant bit."\ \ \ \ \ \ \ \ *) \)\(
\)\( (*******************************************************) \)\(
\)\(BBSrandomBit[
        n_Integer] := 
Module[{b, 
          blist = {}}, 
Do[
x = 
            PowerMod[x, 2, m]; 
b = 
            Mod[x, 2]; 
AppendTo[blist, b], {i, 
            n}]; 
Return[
          blist]
]\)\)\)

\(\(\( (*******************************************************) \)\(
\)\( (*\ The\ BBS\ generator, \ 
      more\ efficient\ \(\(version\)\(.\)\)\ \ \ \ \ \ \ \ \ \ *) \)\(
\)\( (*\ 
      Generate\ a\ list\ of\ n\ pseudorandom\ bytes\ by\ n\ calls\ \ *) \)\(
\)\( (*\ \ "square x, extract 8 lowest significant bits."\ \ \ \ \ *) \)\(
\)\( (*******************************************************) \)\(
\)\(BBSrandomByte[
        n_Integer] := 
Module[{b, 
          blist = {}}, 
Do[
x = 
            PowerMod[x, 2, m]; 
b = 
            Mod[x, 256]; 
AppendTo[blist, 
            b], 
{i, n}]; 
Return[
          blist]
]\)\)\)

m  =  
      BlumNumber[1025]

m0 = N[Ceiling[Sqrt[m]]]

x = 
      Random[Integer, {m0, m - m0}]

ll = BBSrandomBit[1024]

DensityPlot[ll[[32*i + j]], {j, 1, 32},  {i, 0, 31},  
      PlotPoints  ->  32]

Export["C:\Temp\BBS.eps", %]

by = BBSrandomByte[128]