(*******************************************************)
(* 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]