From mboxrd@z Thu Jan  1 00:00:00 1970
X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me
X-Spam-Level: 
X-Spam-Status: No, score=0.6 required=5.0 tests=BAYES_05,INVALID_DATE
	autolearn=no autolearn_force=no version=3.4.4
Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!sol.ctr.columbia.edu!emory!gatech!mcnc!uvaarpa!software.org!blakemor
From: blakemor@software.software.org (Alex Blakemore)
Newsgroups: comp.lang.ada
Subject: Re: looking for Ada random number generator (source posted)
Summary: George Marsaglia's universal random number generator
Keywords: Ada random variates source
Message-ID: <1991Mar21.141116.23895@software.org>
Date: 21 Mar 91 14:11:16 GMT
Sender: blakemor@software.org (Alex Blakemore)
Organization: Software Productivity Consortium, Herndon, Virginia
List-Id: <comp.lang.ada>

-- > 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