(*******************************************************) (* 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 = 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] ] (*----------Example application-----------------------*) (* 1. Generate module. *) m = BlumNumber[1025] (* 2. Lower bound for start value *) m0 = N[Ceiling[Sqrt[m]]] (* 3. Generate random start value. *) (* WARNING: Because Mathematica's internal pseudo *) (* random generator is used, the start value is not *) (* true random as it should be. *) x = Random[Integer, {m0, m-m0}] (* Now generate BBS random Bits *) ll = BBSrandomBit[1024] DensityPlot[ll[[32*i+j]], {j,1,32}, {i,0,31}, PlotPoints -> 32] (* Alternatively generate BBS random bytes *) by = BBSrandomByte[128]