From: blakemor@software.software.org (Alex Blakemore)
Subject: Re: looking for Ada random number generator (source posted)
Date: 21 Mar 91 14:11:16 GMT [thread overview]
Message-ID: <1991Mar21.141116.23895@software.org> (raw)
-- > 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 22070
next reply other threads:[~1991-03-21 14:11 UTC|newest]
Thread overview: 2+ messages / expand[flat|nested] mbox.gz Atom feed top
1991-03-21 14:11 Alex Blakemore [this message]
1991-03-21 16:35 ` looking for Ada random number generator (source posted) Doug Smith
replies disabled
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox