NetNews Usenet Archive 1992 #31/NN_1992_31.iso/spool/comp/lang/fortran/4820

Text File | 1992-12-21 | 47.3 KB | 1,038 lines

Newsgroups: comp.lang.fortran
Path: sparky!uunet!spool.mu.edu!umn.edu!msi.umn.edu!tregilli
From: [e-mail] (Ian Lee Tregillis [Astro])
Subject: Compilation of RNG responses
Message-ID: <[e-mail]>
Sender: [e-mail] (Usenet News Administration)
Nntp-Posting-Host: s8.msi.umn.edu
Organization: University of Minnesota
Date: Mon, 21 Dec 1992 16:32:32 GMT
Lines: 1026

Greetings, all-
    Here it is, finally.  Though not as promptly as I should have,
I've finally posted a summary of the responses I received regarding
my request for a Fortran random-number generator.  I've left out
the requests for me to share all of this stuff;  I will mail those people
directly.  Also, I edited out extraneous material here and there to
cut down on length;  deleted material is denoted with '[...]'.  Finally,
messages are partitioned with a line of asterisks for ease in reading.
    Again, thanks to everybody who responded!

Ian "Magnetic Moment" Tregillis
[e-mail]
[e-mail]

>>Begin included text<<
*************************************************************************
Date: Mon, 16 Dec 1991 23:11:56
From: "Vladimir Vasilchenko" <[e-mail]>
Organization: Mech & Appl Math Res Institute

Hi!
Is it possible for You to test programs from the IMSL (Statistics) Library.
This FORTRAN programs library I think is excellent tool for the specialists
in statistics ( random number generators and so on .... ). Also I have some
news about widely available statistics fortran library in Japan. My friend
received so many programs from this library via Internet network. Moreover,
You can try to find such programs in NAG library ( Oxford ), Cambridge fortran
library ( Numerical Recipes ) and in the book "Computer methods of mathematical
calculations" by Malcolm, Moler, Forsait ( Unfortunately, I don't remember
the English form of the second names of these authors and here I retranslate
ones from the Russian translation of this book. Sorry!).
Best Regards,
Vladimir V.Vasil'chenko, Russia, 16 December 1992.
***************************************************************************
Date: Mon, 14 Dec 92 14:23:09 -0500
From: "John D. McCalpin" <[e-mail]>
Organization: College of Marine Studies, U. Del.

[...]
Try anonymous ftp to perelandra.cms.udel.edu.
Look in pub/Math/random/.
--
John D. McCalpin                        [e-mail]
Assistant Professor                     [e-mail]
College of Marine Studies, U. Del.      [e-mail]
***************************************************************************
Date: Mon, 14 Dec 92 16:21:17 CST
From: [e-mail] (Ron Shepard)

Ian,  Look in "Computers in Physics, Vol. 6 No. 5 Sep/Oct 1992"
on pp. 522-524.  At the top of page 523, there is a 12-line
function ran0() that will probably work for your simple case.

Regards, -ron shepard
****************************************************************************
Date: Mon, 14 Dec 92 17:28:53 PST
Cc: [e-mail] (david hall;WR5756)

THE source now is Numerical Recipes, second edition.  Skip the first edition;
the random number section has been totally revised.  Or, the equivalent is in
Computers in Physics, 6(5), Seot-Oct 1992, Numerical Recipes section, "Portable
Random Number Generators," pp. 522-524.     I have lots of other references,
but these are handy, and in FORTRAN.

Press, William H., Teukolsky, Saul A., Vetterling, William T., Flannery, Brian P.
1992.  Numerical Recipes in FORTRAN: the art of scientific computing. 2nd ed.
Cambridge University Press.  QA297.N866

Would like to see what other references you are given...

David Hall - Programmer/Analyst (Masters thesis on RNG's).
*****************************************************************************
Date: Mon, 14 Dec 92 22:40:20 CST
From: Jeff Haferman <[e-mail]>

Actually, most "real" programmer don't always learn how to make
"decent" random number generators.  To find really decent fortran
source, find a copy of the newest edition of _Numerical_Recipes_in_
Fortran_.  The old edition has some random number generators, but
the latest edition has improved code.  Also, if you can find the
magazine "Computers in Physics", a month or two ago the authors
of the aforementioned book wrote a wonderful column on random
number generation, complete with Fortran source (very short).
[...]
Jeff Haferman
*****************************************************************************
Date: Mon, 14 Dec 92 23:29:51 CST
From: Bradbury Franklin <[e-mail]>

Mr. Tregillis,

     The following subroutine takes the variance n and the global index J
     as input and outputs m a normal random variable with mean 0 and variance
     n.  Note that J is updated with each call of wiener, and the _main program
     variable 'actual' parameter passed in the call of wienar should not be
     touched by anything else in the program-- it is not a dummy variable.
     (Any initial integral value for J will do. This value is the random seed.)
     Oh, and this will give 2**31-2 or so different numbers before cycling.

   Note that this is a rather poor random generator in that it exhibits
   clear lattice structure.   I have some better ones in C.  I will write it
   in Fortran someday not too long from now...

   Also, see L'Ecuyer, P. "Efficient and Portable Combined Random Number
   Generators" Communications of the Assn. for Comptng Machin. 31 (1988) pp.742-
   774.

    SUBROUTINE wiener(m,n,J)
       double precision m,n,a(2)
       integer J,Q,R,S,U,V
    R=7**5
    C=2**31-1
    do 120 S=1,2
    if (J.GT.2**31/R) then
      Q=J*R
    else
      V=MOD(J,127773)
      U=(J-V)/127773
      Q=-2836*V+R*U
    end if
    if (MOD(Q,C).LT.0) then
     J=-MOD(Q,C)
    else
     J=MOD(Q,C)
    end if
    a(S)=float(J)/float(C)
120     continue
       m=sqrt(-2*n*dlog(a(1)))*dsin(6.2832*a(2))
       end

Brad Franklin
[e-mail]
***************************************************************************
Date: 15 Dec 92 14:35:33+0800
From: Peter Hermann <[e-mail]>

Ian,

[...]
Have a look at Communications of the acm -- Nov92.

Regards,
Peter
*****************************************************************************
Date: Tue, 15 Dec 92 14:53:56 +0100
From: [e-mail] (Reinhard     Drube          E1
      )

Ian,
I have one more hint for you to test the selected random number
generator: Try to plot the produced numbers vs. the last number.
On bad RND-generators you will see dependencies of the numbers
in forms of circles, stripes and so on. A good generator will only
produce random noise!  Again good luck!     Reinhard
*****************************************************************************
Date: Tue, 15 Dec 92 10:22:42 -0500
From: Ji Wang <[e-mail]>
Organization: Princeton University

Well, the easiest way is to check a library like NAG or IMSL, which
include a lot of subroutines for this purpose.  It is really easy
to call.  And, check books on FORTRAN, you may find generator in many
books.  Last, you definity have it n Numerical Recipes.

Ji Wang
*****************************************************************************
Date: Tue, 15 Dec 92 10:51:42 EST
From: [e-mail]

I have enclosed a posting I saw a while ago about a library of
random number generators from statlib.  Both statlib and netlib
are vast repositories of code, and are really worth checking out.

netlib:
    email [e-mail] or [e-mail]
        with body of message;
            send index
or
    ftp  to research.att.com, login as netlib, and give e-mail
    address as password.

statlib instructions given below.

Erik Iverson
[e-mail]

Begin included text...

                      Fortran version of RANLIB


RANLIB is a collection of  Fortran routines that provide generators of
random numbers from a variety of distributions.  RANLIB uses published
algorithms, where available  (literature citations are included in the
documentation).

RANLIB  is  available  via  anonymous   FTP from  odin.mda.uth.tmc.edu
(129.106.3.17).  It is
                        /pub/unix/ranlib.f.tar.Z

RANLIB   is  also   available   by   electronic   mail (or  ftp)  from
[e-mail].   The  file 'ranlib'  contains   an expanded
description of the package; the file, 'ranlib.shar' contains  the full
code and documentation.  To obtain both,  send mail with the following
two lines to [e-mail].
         send ranlib from general
         send ranlib.shar from general
Those who are unfamiliar with statlib might want to add a third line
to obtain an introductory document.
         send index

                    Summary of Ranlib Capabilities

The bottom level routines provide 32 virtual random number generators.
Each generator can provide 1,048,576 blocks of numbers, and each block
is of length 1,073,741,824.  Any generator can be set to the beginning
or end  of the current  block or to  its starting value.  Packaging is
provided   so  that  if  these capabilities  are not  needed, a single
generator with period 2.3 X 10^18 is seen.

Using this base, routines are provided that return:
    (1)  Beta random deviates
    (2)  Chi-square random deviates
    (3)  Exponential random deviates
    (4)  F random deviates
    (5)  Gamma random deviates
    (6)  Multivariate normal random deviates (mean and covariance
         matrix specified)
    (7)  Noncentral chi-square random deviates
    (8)  Noncentral F random deviates
    (9)  Univariate normal random deviates
    (10) Random permutations of an integer array
    (11) Real uniform random deviates between specified limits
    (12) Binomial random deviates
    (13) Poisson random deviates
    (14) Integer uniform deviates between specified limits
    (15) Seeds for the random number generator calculated from a
         character string
*****************************************************************************
Date: Tue, 15 Dec 92 11:45:00 -0800
From: [e-mail] (Luke Dones)

Hi --

 I read two of your recent postings on random number generators and was
wondering if you'd mind sharing the code/lore you'd collected. I recently had
the unpleasant experience of reading an admission in Computers in Physics
(Sep/Oct '92) by Press and Teukolsky that some of the random number generators
in the first edition of Numerical Recipes, which I had used
extensively, had problems (despite their previous lofty claims for how good
they were!) Anyway, in CIP they give three new, supposedly improved codes,
which I enclose with this message (I have verified that they run, but have
not tested them extensively). Even though they are now betting $1000 on the
reliability of their new code, I still worry, so I would like to try some other
routines.

Another possibly useful article appears in Dr. Dobb's Journal, Feb. 1992, pp.
34-40. This article gives code (in C) for random number generators with periods
of ~4E9 and 1E27, if the author is to be believed....

Regards,
 Luke Dones
-----------------------------------------------------------------------------

      FUNCTION RAN0(IDUM)
c    new routine: Computers in Physics 6, 522-524 (Sep/Oct 1992)
c    "minimal standard" with period of 2**31 -2 ~ 2.1E9,
c    but has serial correlations
c    set idum = any integer except mask to initialize sequence
      INTEGER IDUM,IA,IM,IQ,IR,MASK
      REAL RAN0,AM
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,
     *    IQ=127773,IR=2836,MASK=123459876)
      INTEGER K
      IDUM=IEOR(IDUM,MASK)
      K=IDUM/IQ
      IDUM=IA*(IDUM-K*IQ)-IR*K
      IF (IDUM.LT.0) IDUM=IDUM+IM
      RAN0=AM*IDUM
      IDUM=IEOR(IDUM,MASK)
      RETURN
      END

      FUNCTION RAN1(IDUM)
c    new routine: minimal standard, shuffled
c    claimed to be adequate for up to roughly 10**8 random deviates
c     (5% of period)
c    set idum = any negative integer to initialize sequence
c    rnmx = "largest floating point value that is less than 1" --
c     this looks dubious
      INTEGER IDUM,IA,IM,IQ,IR,NTAB, ndiv
      REAL RAN1,AM, eps, rnmx
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *    NTAB=32,ndiv=1+(im-1)/ntab, eps = 1.2e-7, rnmx = 1. - eps)
      INTEGER J,K, iv(ntab), iy
      save iv, iy
      data iv /NTAB*0/, iy /0/
      if(idum.le.0 .or. iy.eq.0) then
        IDUM=MAX(-IDUM,1)
        DO J=NTAB+8,1,-1
          K=IDUM/IQ
          IDUM=IA*(IDUM-K*IQ)-IR*K
          IF (IDUM.LT.0) IDUM=IDUM+IM
          if(j.le.NTAB)iv(j)=idum
    enddo
        iy=iv(1)
      ENDIF
        K=IDUM/IQ
        IDUM=IA*(IDUM-K*IQ)-IR*K
        IF (IDUM.LT.0) IDUM=IDUM+IM
        j = 1+iy/ndiv
        iy=iv(j)
        iv(j) = idum
        ran1 = min(am*iy,rnmx)
      RETURN
      END

      FUNCTION RAN2(IDUM)
c    Long-period (> 2E18) random number generator of L'Ecuyer
c     with Bays-Durham shuffle
c    set idum = any negative integer to initialize sequence
      INTEGER IDUM,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB, ndiv
      REAL RAN2,AM,eps,rnmx
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *    IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,
     *    IR2=3791,NTAB=32,ndiv=1+imm1/ntab,eps=1.2e-7,
     *    rnmx=1.-eps)
      INTEGER IDUM2,J,K,IV(NTAB),iy
      SAVE IV,IY,IDUM2
      DATA IDUM2/123456789/, IV/NTAB*0/, IY/0/
      IF (IDUM.LE.0) THEN
        IDUM=MAX(-IDUM,1)
        IDUM2=IDUM
        DO J=NTAB+8,1,-1
          K=IDUM/IQ1
          IDUM=IA1*(IDUM-K*IQ1)-K*IR1
          IF (IDUM.LT.0) IDUM=IDUM+IM1
          if(j .le. ntab)IV(J)=IDUM
    enddo
        IY=IV(1)
      ENDIF
        K=IDUM/IQ1
        IDUM=IA1*(IDUM-K*IQ1)-K*IR1
        IF (IDUM.LT.0) IDUM=IDUM+IM1
        K=IDUM2/IQ2
        IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2
        IF (IDUM2.LT.0) IDUM2=IDUM2+IM2
        J=1+IY/ndiv
        iy=iv(j)-idum2
        IV(J)=IDUM
        IF(iy.LT.1)iy=iy+IMM1
        RAN2=min(am*iy,rnmx)
      RETURN
      END
[...]
****************************************************************************
From: [e-mail]
Date: Tue, 15 Dec 1992 22:36:44 -0500 (EST)

Here is a good random number generator.  Look at the reference!
If you want gaussian, then the simplest way is to add 12 uniforms if
accuracy is not that important.
--
Edward E.L. Mitchell
Mitchell & Gauthier Associates | Internet: mitchell@.mga.com
[...]

------------------------------------zzurv.f------------------------
      subroutine zz urv(seed, output)
c--------------------------------------random number generator based on
c     algorithm proposed by d. h. lehmer. it is known formally as a
c     prime modulus multiplicative linear congruential generator
c     (pmmlcg), or less formally as a lehmer generator. derived from an
c     implementation appearing in: park, stephen k. and miller, keith w.
c     "random number generators: good ones are hard to find", commun.
c     acm 31,10 (october, 1988), pp. 1192 - 1201.
c
c     as described in this paper, this implementation is based on
c     schrage's method: schrage, l. "a more portable fortran random
c     number generator". acm trans. math. softw. 5, 2(june 1979), pp.132
c
c     m        fixed prime modulus equal to (2**31 - 1)
c     a        multiplier that produces a full period generating
c              function f(z)=az mod m
c     q        m div a
c     r        m mod a
c     test     can never take on a value which cannot be represented
c              correctly with 32 bits (including sign)
c     seed     any integer between 1 and (2**31 - 2)
c     output   always REAL not DOUBLEPRECISION
c
      integer        seed    , hi      , test    , a       , q
     .             , r
c
      parameter    ( a = 16807         , m = (2**30 - 1)*2 + 1
     .             , q = m/a           , r = m - q*a
     .             , floatm = m
     .             )
c
      hi      = seed/q
      test    = a*(seed - hi*q) - r*hi
      if (test .gt. 0)then
         seed = test
      else
         seed = test + m
      endif
c-----------------------------return result between zero and one
      output  = seed/floatm
      return
c
      end
------------------------------------zzgrv.f------------------------
      subroutine zz grv(seed, output)
c-----------------------------gaussian random variable. adds twelve
c     uniform random numbers between 0 and 1.0 and subtracts average
c     NOTE. output is always REAL, not DOUBLEPRECISION
c
c     copyright (c) mitchell and gauthier associates (1989)
c     all rights reserved
c
c     inputs
c
c     seed     integer seed that will be changed
c
c     outputs
c
c     output   real random number (sum of 12 uniforms)
c
      integer        seed
c
      grv     = 0.0
      do 110 j = 1, 12
c-----------------------------use uniform random number generator
      call zz urv(seed, urv)
      grv     = grv + urv
  110 continue
c
c-----------------------------subtract average so mean is zero
      output  = grv - 6.0
      return
c
      end
****************************************************************************
Date: Wed, 16 Dec 92 11:39:26 -0700
From: [e-mail] (Lou Glassy)

Ian, this is a random number generator I typed in from
_Numerical Recipes_.  It's reasonable, and has a period
of about 700 000 or so, ie the output numbers will repeat
every 700 000 numbers or so.  This routine is standard
Fortran-77, and should compile on any decent f77 compiler.
I have tried it on MSDOS (BCF77), VMS/VAX FORTRAN, DEC/MIPS Fortran,
HPUX Fortran, and IBM AIX Fortran, and it has worked on these
OS/compiler combinations.
[...]
C=== delete everything above this line.

        REAL FUNCTION RAN2(IDUM)
C returns a uniform random deviate between 0.0 and 1.0.
C set idum to any negative value to initialize or
C reinitialize the sequence.
C source:  Numerical Recipes: The Art of Scientific Computing,
C          by W. H. Press, 1989.

C persistent variables
        SAVE IFF, IR, IY

C parameters
C       inout
        INTEGER IDUM

C constants
        INTEGER M, IA, IC
        REAL RM
        PARAMETER (M=714025, IA=1366, IC=150889, RM=1.0/M)

C variables
        INTEGER IFF, IR(97), IY, J
        DATA IFF /0/

C begin
        IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN
           IFF = 1
           IDUM = MOD(IC - IDUM, M)

           DO 11 J = 1, 97
              IDUM = MOD(IA * IDUM + IC, M)
              IR(J) = IDUM
 11        CONTINUE
           IDUM = MOD(IA * IDUM + IC, M)
           IY = IDUM
        ENDIF

C here's where we start except on initialization.
        J = 1 + (97 * IY) / M
        IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE
        IY = IR(J)
        RAN2= IY * RM
        IDUM = MOD(IA * IDUM + IC, M)
        IR(J) = IDUM

        END
***************************************************************************
Date: 18 Dec 92 20:34:42+0800
From: Peter Hermann <[e-mail]>

Ian,
it happened by chance, that the following mail
arrived to me, just in time. Have fun and insight.
merry christmas,
Peter
[...]
1.  Introduction

In designing a facility for the generation of pseudo-random numbers (hereafter
called simply "random numbers") in Ada 9X, we have sought to balance utility
and simplicity.  We do not offer all the bells and whistles someone might want
or need, but we do offer the capabilities that are difficult to program
correctly from scratch; other capabilities can be composed on top of these
rather easily by the user, if required.

We present in section 2 the specification of a package providing the desired
capabilities.  Section 3 discusses the design decisions leading to this
specification, while section 4 discusses some other alternatives that might
still be considered.  Section 5 discusses properties that the implementation
must satisfy.  We conclude by inviting comments from reviewers.


2.  Specification

The specification of the random number facility for Ada 9X is as follows:

   package Random_Numbers is

      type Seed is <implementation-defined>;  -- but not private
      Seed_Error : exception;                 -- see Set_Seed, below

      type Reset_Integer is
         range <implementation-defined>;      -- see Reset, below

      generic
         type Float_Type is digits <>;
      package Generic_Generator is

         -- Instantiation of this generic package provides a random number
         -- generator in each task; multiple instantiations provide multiple
         -- generators.

         function Random return Float_Type;
           -- Get the "next" random number, uniformly distributed on
           -- [0.0, 1.0), from the current task's generator.

         procedure Reset (I : in Reset_Integer);
           -- Set the state of the current task's generator to an
           -- implementation-defined function of the integer I.

         procedure Reset;
           -- Set the state of the current task's generator to an
           -- implementation-defined function of the current date and time.

         procedure Get_Seed (S : out Seed);
           -- Get the state of the current task's generator.

         procedure Set_Seed (S : in  Seed);
           -- Set the state of the current task's generator to an explicit,
           -- user-provided value; raises Seed_Error if S is unacceptable to
           -- the implementation.

      end Generic_Generator;

   end Random_Numbers;


3.  Design decisions

As in many other languages, only one kind of distribution is provided, namely
the uniform distribution.  More specifically, the random numbers generated by
the facility are uniformly distributed on the semi-open interval [0.0, 1.0),
and are delivered as values of a floating-point type chosen by the user.  Other
intervals are easily obtained by scaling and translation:  if X is uniformly
distributed on [0.0, 1.0), then a*X+b is uniformly distributed on [a, b).
Integer values uniformly distributed on some interval are also easily obtained
by appropriate scaling and conversion to integer, with post-conversion
adjustment if rounding went "the wrong way."

Other distributions, such as the normal (Gaussian), Poisson, binomial, or
exponential distributions, can be obtained from the uniform distribution by
various standard techniques.  We focus our efforts on the capability most
difficult to provide if it is not predefined, namely a reliable uniform random
number generator.

Some algorithms naturally produce 0.0 among the numbers generated, and it is
normal to include this value in the range of the random number generator.
Mathematically, good algorithms exclude 1.0, but care must sometimes be
exercised to avoid rounding a result close to 1.0 up to 1.0.  We could make the
job simpler for implementors by including 1.0 in the specified range of random
numbers but choose not to on account of the benefit to some applications of
knowing that at least one of these extremes is never generated.  For example,
if X is uniformly distributed on [0.0, 1.0), then -log(1.0 - X) is
exponentially distributed with mean 1.0 and standard deviation 1.0, and it
never receives an invalid argument.

A primary goal of the design of the random number facility for Ada 9X is the
preservation of opportunities for innovation on the part of implementors, as
the state of the art advances, while providing for a certain amount of
portability of applications among implementations.  Thus, an algorithm is not
prescribed, nor is the nature of the seed specified.  (The seed, which may be
scalar or composite, embodies the state of a generator.  Technically, we should
reserve "seed" to mean the initial state of a generator, but since a
generator's current state may be though of as its initial state relative to the
future sequence of outputs, we use seed to mean its current state.)  It is our
intention to allow any of the following as implementations of the random number
generator:

   o  good instances of a classical linear congruential generator, whose scalar
      seed is an integer or, in some realizations, a floating-point number;

   o  combination generators, such as those of Wichmann and Hill or L'Ecuyer,
      whose composite seed is a small vector;

   o  the more recent "subtract-with-borrow" or "add-with-carry" Fibonacci
      generators of Marsaglia and Zaman that have phenomenally long periods,
      but whose composite seed is a vector, often of twenty to a hundred
      integer or floating-point components representing the lagged previous
      outputs used in computing future outputs, plus another component
      representing the next borrow or carry.

Some assurance of quality is provided by specifying minimal properties that the
implementation must exhibit (see section 5).  Portable means of obtaining
reproducible sequences of random numbers within a given implementation are
provided, as are portable means of obtaining unique sequences of random numbers
from one run to the next in a given implementation.  (Reproducible sequences
are generally desired for application development and testing, while unique
sequences are often wanted for production use of applications.)  The user must
not expect to be able to obtain the same "reproducible" sequence of random
numbers he or she has been accustomed to in one implementation upon switching
to a different implementation, or when an implementation is upgraded by the
vendor.

Generators have a default (implementation-defined) initial state, which
provides for reproducibility of sequences across runs in simple applications.
Unique sequences in different runs are obtained by replacing the default
initial state with a time-dependent state (by calling the parameterless version
of Reset before generating random numbers).

Normally, one does not need to manipulate the state of a generator (i.e., its
seed) directly, but Get_Seed and Set_Seed are provided for those occasions
when, in more advanced applications, it is desirable to do so.  For example,
the current state of a generator can be obtained by a call to Get_Seed and
saved in a file; the value retrieved from the file can be passed to Set_Seed in
a subsequent run to resume the sequence from the previous stopping point.  This
particular use does not require knowledge of the seed type, but other uses of
Get_Seed and Set_Seed do and are inherently implementation dependent.  One
example of such a non-portable use of Set_Seed is to obtain a reproducible
sequence other than the one provided by default.  Another arises in
applications in which different tasks have separate random number generators
that must produce different, but reproducible, sequences.  The only general
solution to this problem is to allow the user to compute an initial seed for
the generator in each task, perhaps as a function of the current seed in the
task's parent.  An example of a non-portable use of Get_Seed is its use in the
simple printing of the current seed for debugging purposes.

The version of Reset with an integer parameter provides an alternative way to
obtain a reproducible sequence other than the default sequence.  In contrast to
using Set_Seed for that, Reset with an integer parameter is somewhat more
portable though less general.  It is more portable only in the sense that the
source text of some applications, for example those that prompt the user for a
seed value, need not depend on the implementation.  It is less general in the
sense that the possible values of a composite seed are far more numerous than
those of a single integer, so only a subset of the possible states of a
generator may be produced by calling Reset.  Every integer in the range of
Reset_Integer should be acceptable to Reset.

The seed type is a visible type in order to support the convenient computation
of seed values, when required.  Implementations are allowed to define other
types in terms of which Seed is defined, as may be necessary when Seed is
composite; they may also define Seed by a subtype declaration instead of a type
declaration.

There are several different ways to enable separate tasks to have their own
random number generators.  The method we have chosen has some limitations,
discussed along with the alternatives in section 4, but provides within those
limitations for the greatest simplicity of use.  Each instantiation of
Generic_Generator provides each task with its own generator of random numbers,
which is implicitly accessed by the operations Random, Reset, Get_Seed, and
Set_Seed.  (They access the generator of the current task.)  This behavior can
be implemented by appropriate instantiation of the predefined generic package
Task_Identification.Per_Task_Attributes in the body of Generic_Generator (the
seed becomes a per-task attribute).  Those single-task applications requiring
multiple generators and those multitasking applications requiring multiple
generators in some or all tasks can obtain them by instantiating
Random_Numbers.Generic_Generator multiple times.

Some values of type Seed may be inappropriate for certain generators.  In some
cases, inappropriate seed values cannot be excluded simply by the use of range
constraints in the definition of Seed or its components.  The exception
Seed_Error is provided for Set_Seed to use in signaling an attempt to set the
state of a generator to an inappropriate value.  Implementations may prefer to
avoid or reduce the possibility of inappropriate seed values, which arise
largely through the use of floating-point seeds, by defining the seed type to
be an integer or a composite of integer components, even when the internal
state contains floating-point components (having integer or scaled integer
values).

We do not provide a subprogram to fill a vector with random numbers in a single
call.  On scalar computers, such a capability offers a reduction in the
subprogram call overhead.  On vector computers, additional benefits might
accrue by vectorizing the implementation, which is possible for some
algorithms; the primary benefit on such computers, however, is the additional
vectorization of the user's code that might be possible when, in a loop, random
numbers are consumed from a vector rather than by calling a subprogram to
generate each one.  In omitting this capability, we emphasize the minimal
essential functionality of our random number facility.  Omitting it also keeps
the generic formal part simple.


4.  Alternatives

An alternative solution to the problem of obtaining different but reproducible
sequences in different tasks, one that does not require knowledge of the seed
type, is a procedure to "advance" a generator over some presumably large
portion of its full period without generating the intervening elements.  This
is relatively straightforward for linear congruential generators, but difficult
or impossible for some of the other kinds of generators we wish to allow.
Furthermore, obtaining many "different" sequences by systematically advancing
their sequences by given amounts runs the risk of obtaining identical sequences
offset in time, or, in the case of linear congruential generators, more weakly
correlated sequences.  Thus, it seems desirable to give the user full latitude
in deciding how to obtain reproducibly different sequences in different tasks,
even though it means requiring the use of implementation-dependent procedures
for doing so.

We have found in the literature a reference to a library at Queen's University,
Belfast, whose subtract-with-borrow Fibonacci generator has a procedure for
resetting the state of the generator from an integer, each value of which
initiates a very long subsequence of the full period not overlapping with any
other (in any practical sense).  We have not yet seen the theoretical basis for
such a claim, and we do not know how to implement it ourselves.  Our Reset
procedure (with an integer parameter) can provide exactly this functionality
for appropriate generators; for such generators, at least, it provides a simple
method of giving each task its own sequence of random numbers.

The simple printing of the current seed for debugging purposes could be made
easier by something like an Image function that returns the state of a
generator as a string according to some implementation-defined coding scheme.
Indeed, initially we made the seed type private and contemplated Image and
Value functions on seeds.  It does not seem likely that Image and Value would
be suitable replacements for Get_Seed and Set_Seed; using them would, for many
purposes, be inconvenient because of the need to code or decode the strings (in
the case of a composite seed), and they would offer no real escape from the
inherently implementation-defined nature of seeds.  We fear that they would
give a false sense of being able to manipulate seeds portably, given that the
implementation-defined characteristics of the interface would disappear from
the specification and remain hidden away in the semantics.  We therefore omit
Image and Value as not offering sufficient added benefit, given Get_Seed,
Set_Seed, and a visible seed type.

The most difficult problem we faced in designing these facilities was their
"packaging."  We could easily have done something different, and we are still
not completely sure that we shouldn't.

The main attraction of our packaging method is the painless way in which each
task is given its own generator of random numbers.  The user need not create
objects in each task for the storage of the task's generator's state, and no
extra parameters are needed to identify the generator being addressed in the
various operations on generators; in particular, no parameters are passed to
Random, which is therefore as easy to use in separate tasks (where separate
generators are wanted) as in non-tasking applications.  It is anticipated that
implementations will instantiate Task_Identification.Per_Task_Attributes
internally to create a "slot" in each task for the task's random number
generator's seed; this use of Per_Task_Attributes has been mentioned repeatedly
during the design of that facility.

Some small number of applications might require multiple generators per task,
or multiple generators in a single task.  They can be obtained by instantiating
Random_Numbers.Generic_Generator multiple times.  This suffices if a fixed
number of generators is needed, so that they can be addressed without using
indexed names.  (If there are two generators per task, and the two
instantiations of Generic_Generator are called, say, Velocity_Generator and
Mass_Generator, then the random numbers from the two generators could be
obtained by calling Velocity_Generator.Random and Mass_Generator.Random.)  Our
scheme does not support applications requiring an arbitrary number of
generators in a single task or in multiple tasks, where indexing is required.
We do not view this as a serious limitation.

(There may be no overhead for random number generators in tasks that never
invoke a random number operation; it depends on the implementation of
Per_Task_Attributes.  That is, an instantiation of
Random_Numbers.Generic_Generator does not necessarily cause the acquisition of
space for the generator's seed; that happens, at least with the canonical
implementation model of Per_Task_Attributes, only if and when an operation is
first performed on a generator.  There is, however, a little extra overhead for
using random number generators in our proposed scheme than in the alternatives
discussed below.)

The most obvious way to allow the creation of arbitrary numbers of random
number generators without the use of Per_Task_Attributes is to identify a
generator with an object of type Seed.  A generator would be created by
elaborating an object declaration for such an object.  (The declaration of the
type Seed would be changed to specify the initial state of the generator as the
initial value of objects of the type.)  The user would have the responsibility
for creating generators in the quantity required, and in the just the places,
i.e. tasks, required.  Operations on a particular generator would be designated
by passing the seed as an "in" parameter to the operation ("in out" in the case
of Random, which would have to become a procedure, since Random updates the
state).  There would be no need for Get_Seed and Set_Seed operations.

We never contemplated actually making a seed the embodiment of a generator.
Were we to follow something like the approach discussed immediately above, we
would, instead, embody a generator in a separate type.  Our specification might
be as follows:

   with System;
   package Random_Numbers is

      type Generator is limited private;

      type Seed is <implementation-defined>;  -- but not private
      Seed_Error : exception;

      type Reset_Integer is range <implementation-defined>;

      generic
         type Float_Type is digits <>;
      function Generic_Random (G : Generator) return Float_Type;

      procedure Get_Seed (G : in Generator; S : out Seed);
      procedure Set_Seed (G : in Generator; S : in  Seed);
      procedure Reset (G : in Generator; I : in Reset_Integer);
      procedure Reset (G : in Generator);

   private

      type State is <implementation-defined>;  -- typically a record containing
                                               -- a Seed component and maybe
                                               -- other stuff
      type Access_State is access State;
      type Generator is new System.Controlled with record
         S : Access_State := new State'(<implementation-defined>);
      end record;
      procedure Finalize (G : in out Generator);

   end Random_Numbers;

(Deriving Generator from System.Controlled allows finalization of generators
and, thereby, the reclaiming of the storage allocated separately to their
associated states.)

In this alternative, as in the proposal of sections 2 and 3, the state of a
generator is not directly accessible.  That operations are performed on
generators, not seeds, seems like the proper abstraction.  Get_Seed and
Set_Seed play the same role as before in obtaining or changing the internal
state of a generator (which now becomes a parameter to them).  Objects of type
Seed would be declared, as in the original scheme, only if the user needs to
obtain or save the state of a generator.  The validity check performed on a
seed being provided to a generator remains part of the semantics of Set_Seed,
where it is relevant.  If, in contrast, generators were directly associated
with seeds, and created by creating objects of type Seed, and operated upon by
operating on such objects, then Generic_Random would presumably have to
validate the current seed each time a random number is requested (if seed
validity needs to be checked), since the seed might have been changed by direct
manipulation since the instantiation of Generic_Random was last called.  That
would be intolerably expensive.

We would have to state that the generators provided by this alternative scheme
are single-threaded; it would be an error to share them among (call them from)
multiple tasks.  In the scheme proposed in sections 2 and 3, there is no need
to make such a statement, since each task's generator can be called only by its
own task.

In addition to being slightly more general than our original scheme, albeit
slightly harder to use in simple applications, the alternative has two
additional advantages.  One is that it can be implemented in Ada 83 (except for
the reclamation of the storage occupied by a generator's state).  The other is
that it can serve as the foundation of a higher-level abstraction that some
users may wish to implement, using protected types -- that of a multi-threaded
random number generator.  We do not consider the need for random number
generators that can be called by multiple tasks (with appropriate
serialization, of course) important enough to be predefined.  The
nondeterministic nature of (the sequence of) calls to such generators would, in
most cases, eliminate any possibility of reproducible behavior.  Nevertheless,
it would be good to allow the advanced programmer the opportunity, at least, to
create such a capability, if it is needed, out of more basic predefined
components.  With the facility proposed in sections 2 and 3, it appears that
providing such a capability requires the creation of a server task in which
random numbers are generated, and with which other tasks must explicitly
interact.  If control over the identity of generators were provided instead by
the Generator type, the user could embed a Generator object in a protected
object or protected type, as needed.

Are these advantages of the alternative scheme significant enough to warrant
adopting it instead?  Reviewer comment is solicited on this question.

It has been suggested that some combination of both schemes be provided, so
that simple applications can be served simply while not precluding more
advanced applications.  That is, one could provide BOTH the Generator type (and
operations on it) AND the inner generic package Generic_Generator.  (Text_IO
has been cited as a model; there, I/O operations can be performed either on the
"current default input or output file" or on arbitrary files.)  We have not
actually proposed to do so because the resulting facility would have
overlapping or redundant capabilities; most of the random number generators
that would be needed, especially those in simple applications, could be created
in either of two ways, which is confusing and unattractive.

It has been assumed that the need exists for the user to be able to specify the
floating-point subtype of the random numbers that are delivered.  This is the
reason why Random has become a generic function in the alternative proposal.
The analogue of this in section 2 is the inner generic package,
Generic_Generator, that contains the function Random.  There, however,
genericity plays a further role: it is the means by which multiple generators
are obtained in a task (more precisely, in each task).  Of course, the generic
actual subtype used in the instantiation should have a range of at least
0.0 .. 1.0 and should provide "adequate" precision; too little could introduce
repeated values into the output.  It may actually be preferable to standardize
the result subtype of Generic_Random as Float, or let it be an
implementation-defined floating-point subtype.  If this were to be done, then
Generic_Random could revert to an ordinary function in the alternative
proposal, and Generic_Generator could become a parameterless generic package in
the original proposal; the user would generally have to convert the random
numbers obtained to some application-dependent type.  (Evidence that this may
be appropriate comes from our experimentation with different generator
implementations; we observed that the conversion from implementation-dependent
types appropriate to the algorithm to the user's type, denoted by Float_Type,
occurs as the final step.)


5.  Properties required of implementations

Random number generators should have long periods.  It will be specified that
implementations must exhibit a period of at least 2**31-2.  (Unfortunately, it
is not practical for this property to be verified by a validation suite.)  An
ultra-fast computer may be able to process applications for which this period
is inadequate, in which case the vendor should consider using a
subtract-with-borrow or add-with-carry Fibonacci generator.

There needs to be adequate assurance that the version of Reset that uses the
current date and time will provide a unique sequence on two widely separated
runs.  It will be specified that this procedure must reset the generator to
different states when the interval between calls is at least one second and not
greater than fifty years; a failure of this uniqueness property is allowed only
when two calls are separated by exactly one hour and a return from daylight
savings time to standard time intervenes.  Furthermore, the function that maps
the current date and time into generator states should be a wildly varying
function (i.e., one whose sensitivity to the components of the current date and
time is greatest with respect to the second and least with respect to the
year).  By the same token, the version of Reset that takes a user-supplied
integer should map consecutive integers into generator states that initiate
significantly different, and ideally uncorrelated, sequences.

The most important properties that implementations must exhibit are those
related to the "randomness" (unpredictability to one not familiar with the
algorithm) and uniformity of the generated numbers.  The general approach to
uniformity testing is to formulate an experiment based on random sampling for
which the expected outcome for a truly uniform random sample can be calculated,
perform the experiment using the random number generator, calculate the
chi-square value for the experimental observations, and then check that it is
less than the chi-square value expected from theory (with some confidence like
95%).  Exact criteria are still being investigated, as are appropriate
experiments.  The latter can be as simple as distributing the random numbers
(or tuples of them) into cells of even width, which would be filled in equal
number by truly uniform random numbers, or they can be more sophisticated.  One
experiment that looks promising is the playing of a large number of "craps"
games, with rolls of the dice simulated by the random number generator, for
which measurements of wins and losses, game length, and the length of "passes"
(runs of consecutive wins) can be tested against theoretical data for real
dice.  This experiment is useful because correlations many dice rolls apart can
affect the duration of a game; it can be combined with simple equidistribution
tests on the outcomes of the dice rolls.  The literature describes a variety of
experiments that have been used with some success in the past.


Request for Review

Comments are solicited on any part of this proposal, but especially on the
following questions:

   o  Should the approach be like that outlined in sections 2 and 3 (in which
      the implementation is expected to use Per_Task_Attributes to give each
      task a random number generator for each instantiation of
      Generic_Generator), or should it be more like the alternative discussed
      in section 4 (in which the user creates objects of type Generator in the
      numbers and places desired), or should it be a combination of the two
      approaches offering both simplified use and general use?

   o  Was the decision against providing vector output from the random number
      generator correctly made and justified?

   o  Was the decision against making the seed type private, and against
      supporting a private type with Image and Value functions, correctly made
      and justified?

   o  Should the floating-point type of the random numbers be specified as
      Float, or specified by the implementation, instead of being selected by
      the user?


References

Publications that have proven useful during the preparation of this LSN include
the following:

   D. E. Knuth, "Seminumerical Algorithms" (vol. 2 of The Art of Computer
      Programming), Addison-Wesley, 1989.

   F. James, "A Review of Pseudorandom Number Generators," Computer Physics
      Communications 60:329-344, North-Holland, 1990.

   S. L. Anderson, "Random Number Generators on Vector Supercomputers and Other
      Advanced Architectures," SIAM Review 32(2):221-251, June 1990.

   G. Marsaglia and A. Zaman, "A New Class of Random Number Generators," Annals
      of Applied Probability 1(3):462-480, August 1991.

   Ada BCSLIB (Mathematical/Statistical Library), Report 20462-0519-R2, Boeing
      Computer Services, Seattle, 1990.

   B. A. Wichmann and I. D. Hill, "A Pseudo-Random Number Generator," Report
      DITC 6/82, National Physical Laboratory, Teddington, UK, June 1982.

   B. A. Wichmann and J. G. J. Meijerink, Report DITC 39/84, National Physical
      Laboratory, Teddington, UK, March 1984.

   P. Lecuyer, "Efficient and Portable Combined Random Number Generators,"
      CACM 31(6):742-749,774, June 1988.
*****************************************************************************
>>End included text<<