blakemor@software.software.org (Alex Blakemore) (03/21/91)
-- > CML8@ROBIN.CS.UOFS.EDU asks (UNIX)
-- > One of the teachers in my department is looking for an random number generator
-- > written in Ada. He only needs it for integers, but a generic one would be fine too.
-- This algorithm appeared in Ada letters a few years ago, and one of
-- our word processors was kind enough to retype it (before OCR).
-- You may want to read the article that came with it and perform the
-- the tests it describes. I think it is claimed to give the same results
-- on any platform.
-- George Marsaglia's (FSU) universal random number generator algorithm
-- from Ada Letters v. VIII, n. 2 1988
generic
type real is digits <>;
package random_number is
m1 : constant := 179;
m2 : constant := m1 - 10;
subtype seed_range_1 is integer range 1 .. m1 - 1;
subtype seed_range_2 is integer range 1 .. m2 - 1;
default_i : constant := 12;
default_j : constant := 34;
default_k : constant := 56;
default_l : constant := 78;
procedure start (new_i : seed_range_1 := default_i;
new_j : seed_range_1 := default_j;
new_k : seed_range_1 := default_k;
new_l : seed_range_2 := default_l);
function next return real;
end random_number;
package body random_number is
m3 : constant := 97;
init_c : constant := 362436.0 / 16777216.0;
cd : constant := 7654321.0 / 16777216.0;
cm : constant := 16777213.0 / 16777216.0;
subtype range_1 is integer range 0 .. m1 - 1;
subtype range_2 is integer range 0 .. m2 - 1;
subtype range_3 is integer range 1 .. m3;
-- i, j, k must be in the range 1 to 178, and not all 1 while l may be an integer
-- ranging from 0 to 168. if nonnegative integer values are assigned to i, j, k or l
-- outside the specified ranges, the generator will still perform satisfactorily but
-- not product exactly the same bit patterns on different computers.
i, j, k : range_1;
ni, nj : integer;
l : range_2;
c : real;
u : array (range_3) of real;
-- this procedure initializes the table for
-- f (m3, range_3'last / 3 + 1, - mod 1.0), a lagged-fibonacci sequence
-- generator, and produces values for the arithmetic sequence.
procedure start (new_i : seed_range_1 := default_i;
new_j : seed_range_1 := default_j;
new_k : seed_range_1 := default_k;
new_l : seed_range_2 := default_l) is
s, t : real;
m : range_1;
begin -- start
i := new_i;
j := new_j;
k := new_k;
l := new_l;
ni := range_3'last;
nj := (range_3'last / 3) + 1;
c := init_c;
for ii in range_3 loop
s := 0.0;
t := 0.5;
for jj in 1 .. 24 loop
m := (((j * i) mod m1) * k) mod m1;
i := j;
j := k;
k := m;
l := (53 * l + 1) mod m2;
if (((l * m) mod 64) >= 32) then
s := s + t;
end if;
t := 0.5 * t;
end loop;
u (ii) := s;
end loop;
end start;
function next return real is
temp : real;
begin -- next
temp := u (ni) - u (nj);
if temp < 0.0 then
temp := temp + 1.0;
end if;
u (ni) := temp;
ni := ni - 1;
if ni = 0 then
ni := range_3'last;
end if;
nj := nj - 1;
if nj = 0 then
nj := range_3'last;
end if;
c := c - cd;
if c < 0.0 then
c := c + cm;
end if;
temp := temp - c;
if temp < 0.0 then
temp := temp + 1.0;
end if;
return temp;
end next;
begin -- random_number
-- initialize table u.
start;
end random_number;
--
---------------------------------------------------------------------
Alex Blakemore blakemore@software.org (703) 742-7125
Software Productivity Consortium 2214 Rock Hill Rd, Herndon VA 22070smithd@software.org (Doug Smith) (03/22/91)
> -- > CML8@ROBIN.CS.UOFS.EDU asks (UNIX) > -- > One of the teachers in my department is looking for an random number generator > -- > written in Ada. He only needs it for integers, but a generic one would be fine too. > I have been using this generator extensively and it seems to work fine. It is however intended for 32 bit machines, as indicated by the CACM articles referenced. I highly recommend these concise articles which should answer almost any question you might have. If you don't like the generic form of the random number generator, you can change it easily enough. I find that having a different generator instantiated in carefully selected sections of code can maintain a deterministic behavior despite debugging and enhancement. **************************** Cut Here ************************************** -- CACM, October 1988, pp. 1192-1201. -- CACM, January 1990, pp. 87-88. -- Implemented in Ada by D. Douglas Smith generic -- To be a full period generator, Multiplier must be a primitive root -- of the Prime_Modulus. This is not checked for, and if not a -- primitive root, then the Multiplier will generate only a subset -- (partition) of the Positive numbers depending on the Initial -- seed. -- Also, 2 <= Multiplier < Prime_Modulus and will otherwise raise -- Program_Error during elaboration; -- This "efficient" implementation requires -- (Prime_Modulus mod Multiplier) < (Prime_Modulus / Multiplier) -- and will otherwise raise Program_Error during elaboration; Prime_Modulus : in Positive := Positive'Last; -- 2**31 -1 = Mersenne prime; Multiplier : in Positive := 7**5; Initial_Seed : in Positive := 40; package Random is -- Note that 0 < Next_Positive < Prime_Modulus. function Next_Positive return Positive; -- (Multiplier * Seed) mod Prime_Modulus end Random; package body Random is a : constant Integer := Multiplier; m : constant Integer := Prime_Modulus; q : constant Integer := m/a; r : constant Integer := m mod a; seed : Integer := Initial_Seed; function Next_Positive return Positive is begin seed := a *(seed mod q) - r*(seed / q); if seed <= 0 then seed := seed + m; end if; return seed; end Next_Positive; begin if Multiplier < 2 or else Multiplier >= Prime_Modulus or else (Prime_Modulus mod Multiplier) >= (Prime_Modulus / Multiplier) then -- This implementation will not work. raise Program_Error; end if; end Random; **************************** Cut Here **************************************