Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/unix/volume28/mrandom-3.0/part03

Text File | 1994-05-06 | 98KB | 2,893 lines

Newsgroups: comp.sources.unix
From: [e-mail] (Clark D. Thomborson)
Subject: v28i029: mrandom-3.0 - random number generator with persistent state, Part03/06
References: <[e-mail]>
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (Clark D. Thomborson)
Posting-Number: Volume 28, Issue 29
Archive-Name: mrandom-3.0/part03

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 3 (of 6)."
# Contents:  doc/mrandom.txt src/mrandom.c
# Wrapped by [e-mail] on Fri May  6 21:42:55 1994
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'doc/mrandom.txt' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'doc/mrandom.txt'\"
else
echo shar: Extracting \"'doc/mrandom.txt'\" \(50066 characters\)
sed "s/^X//" >'doc/mrandom.txt' <<'END_OF_FILE'
X        Massachusetts Institute of Technology
X               77 Massachusetts Avenue
X                Cambridge, MA 02139
X
X
X              mrandom 3.0 User's Manual
X
X
X                by
X
X
X              Robert Plotkin
X
X
X
X    Department of Electrical Engineering and Computer Science
X
X                 May 1993
X
X
X1 Introduction
X--------------
X--------------
X
XThe mrandom package contains a library of routines for using random
Xnumber generators (RNGs) in C in the Unix 4.3bsd environment.  This
XUser's Manual is designed as a guide for programmers who wish to use the
Xmrandom package within their own programs.  The current version of
Xmrandom (version 3.0) is a major rewrite of mrandom version 2.1,
Xreleased by Clark Thomborson in July 1992 (available via anonymous
Xftp from theory.lcs.mit.edu).  The package now provides:
X
X* A standardized interface to many simultaneously-active RNGs,
X* A standardized, and unbiased, method for generating random
Xintegers in the range 0..m-1,
X* A standardized method for generating floating point numbers
Xuniformly distributed in [0.0, 1.0),
X* Two standardized methods for generating pseudorandom bit
Xstreams,
X* Time-efficient vectorized calls, returning multiple uniform
Xvariates,
X* Buffered and unbuffered calls for efficient generation of
Xpseudorandom generates in both large and small quantities,
X* The ability to ``split'' RNGs to produce parallel output streams
Xusing the ``leapfrog'' method,
X* A shorthand notation for completely specifying the
Xalgorithm and current state of an RNG, in an 80-character human-
Xreadable ASCII string,
X* A method for reconstructing an RNG state from its shorthand
Xnotation,
X* A standardized method for adding new RNGs to the package,
Xand
X* A file-I/O interface allowing fast saves and restarts of RNG state
Xvectors.
X
XYou can obtain the complete distribution of mrandom by anonymous ftp
Xfrom theory.lcs.mit.edu.  Take the file `mrandom.tar.Z' from the
Xdirectory /pub/cthombor/Mrandom/V3.0.  The distribution contains all
Xsource code, instructions, and this manual.
X
XIn addition, the mrandom package includes an unsupported set of routines
Xfor testing the statistical properties of the RNGs in the package.  The
Xroutines are packaged in an executable called mrtest.  Information about
Xmrtest is available in the doc/mrtest directory of the distribution.
XAlthough mrtest is included in the current distribution, it is not
Xsupported.
X
XQuestions, comments, bug reports, etc. should be sent to Clark
XThomborson at [e-mail].
X
X2 Files in the Distribution Directory
X-------------------------------------
X-------------------------------------
X
XThe mrandom source code distribution includes the following
Xfiles:
X
Xmakefile
X    The makefile for creating the mrtest program and the
X    mrandom.a library.
XREADME
X    General information about the mrandom package, including
X    changes to the last version.
Xmrandom.c mrandom.h
X    The source and header files for the main mrandom module.
Xbentley.c bentley.h
X    The source and header files for Bentley's version of the
X    generator described in Knuth Vol 2, Section 3.6.
Xpcrand.c pcrand.h
X    The source and header files for the Portable Combined RNG.
Xran0.c ran0.h ran1.c ran1.h ran2.c ran2.h
X    The source and header files for Press and Teukolsky's ran0,
X    ran1, and ran2.
Xultra.c ultra.h
X    The source and header files for Marsaglia's Ultra generator.
Xmrtest.c
X    The mrtest source file.
Xxsq.c xsq.h
X    Code used by mrtest.
Xrngs.h
X    The header file for the UNIX RNGs and the trivial RNG.
Xnewrng.c newrng.h
X    Source and header file templates for a new RNG.
Xmrandom.3
X    The man pages for mrandom.
Xmrtest.1
X    The man pages for mrtest.
Xscript
X    A test script for mrtest.
Xmrandom.tex
X    The latex source for this manual.
Xlatexinfo.sty
X    The style file needed to latex this manual.
Xmrandom.txt
X    Plain ASCII text version of this manual.
Xmrandom.ps
X    PostScript version of this manual.
X
X3 Installing mrandom
X--------------------
X--------------------
X
XPreparing the mrandom package for use by other programs is
Xsimple.  Merely position yourself in the directory which contains the
Xmrandom files and type:
X
Xmake all
X
XThis will compile all necessary source files and create the mrandom.a
Xlibrary, as well as the mrtest binary executable.
X
XYou can also make either the mrandom.a library or the mrtest executable
Xby typing:
X
Xmake mrandom.a
X
Xor
X
Xmake mrtest
X
Xrespectively.
X
XTo save disk space, various intermediate object files can be removed
Xwith:
X
Xmake clean
X
XThe system can be restored to its original state with:
X
Xmake realclean
X
XThe mrandom package is written in ANSI C, and should be
Xeasily portable to any UNIX-based system supporting a C language
Xcompiler.  However, when compiling on a new system, confirm that
Xthe long type is represented in 32 bits.
X
X4 What is an RNG?
X-----------------
X-----------------
X
XWithin the mrandom package, an RNG is represented by the RNGdata
Xstructure.  An RNG has several abstract characteristics which are
Xdescribed in this section.  For a more detailed description of the
Xrepresentation of the RNGdata structure in C, see
XAppendix A.
X
XAssociated with every RNG is:
X
X* an algorithm number, which determines which algorithm the RNG
Xuses to produce pseudorandom generates.  For descriptions of the
Xalgorithms which are currently installed in the package, see Appendix B.
X* an ``mrandom algorithm number,'' which determines which
Xalgorithm the mrandomrv routine will use to produce
Xrestricted-range integers when called with the RNG.  (See the
Xdescription of mrandomrv in Section 5.4.2 for more
Xinformation.)
X* a state vector, which contains the current state of the RNG.
X* the size of the state vector.
X* a seed vector, which contains the seeds which were used to
Xseed the RNG.
X* the size of the seed vector.
X* a count of the number of times the RNG has been called since
Xit was initialized.
X* two buffers
X    - a bit buffer for bit generates
X    - a main buffer for all other generates.
XSee Section 5.4.2 for a description of how the two buffers work.
X* a split value, which determines how many generates to skip
Xbetween values which are returned from the RNG.  For a detailed
Xdescription of the split value, see Section 5.4.7.
X* a string containing the human-readable name of the RNG.
X* the range of the RNG, such that the RNG is capable of
Xproducing integers with a maximum value of range-1.
X
X5 Using the mrandom package
X---------------------------
X---------------------------
X
X5.1 Overview of the Library
X
XThis section describes the mrandom library.  The procedures
Xin the library are gathered into the following groups:
X
X* functions for initializing and de-initializing (killing) RNGs
X* functions for returning generates from RNGs
X* functions for saving and restarting RNGs
X* a function for seeding RNGs
X* a function for checking the integrity of RNG state vectors
X* a function for producing a human-readable ASCII description
Xof an RNG
X* functions for examining and modifying characteristics of RNGs
X
X5.2 Return codes from mrandom routines
X
XMost mrandom routines follow one of two conventions for return
Xvalues.
X
XFor those routines which return pointers to an RNGdata
Xstructure:
X
X* A valid pointer to an RNGdata structure is returned upon success.
X* A null pointer is returned upon failure.
X
XFor routines which produce vectors of pseudorandom generates:
X
X* The first cell of the vector is returned upon success.
X* The program aborts upon error.
X
XFor all other routines:
X
X* A 0 or -1 is returned upon failure.
X
XDetails are provided in the description of the mrandom library
Xin Section 5.4.
X
X5.3 Linking
X
XIn order to use the mrandom library of routines in your programs
Xyou must:
X
X* Link the mrandom.a library with your program.  This will
Xtypically involve merely including the library on your compiler's
Xcommand line, such as:
X
Xcc myprog.c mrandom.a
X
XSince mrandom uses some UNIX mathematical functions, you may also need
Xto link a math library with your program, as in the following:
X
Xcc myprog.c mrandom.a -lm
X
XCheck the documentation for your C compiler for details.
X* Include the following line at the top of your source file:
X
X#include "mrandom.h"
X
X5.4 The mrandom library
X
X5.4.1 Initialization and De-Initialization
X
XIn order to use an RNG, it must first be initialized.  This is
Xaccomplished by first declaring a pointer to an RNGdata structure, and
Xthen calling init_rng, which allocates memory for the RNG and readies it
Xfor use by the other routines in the package.
X
X    RNGdata *init_rng(alg,mrandom_alg,seed,count1,count2,bufsize)
X    long alg;
X    long mrandom_alg;
X    long *seed;
X    long count1, count2;
X    long bufsize;
X
X    int kill_rng(rng)
X    RNGdata *rng;
X
Xinit_rng returns a pointer to an initialized RNG.  A pointer returned by
Xinit_rng is valid for use by all other mrandom routines.
X
Xalg is the number of the algorithm to be used by the RNG.  (See Appendix
XB.)
X
Xmrandom_alg is the algorithm to be used by mrandomrv when called with
Xthe RNG.  (See Section 5.4.2.)
X
Xseed is a pointer to a seed vector to be used to seed the RNG.
X(See Section 5.4.4.)
X
Xcount1 and count2 determine the number of initial values to be generated
Xby the RNG and then discarded, according to the formula:
X
Xnumber to discard = count1 + BILLION*count2,
X
Xwhere BILLION is defined in 'mrandom.h' as the decimal constant
X1000000000.
X
Xbufsize is the size of the RNG's main buffer.  A non-positive value of
Xbufsize will be interpreted as a value of 1.  (See Section 5.4.2 for
Xmore information on buffering.)
X
Xkill_rng destroys the RNG, making it invalid for use.  This procedure
Xde-allocates the space used by the RNG, and should therefore be used to
Xkill RNGs which will no longer be used.
X
XDo *not* use an RNGdata pointer which points to an active RNG
Xto store the return value of init_rng.  In order to initialize an
XRNG, you should either:
X
X* Declare a new RNGdata pointer, and then use it to store the
Xreturn value of init_rng, as shown in Figure 1.
X
X-----------------------------------------------
XRNGdata *rng;
Xlong seed[1];
X
Xrng=init_rng(2, 0, seed, 10000, 0, 8192);
X
XFigure 1: Proper initialization of an RNG
X-----------------------------------------------
X
X* Use kill_rng to de-initialize an RNGdata pointer which points to an
Xactive RNG, and then use that pointer to store to the return value of
Xinit_rng, as shown in Figure 2.  Figure 3 shows how an RNG should *not
Xbe re-initialized.
X
X-----------------------------------------------
XRNGdata *myrng;
Xlong seed[1];
X
Xseed[0]=12345;
Xmyrng=init_rng(2, 0, seed, 10000, 0, 8192);
Xkill_rng(myrng);
Xmyrng=init_rng(3, 0, seed, 5000, 0, 2048);
X
XFigure 2: Proper re-initialization of an RNG
X-----------------------------------------------
X
X-----------------------------------------------
XRNGdata *myrng;
Xlong seed[1];
X
Xseed[0]=12345;
Xmyrng=init_rng(2, 0, seed, 10000, 0, 8192);
Xmyrng=init_rng(3, 0, seed, 5000, 0, 2048);
X
XFigure 3: Improper re-initialization of an RNG
X-----------------------------------------------
X
X5.4.2 Procedures for Generating Pseudorandom Numbers
X
XAn initialized RNG can be used to generate pseudorandom numbers by using
Xa variety of routines described in this section.
X
XRoutines are provided for producing generates of four types:
X
X* double precision floating point generates in the range [0,1)
X* single precision floating point generates in the range [0,1)
X* long integer generates in the range 0..r-1, where r is the range of
Xthe RNG being used to produce generates
X* long integer generates in the range 0..m-1, for any
X1 <= m < range_rng(rng)
X
XNote that although both single and double precision floats can be
Xreturned by the generate-producing routines, the actual precision of the
Xgenerates produced is determined by the precision of the underlying
Xgenerator being used.  In other words, the difference between routines
Xwhich return generates of type double and those which return generates
Xof type float is merely in the ``packaging'' of the generates, not in
Xthe precision they provide.  Information about the precision of the RNGs
Xcurrently installed in the package is in Appendix B; such information
Xcan also be obtained at run-time through the range_rng routine (see
XSection 5.4.7).  The current version of the package only supports RNGs
Xwith precisions of no more than 32 bits.
X
XBoth buffered and unbuffered routines are provided.  Unbuffered routines
Xcall the underlying RNG only as many times as are needed to produce the
Xrequested number of generates, while buffered routines maintain buffers
Xof generates, so that generates may be produced efficiently even when
Xrequested in small quantities.  Roughly, buffered routines are
Xpreferable when generates are requested one at a time or in small
Xquantities, while unbuffered routines are preferable when generates are
Xrequested in large quantities.  Some other differences between buffered
Xand unbuffered routines are discussed later in this section.  The size
Xof the buffer used by an RNG is determined at the time of the RNG's
Xinitialization; effective buffer sizes will vary from application to
Xapplication.
X
XThe name of a routine denotes the type of the value which the routine
Xreturns and whether the routine is buffered or unbuffered. The first
Xletter of a routine denotes the type of value which it returns: ``d''
Xfor double precision and ``f'' for single precision floating point in
Xthe range [0,1); ``l'' for long integer in the range
X0..(range_rng(rng)-1), and ``b'' for bit (either a 0 or a 1).  If the
Xsecond letter of the routine's name is an ``x'', then the routine is
Xunbuffered.  Otherwise, the routine is buffered.
X
XFor convenience in user programming, we also provide a number of macros
Xthat supply default parameter values.  The last two letters of all our
Xfundamental routines is ``rv''.  This means that they must be provided
Xwith both a pointer to an RNGdata structure and a vector to fill with
Xgenerates from the RNG.  Macros whose names do not contain an ``r'' have
Xthe RNGdata pointer omitted from their parameter list; they use the
Xmost-recently initialized or restarted RNG to produce generates.  Macros
Xwhose names do not contain a ``v'' have the vector and number of
Xgenerates omitted from their parameter list; they produce and return a
Xsingle generate.
X
XAll generating routines abort with a message to stderr if called with an
Xinvalid RNGdata pointer.
X
XBuffering
X
XThe operation of the buffered routines and their interaction with the
Xunbuffered routines requires some elaboration.  Each RNG maintains two
Xseparate buffers: one for buffering bit values (the ``bit buffer''), and
Xone for buffering all other values (the ``main buffer'').  The size of
Xthe main buffer of an RNG is determined at the time of the RNG's
Xinitialization, while the size of the bit buffer is currently fixed at
X32 bits.
X
XConsider a freshly-initialized RNG with a main buffer size of 1000.  A
Xrequest is made for a single generate of type long.  The RNG's buffer
Xgets filled with 1000 generates, and the first such generate is returned
Xto the user.  So the buffer now contains 999 generates.  If another
Xgenerate of type long, double, or float, is requested, a generate will
Xbe pulled from the buffer and returned to the user after being converted
Xto the proper type.  If the user continues to request generates in this
Xway and the main buffer becomes depleted, it will be filled again with
X1000 generates, and so on.
X
XThe unbuffered routines do not interfere with either of the RNG's
Xbuffers.  Again, consider our RNG, with its buffer filled with 1000
Xgenerates.  The user now makes a request for a single unbuffered
Xgenerate.  The underlying RNG will then be called once, returning a
Xsingle generate, leaving our buffer of 1000 generates untouched, and
Xstill ready to be accessed by the buffered routines.  If, in a
Xparticular application, it is necessary to always use consecutive
Xgenerates from an RNG, then that RNG should always be called *either*
Xwith buffered or unbuffered routines, but not with a combination of
Xboth.
X
XThe bit buffer of an RNG operates similarly to the main buffer, with one
Xkey difference: the bit buffer is filled by sequentially retrieving
Xgenerates from the main buffer.  Once again, consider our RNG, with its
Xbuffer filled with 1000 generates, and with its bit buffer empty.  A
Xsingle bit is then requested.  Thirty-two generates will be pulled from
Xthe main buffer, transformed into thirty-two one-bit 0-1 values, then
Xstored in the bit buffer.  (For users who are more concerned with speed
Xthan accuracy, we also provide a ``fast'' bit-buffer call, in which a
Xsingle 32-bit generate from the main buffer is transformed into
Xthirty-two 0-1 variates.  See the descriptions of bxrandom_f and
Xbrandom_f, below.)
X
XUnbuffered and buffered calling procedures
X
X    double dxrandomrv(rng, n, v)
X    double drandomrv(rng, n, v)
X    RNGdata *rng;
X    long n;
X    double v[];
X
X    float fxrandomrv(rng, n, v)
X    float frandomrv(rng, n, v)
X    RNGdata *rng;
X    long n;
X    float v[];
X
X    long lxrandomrv(rng, n, v)
X    long lrandomrv(rng, n, v)
X    RNGdata *rng;
X    long n;
X    long v[];
X
X    int bxrandomrv(rng, n, v)
X    int brandomrv(rng, n, v)
X    RNGdata *rng;
X    long n;
X    int v[];
X
X    int bxrandomrv_f(rng, n, v)
X    int brandomrv_f(rng, n, v)
X    RNGdata *rng;
X    long n;
X    double v[];
X
XThese routines fill the vector v with n generates from rng, and return
Xthe first generate produced (i.e. v[0]).
X
XIf rng is a null pointer, then the most-recently initialized or
Xrestarted RNG is used to produce generates.  If n is 0, then the v
Xparameter need not be provided, and a single generate is produced and
Xreturned.
X
Xbxrandomrv uses one generate from rng to generate each bit.  This
Xroutine is slower than bxrandomrv_f, but returns bits of higher quality.
X
Xbxrandomrv_f uses each generate from rng to produce 32 bits.  Therefore,
Xrequests for bits in other than multiples of 32 will result in bits from
Xthe stream being ``lost'' between calls.  The routine returns -1 if it
Xis called with an RNG whose range is not 2^32.  This routine is faster
Xthan bxrandomrv, but we do not recommend its use, since we know of no
Xone who has rigorously tested such a bit stream.  We gain confidence in
Xour slower bxrandomrv bit stream, in comparison, every time the
Xunderlying generator passes a test sensitive to correlations in the
Xleading digit of its floating point output, or to the most significant
Xbit of its fixed point output.
X
Xbrandomrv_f, the buffered version of bxrandomrv_f, does not exhibit the
Xsame property of ``losing bits'' as does bxrandomrv_f, since bits which
Xare not used in one call to brandomrv_f are stored in the bit buffer and
Xare available for use upon future calls.
X
X    int flush_rng(rng)
X    RNGdata *rng;
X
X    Flushes both of the RNG's buffers.
X
XProcedure for generating integers in a restricted range
X
X    long mrandomrv(rng, m, n, v)
X    RNGdata *rng;
X    long m,n,v[];
X
Xmrandomrv fills the vector v with n generates in the range 0..m-1 using
Xrng, where 1 <= m <= range_rng(rng).  If range_rng(rng) < m, the program
Xaborts with an error.
X
XThe algorithm used by mrandomrv to fill v is set by init_rng or by
Xmralg_rng.  (See Section 5.4.1 and Section 5.4.7.)
X
XAlgorithm 0 is Thomborson's unbiased method, which produces unbiased
Xlong integers in the range [0..m).  The algorithm discards any outputs
Xfrom rng which are larger than r-(r mod m), where r is equal to
Xrange_rng(rng).  At worst, this code will discard (on long-term average)
Xat most one value of r for every one that is useful.  This worst case is
Xonly encountered for extremely large m; for fixed and moderate m, this
Xcode will rarely discard a value, and thus will run essentially as fast
Xas algorithm 1.  When the value of m changes on each call to mrandom,
Xhowever, this code is slower than algorithm 1, due to the necessity of
Xrecomputing r-(r mod m).
X
XThe program aborts with an error message to stderr if rng is behaving so
Xnon-randomly that Algorithm 0 must make an excessive number of calls to
Xrng in order to produce the requested number of generates.
X
XThe program aborts with an error to stderr if mrandomrv is asked to use
XAlgorithm 0 with a value of m for which m > range_rng(rng).
X
X
XAlgorithm 1 is the standard (long)(m*dxrandomr(rng)).  This algorithm
Xmay be biased: for large m, some results may be be more likely than
Xothers.  The bias is (r mod m)/m, which is upper-bounded by 0.1% if m is
Xless than a million and the range r of rng is at least a billion.
X
XWe do not support, and indeed we actively discourage, generating
Xrestricted-range integers with lrandomr(rng)%m.  Many RNGs have poor
Xbehavior under this transformation, most noticeably when m is a power of
X2.  When m is not a power of 2, fixed-point division required by an
X``%'' operation is time-consuming on many workstations.
X
XNOTES
XThe mrandomrv procedure is capable of generating long integers in the
Xfull range of any RNG for which 1 <= range_rng(rng) <= 2^32.  In order
Xto accomplish this, with the parameter m a signed long integer, the
Xfollowing mapping is used:
X
XRange(mrandom(m)) = 0..m-1        if 1 <= m < 2^31
X            0.. 2^32-1      if m=0
X            0..(2^31-m-1)    if -2^31 <= m < 0
X
XMacros are defined for easy calling of mrandomrv with various default
Xparameters.  See Section 5.4.2 for the naming conventions followed by
Xthe macros.
X
X5.4.3 Saving and restarting RNGs
X
XRNGs may be saved to disk files in a human-readable ASCII format and
Xlater restarted.  RNG buffers are not saved, and therefore all restarted
XRNGs have empty buffers, and any data remaining in an RNG's buffer at
Xthe time of its state-save will *not* be reconstructed.
X
X    int save_rng(rng, filename)
X    RNGdata *rng;
X    char *filename;
X
X    RNGdata *restart_rng(filename)
X    char *filename;
X
Xsave_rng saves rng to the ASCII file named filename.
X
Xrestart_rng restarts an RNG from a previously saved statefile.  The
Xrestarted RNG will begin where the saved RNG ``left off.''  As with
Xinit_rng, the RNGdata pointer used to store the restarted RNG *must* be
Xeither a freshly declared pointer or a pointer to a freshly killed RNG
X(see Section 5.4.1).
X
XRNGs may store their state and seed vectors in any of a number of
Xformats, and this is reflected in the format of the state file.
XFigure 4 shows a sample state file of an RNG using the Knuth/Bentley
Xlagged-Fibonnacci generator prand (see Appendix B), which stores its
Xstate and seeds as 32-bit long ints.  Figure 5 shows a sample state
Xfile of an RNG using 4.3bsd nrand48, which stores its state and seeds
Xas 16-bit ints.
X
X--------------------------------------
XRNG statefile for algorithm 2, (Knuth/Bentley prand: lagged Fibbonacci)
XBuffer size = 1024 bytes
XInitial seed table =
X   00000927
XNumber of calls to underlying RNG after seeding = 0 billion + 2000
XNext value in this pseudorandom sequence = 0b64d0ea
XThis RNG returns every 1 generates
XThis RNG uses mrandom algorithm 0
XRNG state table =
X   0911c27a 10641ca0 2ba00807 1aabed0a
X   273ff367 1ab88564 2ae76a9e 2a7e6bc0
X   35c7568e 201b6b04 3ad90695 303208b2
X   1e718896 054c9886 00e8c93f 130a41cb
X   11de97bf 0da54e15 2f4fcca0 0ebb1f70
X   01c195c3 3283980e 37dee108 0893a89b
X   326849b0 167bb45e 19cc9765 33d97b51
X   36b425d1 35704e34 29a638ca 280a086f
X   11dfa5d6 14dcbcc4 2610bdf4 02534109
X   2817daf4 0bcf76ab 19b0a07d 0eebf7f6
X   113c003e 31b996b0 12bab234 05eddb36
X   1ed71381 377742a3 3878e079 2668c922
X   22cc8033 22368c85 18e960ea 2002b06f
X   22ff23e8 251187dc 340c3dcd 00000023
X   00000004
X
XFigure 4: A sample RNG state file for the Knuth/Bentley prand().
X--------------------------------------
X
X--------------------------------------
XRNG statefile for algorithm 4, (4.3bsd nrand48.c: 48-bit multiplicative)
XBuffer size = 8192 bytes
XInitial seed table =
X   0096   b43f   0034   bf15
X
XNumber of calls to underlying RNG after seeding = 0 billion + 11000
XNext value in this pseudorandom sequence = 04a3689e
XThis RNG returns every 1 generates
XThis RNG uses mrandom algorithm 0
XRNG state table =
X   07c5   8f2d   0000   a7d6
X
XFigure 5: A sample RNG state table for nrand48
X--------------------------------------
X
XA few examples of how to save and restart RNGs are displayed in Figure
X6.
X
X--------------------------------------
X/* Proper way to re-initialize an active RNG */
Xmrandom(rng,10,n,v);
Xkill_rng(rng);
Xrng=restart_rng("mystatefile");
X
X/* Proper way to restart an inactive RNG */
XRNGdata *rng;
Xrng=restart_rng("mystatefile");
X
X/* Improper way to restart an active RNG */
Xmrandom(rng,10,n,v);
Xrng=restart_rng("mystatefile");
X
XFigure 6: Examples of saving and restarting RNGs
X--------------------------------------
X
X5.4.4 Seeding
X
XEach RNG is initially seeded during initialization by init_rng.  An RNG
Xmay also be reseeded at any time after initialization.
X
X    void seed_rng(rng,seed)
X     RNGdata *rng;
X     long *seed;
X
Xseed_rng seeds rng with the seed table pointed to by seed.  The RNG's
Xcounter is reset to 0 and its buffers are flushed.
X
X5.4.5 Checking RNG integrity
X
XAn RNG can be checked to see if it has been corrupted or is otherwise
Xnot in proper condition for use.
X
X    int check_rng(rng);
X    RNGdata *rng;
X
Xcheck_rng checks the integrity of the RNG, in order to determine whether
Xit can be used by the other mrandom library routines.
X
X5.4.6 Obtaining a human-readable description of the RNG
X
X    char *describe_rng(rng,rngid)
X    RNGdata *rng;
X    char rngid[RNGIDSTRLEN];
X
Xdescribe_rng places a human-readable description of rng in
Xthe string rngid.  The string has the following format:
X
XRNG state identifier is (alg, mralg: seed1, seed2; count1,count2;
Xbufsize, split)
X
Xwhere
X
X* alg is the number of the algorithm used by rng to generate
Xpseudorandom numbers.  (See Appendix B.)
X* mralg is the number of the algorithm used by
Xmrandomrv when called with rng.  (See Section 5.4.2.)
X* seed1 and seed2 are the first and second entries in trng's seed table.
XIf rng's seed table has more than two entries, only the first two are
Xincluded in its description.  (See Section 5.4.4.)
X* count1 and count2 represent rng's counter.  (See Section 5.4.1.)
X* bufsize is the number of entries in rng's buffer.  (See Section
X5.4.2.)
X* split is rng's current split value.  (See Section 5.4.7.)
X
Xdescribe_rng exits with a message to stderr if called with an invalid
XRNGdata pointer.
X
X5.4.7 Procedures for examining and modifying RNG parameters
X
XProcedures are available for examining and modifying an RNG's parameters
Xonce it has been initialized.
X
X    int mralg_rng(rng, new_value)
X    RNGdata *rng;
X    long new_value;
X
X    int split_rng(rng, new_value)
X    RNGdata *rng;
X    long new_value;
X
X    double range_rng(rng)
X    RNGdata *rng;
X
Xmralg_rng sets rng's mrandom algorithm number (See Section 5.4.2 for
Xinformation on mrandom algorithm numbers).  It returns 0 if new_value is
Xan invalid value.
X
Xsplit_rng sets the split value of rng.  It returns 0 if new_value < 0.
XAn RNG's split value is set to SPLIT_DEF upon initialization.  SPLIT_DEF
Xis #defined in 'mrandom.h', and currently has a value of 0.
X
XThe function of the split value is to simulate one ``branch'' of a
Xgenerator which has been ``split'' into two or more generators.  This is
Xbest illustrated with an example.  Consider an (apparently non-random)
XRNG which returns the raw sequence:
X
X0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...
X
XThe split value indicates how many elements of the sequence to skip
Xbetween generates.  For example, if our sample RNG were given a split
Xvalue of 1 immediately after initialization, it would then return the following
Xsequence:
X
X0 2 4 6 8 10 ...
X
XA split value of 2 after initialization would produce:
X
X0 3 6 9 12 15 ...
X
XAn RNG may be split at any time after its initialization.  So, for
Xexample, our sample RNG might be initialized and then made to generate
Xthe following values:
X
X0 1 2 3 4 5
X
Xbefore being split with a split value of 3, producing the following
Xgenerates:
X
X6 10 14 18 22 ...
X
XSplitting can be used to create several ``leapfrogged'' RNGs from one
XRNG, as shown in Figure 7.
X
X------------------------------------
XRNGdata *rngs[10];
Xlong seed;
Xint i;
X
Xseed=12345;
Xfor (i=0; i<10; i++) {
X  /* RNG #i gets cycled i times */
X  rng[i]=init_rng(2,0,&seed,i,0,1024);
X  split(rng[i],9);
X}
X
XFigure 7: Creating ``leapfrogged'' RNGs.
X------------------------------------
X
XThis operation may be useful in parallel codes, as in testing an RNG for
Xlong-range correlations.  Unfortunately, our current implementation is
Xinefficient for leapfrogging large numbers of RNGs.  A more efficient
Xmethod may be included in a future version of mrandom.
X
XThe split value affects all of the pseudorandom number generating
Xroutines (See Section 5.4.2).
X
Xrange_rng returns the range of rng.
X
XNOTES
XThese procedures exit with a message to stderr if rng does not
Xpoint to a valid RNGdata structure.
X
X6 Adding new RNGs to the Package
X--------------------------------
X--------------------------------
X
XThis section is designed for the programmer who wishes to add new RNGs
Xto the mrandom package.  The section first describes the routines which
Xmust be provided by the programmer to serve as an interface between the
XRNG and the mrandom package.  It then describes the RNG parameters which
Xmust be defined in a header file for the RNG, and how these parameters
Xare to be incorporated into the code for mrandom itself.  Finally, it
Xdescribes how to remake the mrandom package after adding or modifying
XRNGs.  Appendix B contains descriptions of the ten RNGs installed in the
Xcurrent version of the package.
X
X6.1 Routines Provided by the Programmer
X
XThis section describes routines which must be provided by each RNG in
Xthe package.
X
XThe routines provided by the programmer must manipulate an
XRNGdata structure.  In order to facilitate this, two macros are
Xavailable for accessing an RNG's state and seed vectors:
X
X* RNGstate refers to the RNG's state vector.
X* RNGseed refers to the RNG's seed vector.
X
XThese vectors are one-dimensional C arrays, e.g.  RNGstate[0] is the
Xfirst element in the RNG's state vector, RNGstate[1] is the second
Xelement, etc.  In order to make use of these macros, the name of the
XRNGdata pointer in your routines' parameter lists *must be rng, as shown
Xin the examples in this section.
X
XThe names used as examples in this section begin with ``myrng'';
Xhowever, there are no restrictions on naming of routines provided by an
XRNG.  However, for ease of readability and consistency, we suggest that
Xthe naming conventions used in this section be followed.
X
XThe routines described in this section should be included in a single .c
Xfile.  A template for such a source file, called `newrng.c', is included
Xin the distribution and displayed in Figure 8.
X
XRemember that all RNG state information must be included in the RNGstate
Xfield of the RNGdata structure.  In particular, do *not use global or
Xstatic variables to hold RNG state information.  Doing so will make it
Ximpossible to run several instantiations of your RNG simulataneously.
X
X------------------------------------
X/* newrng.c */
X/* RNG source file template */
X/* Robert Plotkin
X/* 5/3/93 */
X
X#include "newrng.h"
X
X/* Generating procedure */
X/* Only one of the following two procedures should be */
X/* defined, depending on the kind of value that */
X/* your RNG returns */
X
Xlong newrng_lgen(rng)
XRNGdata *rng;
X{
X/* Your generating procedure goes here */
X}
X
Xdouble newrng_dgen(rng)
XRNGdata *rng;
X{
X/* Your generating procedure goes here */
X}
X
X/* Seeding procedure */
Xvoid newrng_seed(rng,seed)
XRNGdata *rng;
Xlong *seed;
X{
X/* Your seeding procedure goes here */
X}
X
X/* Checking procedure */
Xint newrng_check(rng)
XRNGdata *rng;
X{
X/* Your checking procedure goes here */
X}
X
XFigure 8: RNG source file template
X------------------------------------
X
X6.1.3 Seeding Routine
X
Xvoid myrng_seed(rng, seed)
XRNGdata *rng;
Xlong *seed;
X
XThis procedure is used for seeding the RNG.  The interpretation of the
Xseed parameter is left entirely to the programmer.  It may, for example,
Xpoint to a single integer or to an array of 5,000 integers.  One of the
XRNGs currently installed in the package interprets the seed parameter as
Xpointing to three 16-bit integers.  Obviously, RNGs which are capable of
Xbeing seeded with a variable number of seeds need to be passed a seed
Xpointer which contains adequate information about the number of seeds to
Xwhich it points.
X
XAlthough the seeding procedure is passed an entire RNGdata
Xstructure as a parameter, it should only manipulate the RNGstate
Xfield of that structure.  (See Appendix A for information
Xon the RNGdata structure.)  Many RNG seeding procedures will simply
Xcopy the seed parameter into RNGstate, as shown in Figure 9.
X
X------------------------------------
Xvoid myrng_seed(rng, seed)
XRNGdata *rng;
Xlong *seed;
X{
XRNGstate[0]=seed[0];
XRNGstate[1]=seed[1]; /* This RNG uses two long seeds */
X}
X
XFigure 9: A sample seeding procedure
X------------------------------------
X
XOther seeding procedures may fill RNGstate with the results
Xof some complicated function performed on the initial seed table.
X
X6.1.3 Pseudorandom Number Generating Procedure
X
Xlong myrng_lgen(rng)
XRNGdata *rng;
X
Xdouble myrng_dgen(rng)
XRNGdata *rng;
X
XThe programmer must provide *one* procedure, matching one of the two
Xprototypes given above, which returns a single generate from the RNG.
XThe routine may return either:
X
X* a double precision floating point number in the range [0,1), or
X* a long (32-bit) integer in the range 0..range_rng(rng)-1
X
XIt is pointless for the programmer to provide procedures of both types.
XIf this is done, only one of them will be accessible by any user code,
Xdepending on the value given to RNGreturns, RNGdgen, and RNGlgen (see
XSection 6.2).
X
XAlthough the generating procedure is passed an entire RNGdata
Xstructure, it should only manipulate the RNGstate field of the
Xstructure.  A sample generating procedure is displayed in Figure 10.
X
X------------------------------------
Xlong myrng_lgen(rng)
XRNGdata *rng;
X{
XRNGstate[0]*=12345+6789;
Xreturn(RNGstate[0]);
X}
X
XFigure 10: A sample generating procedure
X------------------------------------
X
X6.1.4 RNG Checking Procedure
X
Xint myrng_check(rng)
XRNGdata *rng;
X
XThe programmer must provide a procedure to check the integrity of the
XRNGdata structure.  The procedure returns a value of 1 if the RNG is fit
Xfor use, and returns 0 otherwise.  The coding of the procedure is
Xentirely RNG-specific, and may be extremely simple or extremely
Xcomplicated, depending on the nature of the RNG and the extent of
Xintegrity desired.  On one extreme is the procedure which always
Xdeclares success, and on the other extreme is the perfect (and slow)
Xprocedure which creates a new RNG, seeds it with the seeds of the RNG to
Xbe checked, cycles it through the number of generates which were
Xproduced by the RNG to be checked, and compares the state tables of the
Xtwo RNGs.  Clearly, the procedure should not modify the RNG in any way.
XWhen writing a checking procedure it might be useful to examine those
Xincluded in the existing package.
X
X6.2 RNG Header Files
X
XFor each RNG included in the package, there must be a corresponding
Xheader file.  The header file contains information about the RNG which
Xis used by the mrandom library routines.  This section describes the
Xinformation contained in RNG header files, and describes how to use such
Xheader files to incorporate new RNGs into the mrandom package.  A
Xtemplate for such a header file, called `newrng.h', is included in
Xthe distribution and displayed in Figure 11.
X
X------------------------------------
X/* newrng.h */
X/* RNG header file template */
X/* Robert Plotkin
X/* 5/3/93 */
X
X#include "mrandom.h"
X
X/* Information for mrandom */
X#define RNGstatesize_n
X#define RNGseedsize_2
X#define RNGrange_2
X#define RNGname_2
X#define RNGreturns_2
X#define RNGstatetype_2
X#define RNGdgen_2
X#define RNGlgen_2
X#define RNGseed_2
X#define RNGcheck_2
X
X/* mrandom interface routines */
Xlong newrng_gen(/* RNGdata * */);
Xvoid newrng_seed(/* RNGdata *, long * */);
Xint newrng_check(/* RNGdata * */);
X
XFigure 11: RNG header file template
X------------------------------------
X
X6.2.1 Information in Header Files
X
XEach header file should begin with the following directives:
X
X#ifndef MRANDOM
X#include "mrandom.h"
X#endif
X
XDefinition of RNG parameters
X
XThe next set of lines in the file should contain #define statements
Xwhich assign values to the RNG's parameters.  The names used in the
X#define statements must contain the RNG's number.  There are currently
Xten RNGs included in the package, labeled 0 through 9.  The next RNG
Xincluded in the package should be labeled number 10, and so on.  The
X``n'' in each parameter name in the following list should be interpreted
Xas the RNG's number.
X
XRNGname_n
X    A string constant containing the name of the RNG, terminated
X    with a newline character.
XRNGstatesize_n
X    The number of entries in the RNG's state table.  Each entry is a
X    (32-bit) long.  If the RNG is capable of using state tables of
X    varying sizes, RNGstatesize_n should be defined as the maximum
X    possible size.
XRNGseedsize_n
X    The number of entries in the RNG's seed table.  Each entry is a
X    (32-bit) long.  If the RNG is capable of using seed tables of
X    varying sizes, RNGseedsize_n should be defined as the maximum
X    possible size.
XRNGrange_n
X    The range of the RNG, expressed as a double precision floating
X    point number.  The range of the RNG is one more than the maximum
X    value the RNG is capable of generating.  For RNGs which produce
X    double precision generates with a precision of p (i.e. in the
X    range [0,(RNGrange-1.0)/(1<<p)), RNGrange should be defined as
X    2^p.  For example, an RNG which produces 8-byte IEEE
X    floating-point generates using single-precision IEEE arithmetic
X    (24-bit mantissas) has a range of 16777216.0.
X
XRNGreturns_n
X    A number signifying the type of the generate returned by the
X    RNG.  An RNG can return a value of one of two types:
X        * a long in the range 0..RNGrange-1
X        * a double in the range [0,1)
X    RNGs which return values of type long and double return
X    types RET_LONG and RET_DOUBLE, respectively, as defined in
X    mrandom.h.
XRNGstatetype
X    A number signifying the interpretation of the values stored in
X    the RNG's state and seed vectors.  This value is used by the
X    routines that read and write the ASCII state files, thereby
X    allowing portability of state files across machines with
X    different byte orderings (see Section 5.4.3).  The following
X    values are currently supported:
X
X            Value        Type
X            ---------------------------
X            STATE_CHAR    8-bit character
X            STATE_INT    16-bit integer
X            STATE_LONG    32-bit long integer
X
X    The values of STATE_FLOAT (IEEE-standard 32-bit float) and
X    STATE_DOUBLE (IEEE-standard 64-bit float) are not currently
X    supported and are reserved for future use.
XRNGdgen_n and RNGlgen_n
X    The label of the procedure to be used for generating
X    pseudorandom numbers.  If the RNG returns doubles, then RNG_dgen
X    should be defined as the label of the RNG generating procedure,
X    and RNG_lgen should be defined as 0.  If the RNG returns longs,
X    then RNG_lgen should be defined as the label of the RNG
X    generating procedure, and RNG_dgen should be defined as 0.
XRNGseed_n
X    The label of the procedure to be used for seeding the RNG.
XRNGcheck_n
X    The label of the procedure to be used for checking the integrity
Xof the RNG.
X
XProcedure prototypes
X
XFinally, the header file must contain function prototypes for
Xthe three procedures provided by the RNG, so that the procedures
Xcan be accessed by the main mrandom code.  For example:
X
Xlong myrng_gen();
Xvoid myrng_seed();
Xint myrng_check();
X
X6.3 Modifying the mrandom code
X
XOnly a few lines of mrandom.h and mrandom.c need to be modified when
Xadding a new RNG to the package.
X
X* The number of RNGs currently installed in the package is defined as
XNUM_RNGS in `mrandom.h'.  The current value is 10.  This value should be
Xincremented when a new RNG is added to the package.
X* The header file for the new RNG needs to be #included in mrandom.c.
XThe #include directive should be included in the section marked by the
Xcomment ``Header files for RNGs currently included in package.''
X* Several additions need to be made in mrandom.c in the section marked
Xby the comment ``Arrays to hold information about RNGs.''  This section
Xof the code declares and initializes several arrays which hold
Xinformation about the RNGs included in the package.  When installing a
Xnew RNG, the appropriate #defined values need to be inserted at the end
Xof each initialization list.  For example, the declaration of RNGname_a
Xcurrently reads:
X
Xchar RNGname_a[NUM_RNGS][RNGIDSTRLEN]={RNGname_0, RNGname_1,
X      RNGname_2, RNGname_3, RNGname_4, RNGname_5, RNGname_6,
X      RNGname_7, RNGname_8, RNGname_9};
X
XAfter adding a new RNG to the package, this declaration would read:
X
Xchar RNGname_a[NUM_RNGS][RNGIDSTRLEN]={RNGname_0, RNGname_1,
X      RNGname_2, RNGname_3, RNGname_4, RNGname_5, RNGname_6,
X      RNGname_7, RNGname_8, RNGname_9,
X      /* RNG #10 added -> */ RNGname_10};
X
XThe arrays statesize_a, seedsize_a, range_a, returns_a, statetype_a,
Xseed_a, dgen_a, lgen_a, and check_a need to be similarly modified.
X
X6.4 Remaking the mrandom Package
X
XOnce you have added an RNG to the package as described in the
Xprevious sections, you will need to remake the mrandom package.  To do
Xthis:
X
X
X* Make sure that all of the files for the mrandom package, include the
Xsource and header files for your new RNG, are in the same directory.
X* Include the names of your header, source, and object files in makefile
Xon the lines labeled INCS, SRCS, and OBJS, respectively, as show in
XFigure 12.
X* Follow the instructions for making the mrandom package, as
Xdescribed in Section 3.
XOnce the package has been remade it will be ready for use, with your new
XRNG, by other programs.
X
X--------------------------------
XINCS = mrandom.h bentley.h pcrand.h ran0.h ran1.h ran2.h ultra.h xsq.h myrng.h
XSRCS = mrtest.c mrandom.c bentley.c pcrand.c ran0.c ran1.c ran2.c ultra.c xsq.c myrng.c
XOBJS = mrandom.o bentley.o pcrand.o ran0.o ran1.o ran2.o ultra.o xsq.o myrng.o
X
XFigure 12: Addition of myrng to makefile
X--------------------------------
X
X
XA The RNGdata Structure
X-----------------------
X-----------------------
X
XA.1 Introduction
X
XThis section describes the representation in C of the RNGdata structure
Xwhich is used by the mrandom package to represent RNGs.  This structure
Xneed never be manipulated by the programmer, except as described in
XSection 6.1.  This section, therefore, is intended for those who are
Xinterested in learning a little more about the inner workings of the
Xmrandom package.
X
XIn order to generate random numbers, the user must first declare a
Xpointer to an RNGdata structure, and use init_rng to allocate space for
Xthe RNG and perform various initialization functions.  The user uses the
XRNG entirely through calls provided by the interface described in
XSection 5.4; i.e. the user should not directly manipulate the RNGdata
Xstructure.
X
XA.2 Inside the Structure
X    The definition of the RNGdata structure is displayed in
XFigure 13.
X
X--------------------------------
Xstruct rngdata {
X    long rngalg;
X    long mrandom_alg;
X    long *rngstate;
X    long *rngseed;
X    long rngcount1;
X    long rngcount2;
X    struct {
X        long size;
X        long nleft;
X        long nbleft;
X        double *dbuf,*dbufp;
X        long *lbuf,*lbufp;
X        int *bbuf,*bbufp;
X    } buffer;
X    long rngnextval;
X    long rngsplit;
X    char rngname[];
X    long rngstatesize;
X    long rngseedsize;
X    long rngrangem1;
X    double rngrange;
X    signed int rngreturns;
X};
Xtypedef struct rngdata RNGdata;
X
XFigure 13: The RNGdata structure
X--------------------------------
X
XDescriptions of its fields are as follows:
X
Xrngalg
X    A number identifying the algorithm to be used by the RNG to
X    produce pseudorandom generates.  Algorithms in the package are
X    numbered sequentially starting with 0; currently there are 10
X    algorithms installed, numbered 0 through 9.  A table of RNGs
X    which are currently installed in the mrandom package, with their
X    corresponding algorithm numbers, is in Appendix B.
Xmrandom_alg
X    The algorithm use by mrandomrv when called with this RNG.
X    See Section 5.4.2 for more on mrandomrv.
Xrngstate
X    A pointer to the RNG's state vector, used to store the current
X    state of the RNG.  See Sections 5.4.3, 6.1, and 6.2.1 for more
X    information on RNG state vectors.
Xrngseed
X    A pointer to the RNG's seed vector.  See Section 5.4.4
X    for more information on RNG seed vectors.
Xrngcount1, rngcount2
X    These two values represent the number of generates the RNG has
X    produced since initialization, according to the formula:
X        rngcount1+rngcount2*BILLION
X
X    where BILLION is defined in `mrandom.h'.  Please note that the
X    value represented by rngcount1 and rngcount2 is the *total*
X    number of generates produced by the RNG since initialization,
X    including those discarded due to splitting of the RNG.  (See
X    Section 5.4.7 for more information about splitting RNGs.)
Xrngnextval
X    The next value to be output from the RNG.  This
X    value is used internally by the mrandom library and is not
X    guaranteed to be accurate.
Xrngsplit
X    Every (split+1)-th generate of the underlying RNG will be
X    returned by the RNG calling procedures.  rngsplit is set to
X    DEF_SPLIT upon initialization of the RNG, as defined in
X    `mrandom.h'.  See Section 5.4.7 for more information about
X    splitting RNGs.
Xbuffer
X    This structure contains information about the RNG's buffer and
X    its bit buffer.  (See Section 5.4.2 for more information on RNG
X    buffers.)  It contains several fields:
X
X    size
X        The number of entries in the RNG's buffer.
X    nleft
X        The number of values left in the RNG's buffer.
X    nbleft
X        The number of values left in the RNG's bit buffer.
X    dbuf, dbufp
X        A pointer to the first entry in the double
X        buffer, and a pointer to the next entry to be retrieved
X        from the double buffer.
X    lbuf, lbufp
X        Same for the long buffer.
X    bbuf, bbufp
X        Same for the bit buffer.
X
XThe remaining values in the RNGdata structure are derived
Xfrom the RNG's header file upon initialization.  For more information on
Xthe values of these fields, see Section 6.2.1.
X
X
XB RNGs Currently Installed in the Package
X-----------------------------------------
X-----------------------------------------
X
XThere are currently ten RNGs installed in the mrandom package.  This
Xappendix provides brief descriptions of each of them.  References are
Xprovided for those who are interested in finding out about the RNGs in
Xmore detail.
X
XRNG algorithm 0: A trivial RNG
X    A trivial RNG is included in the package, primarily for testing
X    purposes.  The generates it produces are not ``random'' in
X    virtually any sense of the word; it simply produces generates
X    from an arithmetical progression determined by its initial
X    seeds.  For example, if it is seeded with 5 and 7, respectively,
X    it will produce the sequence 5, 12, 19, 26, etc.
X
X    This RNG takes two longs as seeds.  It returns generates of type
X    long.
X
XRNG algorithm 1: 4.3bsd random
X    This is UNIX 4.3bsd random.  It is a 31-bit nonlinear
X    additive feedback generator with a range of 2^31 and a period of
X    approximately 16*2^31-1.  It is nominally able to save and
X    restore state, but its state-saving code is buggy.  Therefore,
X    when using random with the mrandom package, no more than one RNG
X    should use random at a time.
X
X    This RNG takes a single long as a seed.  It returns generates of
X    type long.
X
XRNG algorithm 2: the Knuth/Bentley prand
X    This lagged-Fibonacci RNG was introduced by Jon Bentley in his
X    ``Software Exploratorium'' column in Unix Review, Vol. 10, No.
X    6, June 1992, and is based on one first presented in Donald E.
X    Knuth's The Art of Computer Programming , Vol. 2,
X    Addison-Wesley, Reading, Mass., 1981.  It has a range of
X    1,000,000,000.
X
X    This RNG takes a single long as a seed.  It returns generates of
X    type long.
X
XRNG algorithm 3: The Portable Combined RNG
X    This combined prime multiplicative congruential RNG was
X    developed based on algorithms and selections of prime numbers
X    presented in ``Efficient and Portable Combined Random Number
X    Generators,'' Pierre L'Ecuyer, Communications of the ACM, Vol.
X    10, No.  6, June 1992, and ``Random Number Generators: Good Ones
X    are Hard to Find,'' Stephen Park and Keith Miller,
X    Communications of the ACM, Vol.  31, No. 10, October 1992.  It
X    has a range of 2147483561.
X
X    This RNG takes two longs as seeds.  It returns generates of type
X    long.
X
XRNG algorithm 4: 4.3bsd nrand48
X    This is UNIX 4.3bsd nrand48.  It produces generates using a
X    linear congruential algorithm and 48-bit integer arithmetic.  It
X    has a range of 2^31.
X
X    This RNG takes three unsigned shorts as seeds.  They are passed
X    to the seeding procedure as two longs, and are interpreted in
X    the following way:
X
X        * The 16 least significant bits of the second long is
X        the first seed.
X        * The 16 least significant bits of the first long is the
X        second seed.
X        * The 16 most significant bits of the first long is the
X        third seed.
X        * The 16 most significant bits of the second long is
X        ignored.
X
X    This RNG returns generates of type long.
X
XRNG algorithm 5: 4.3bsd rand
X    This is UNIX 4.3bsd rand.  It uses a multiplicative congruential
X    algorithm.  It has a period of 2^32 and a range of 2^31.
X
X    This RNG takes a single long as a seed.  It returns generates of
X    type long.
X
XRNG algorithm 6, 7, and 8: Press and Teukolsky's ran0, ran1, and ran2
X    These three multiplicative congruential RNGs are adapted from
X    those presented in ``Portable Random Number Generators,''
X    William H.  Press and Saul A. Teukolsky, Computers in Physics,
X    Vol. 6, No. 5, Sep/Oct 1992.  They all have a period of
X    2^31-2 and a range of 2^31-1.
X
X    These RNGs take a single long as a seed.  They return generates
X    of type double.
X
XRNG algorithm 9: Marsaglia's Ultra RNG
X
X    We obtained the source code for this generator by anonymous ftp
X    from nic.funit.fi (take the file fsultra.zip from the directory
X    /pub/msdos/science/math/fsultra).  A note in the readme file
X    says: ``To obtain permission to incorporate this program into
X    any commercial product, please contact the authors at the e-mail
X    address given above [[e-mail] or [e-mail]] or
X    at Department of Statistics and Supercomputer Computations
X    Research Institute, Florida State University, Tallahassee, FL
X    32306.''  This RNG is one of those originally presented in ``A
X    New Class of Random Number Generators,'' George Marsaglia and
X    Arif Zaman, The Annals of Applied Probability, Vol.  1, No. 3,
X    1991.  It is a ``subtract-with-borrow'' generator with a range
X    of 2^32 and a staggering period of 10^354.
X
X    This RNG takes two unsigned longs as seeds.  It returns
X    generates of type double.
END_OF_FILE
if test 50066 -ne `wc -c <'doc/mrandom.txt'`; then
    echo shar: \"'doc/mrandom.txt'\" unpacked with wrong size!
fi
# end of 'doc/mrandom.txt'
fi
if test -f 'src/mrandom.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/mrandom.c'\"
else
echo shar: Extracting \"'src/mrandom.c'\" \(42469 characters\)
sed "s/^X//" >'src/mrandom.c' <<'END_OF_FILE'
X/* mrandom.c 3.1 5/28/93 */
X/*
X *        mrandom.c
X *
X *      An interface for random number generators (RNGs)
X *
X *    Original Implementation:
X *        Clark Thomborson, September 1991
X *    Modifications:
X *       - give uniform interface to other rngs
X *       - remove bias from integer generator
X *       - avoid possible infinite loop in mrandom(m)
X *       - allow access to multiple simultaneous rngs
X *       - add interface to nrand48()
X *       - add interface to rand()
X *       - tuned for speed
X *        Clark Thomborson, May 1992
X *
X *      Version 3.0:
X *              Robert Plotkin, May 1993
X *      Modifications include:
X *         - Standardized interface to underlying RNGs
X *         - Added buffered generating routines
X *         - Added interfaces to ran0, ran1, ran2, Ultra RNGs
X *         - Added routines for generating bits
X *         - Added split feature
X *         - Allow choice of algorithms in mrandomrv
X *         - Dynamic allocation of RNGs
X *
X *    This material is based upon work supported by the National
X *    Science Foundation under grant number MIP-9023238.  The
X *    Government has certain rights in this material.
X *
X *    Any opinions, findings, and conclusions or recommendations
X *    expressed in this material are those of the author and do
X *    not necessarily reflect the view of the National Science
X *    Foundation.
X *
X *    This code is neither copyrighted nor patented, although we are
X *    not sure about the status of some of the routines it calls.
X */
X
X#ifndef MRANDOM
X#include "mrandom.h"
X#endif
X
X/*****************/
X/* Include files */
X/*****************/
X#include <malloc.h> /* To allocate space for RNGs*/
X#include <math.h>   /* we need floor() */
X#include <stdio.h>  /* We need FILE */
X#include <strings.h>
X#include <values.h> /* We need MAXLONG */
X
X/*******************************************************/
X/* Header files for RNGs currently included in package */
X/*******************************************************/
X#include "rngs.h"    /* "Built-in" rngs:
X            Trivial RNG:               Algorithm #0
X            4.3bsd random:             Algorithm #1
X            4.3bsd nrand48:            Algorithm #4
X            4.3bsd rand:               Algorithm #5 */
X#include "bentley.h" /* Knuth/Bentley prand:       Algorithm #2 */
X#include "pcrand.h"  /* The portable combined RNG: Algorithm #3 */
X#include "ran0.h"    /* Press & Teukolsky's ran0:  Algorithm #6 */
X#include "ran1.h"    /* Press & Teukolsky's ran1:  Algorithm #7 */
X#include "ran2.h"    /* Press & Teukolsky's ran2:  Algorithm #8 */
X#include "ultra.h"   /* Marsaglia's Ultra RNG:     Algorithm #9 */
X
X/*****************************************/
X/* Arrays to hold information about RNGs */
X/*****************************************/
Xchar RNGname_a[NUM_RNGS][RNGIDSTRLEN]={RNGname_0, RNGname_1,
X                     RNGname_2, RNGname_3,
X                     RNGname_4, RNGname_5,
X                     RNGname_6, RNGname_7,
X                     RNGname_8, RNGname_9};
Xlong statesize_a[NUM_RNGS]={RNGstatesize_0, RNGstatesize_1,
X                  RNGstatesize_2, RNGstatesize_3,
X                  RNGstatesize_4, RNGstatesize_5,
X                  RNGstatesize_6, RNGstatesize_7,
X                  RNGstatesize_8, RNGstatesize_9},
X  seedsize_a[NUM_RNGS]={RNGseedsize_0, RNGseedsize_1, RNGseedsize_2,
X              RNGseedsize_3, RNGseedsize_4, RNGseedsize_5,
X              RNGseedsize_6, RNGseedsize_7, RNGseedsize_8,
X              RNGseedsize_9};
Xdouble range_a[NUM_RNGS]={RNGrange_0, RNGrange_1, RNGrange_2,
X                RNGrange_3, RNGrange_4, RNGrange_5,
X                RNGrange_6, RNGrange_7, RNGrange_8,
X                RNGrange_9};
Xsigned int returns_a[NUM_RNGS]={RNGreturns_0, RNGreturns_1,
X                  RNGreturns_2, RNGreturns_3,
X                  RNGreturns_4, RNGreturns_5,
X                  RNGreturns_6, RNGreturns_7,
X                  RNGreturns_8, RNGreturns_9};
X  int statetype_a[NUM_RNGS]={RNGstatetype_0, RNGstatetype_1,
X                 RNGstatetype_2, RNGstatetype_3,
X                 RNGstatetype_4, RNGstatetype_5,
X                 RNGstatetype_6, RNGstatetype_7,
X                 RNGstatetype_8, RNGstatetype_9};
X
X
Xtypedef void (*SP)(RNGdata *, long *);
Xtypedef double (*DPP)(RNGdata *);
Xtypedef long (*LPP)(RNGdata *);
Xtypedef int (*CP)(RNGdata *);
XSP seed_a[NUM_RNGS] = {RNGseed_0, RNGseed_1, RNGseed_2, RNGseed_3,
X             RNGseed_4, RNGseed_5, RNGseed_6, RNGseed_7,
X             RNGseed_8, RNGseed_9};
XDPP dgen_a[NUM_RNGS] = {RNGdgen_0, RNGdgen_1, RNGdgen_2,
X               RNGdgen_3, RNGdgen_4, RNGdgen_5,
X               RNGdgen_6, RNGdgen_7,  RNGdgen_8,
X               RNGdgen_9};
XLPP lgen_a[NUM_RNGS] = {RNGlgen_0, RNGlgen_1, RNGlgen_2, RNGlgen_3,
X               RNGlgen_4, RNGlgen_5, RNGlgen_6,
X               RNGlgen_7, RNGlgen_8, RNGlgen_9};
XCP check_a[NUM_RNGS] = {RNGcheck_0, RNGcheck_1, RNGcheck_2,
X              RNGcheck_3, RNGcheck_4, RNGcheck_5,
X              RNGcheck_6, RNGcheck_7, RNGcheck_8, RNGcheck_9};
X
X/* The most-recently initialized or restarted RNG */
XRNGdata *mru_rng = 0;
X
X/* Format of the RNG statefile */
X#define RNGfileLINE0   "RNG statefile for algorithm %ld, "
X#define RNGfileLINE1   "Buffer size = %ld bytes\n"
X#define RNGfileLINE2   "Initial seed table =\n"
X#define RNGfileLINE3   \
X "Number of calls to underlying RNG after seeding = %ld billion + %ld\n"
X#define RNGfileLINE4   "Next value in this pseudorandom sequence = \
X%08lx\n"
X#define RNGfileLINE5   "This RNG returns every %ld generates\n"
X#define RNGfileLINE6   "This RNG uses mrandom algorithm %ld\n"
X#define RNGfileLINE7   "RNG state table =\n"
Xchar *RNGfileLINEn[]={"",
X            "   %02lx %02lx %02lx %02lx\n",
X            "   %04lx %04lx %04lx %04lx\n",
X            "",
X            "   %08lx %08lx %08lx %08lx\n"},
X  *RNGfileLINEnn[]={"", "   %02lx", "   %04lx", "", "   %08lx"};
X
X/* Error message for invalid RNGs */
X#define BADRNG_ERR \
X fprintf(stderr, "RNG must be initialized or restarted before use!\n"); \
X fflush(stderr);\
X exit(1)
X
X/**************************************************/
X/* Auxiliary macros for RNG generating procedures */
X/**************************************************/
X/*************************************/
X/* set_rng                           */
X/* Sets rng to point to a valid      */
X/*  RNGdata structure, if any        */
X/*  can be found.                    */
X/* Otherwise it calls                */
X/*  the error routine no_rng().      */
X/*************************************/
X#define set_rng(rng)   if (rng == 0) { \
X                if (mru_rng) { \
X                    rng = mru_rng; \
X                } else { \
X                    BADRNG_ERR; \
X                } \
X            }
X
X/*************************************/
X/* inc_counts                        */
X/* Increment an RNG counter n times  */
X/*************************************/
X/* Using two while loops is necessary.  If RNGcount1 is added to n
X *  before n is reduced to a number less than BILLION, it is
X *  possible that RNGcount1+n will overflow.
X * Possible overflow of count2 is not handled, because on today's
X *  computers it will not occur.
X */
X#define inc_counts(rng,n) \
X  while (n >= BILLION) { \
X    n -= BILLION; \
X    RNGcount2 += 1; \
X  } \
X  RNGcount1 += n; \
X  while (RNGcount1 >= BILLION) { \
X    RNGcount1 -= BILLION; \
X    RNGcount2 += 1; \
X  }
X
X/*************************************/
X/* fill_buffer                       */
X/* Fill an RNG buffer                */
X/*************************************/
Xvoid fill_buffer(rng)
X     RNGdata *rng;
X{
X  long i,j; /* Loop variables */
X
X  inc_counts(rng,RNGbuffer.size);
X
X  if (RNGreturns == RET_LONG) { /* Underlying RNG returns longs */
X    LPP proc;
X    long *temp_l; /* Temporary pointer for efficiency */
X    proc=lgen_a[RNGalg];
X    RNGbuffer.lbufp=temp_l=RNGbuffer.lbuf;
X
X    for (i=0; i<RNGbuffer.size; i++) { /* Fill the buffer */
X      *temp_l++  = proc(rng);
X      for (j=0;j<RNGsplit;j++) /* Handle split */
X    proc(rng);
X    }
X  }
X  else if (RNGreturns == RET_DOUBLE) { /* Underlying RNG returns */
X                                       /* doubles */
X    DPP proc;
X    double *temp_d; /* Temporary pointer for efficiency */
X    proc=dgen_a[RNGalg];
X    RNGbuffer.dbufp=temp_d=RNGbuffer.dbuf;
X
X    for (i=0; i<RNGbuffer.size; i++) { /* Fill the buffer */
X      *temp_d++ = proc(rng);
X      for (j=0;j<RNGsplit;j++) /* Handle split */
X    proc(rng);
X    }
X  }
X  RNGbuffer.nleft=RNGbuffer.size-1; /* Buffer is full */
X}
X
X/*************************************/
X/* bfill_vector                      */
X/* Fill a vector, with buffering,    */
X/*  with a function of               *
X/*  the underlying RNG               */
X/*************************************/
X#define bfill_vector(function) \
X  while (n--) { \
X         if (RNGbuffer.nleft-- == 0) \
X           fill_buffer(rng); \
X         *v++ = function; \
X         }
X
X/*************************************/
X/* random_setup                      */
X/* Handle special case of n == 0     */
X/*************************************/
X#define random_setup()\
X  set_rng(rng); \
X  if (n == 0) { \
X    n=1; \
X    v=vdef; \
X  } \
X  vreturn=v
X
X/*************************************/
X/* randomrv                          */
X/* Kernel of buffered routines       */
X/*************************************/
X/* The order of the cases is important for maximum efficiency.
X * The order is optimized for the case in which a request is made for
X *  the same type as the value returned by the underlying RNG
X */
X#define randomrv(case1, function1, case2, function2)\
X  random_setup(); \
X  if (RNGreturns == case1) \
X    bfill_vector(function1) \
X  else if (RNGreturns == case2) \
X    bfill_vector(function2) \
X  else \
X    while (n--) \
X      *v++ = 0; \
X  return(*vreturn)
X
X/*************************************/
X/* fill_vector                       */
X/* Fill a vector, without buffering, */
X/*  with a function of               */
X/*  the underlying RNG               */
X/*************************************/
X#define fill_vector(function) \
X  while (n--) { \
X         *v++ = function; \
X         for(j=0;j<RNGsplit;j++) \
X           proc(rng); \
X         }
X
X/*************************************/
X/* xrandom_setup                     */
X/* Setup common to all xrandoms      */
X/*************************************/
X#define xrandom_setup() \
X  random_setup(); \
X  inc_counts(rng,n)
X
X/* Auxiliary routines for 32-bit generators */
Xdouble longtodouble(l)
X     long l;
X{
X  if (l > 0)
X    return((double)l);
X  else if (l < 0)
X    return((double)(l & 0x7fffffff) + 2147483648.0); /* 2^31 */
X  else
X    return(4294967296.0); /* 2^32 */
X}
X
X/* Coded as a macro for efficiency */
X#define doubletolong(d) ((long)((unsigned long)(d)))
X
X/*******************************************************************/
X/* Procedures for generating pseudorandom numbers                  */
X/* These procedures write a vector v of length n, by repeating the */
X/* following two steps n times:                                    */
X/*      - place the next generate in v                             */
X/*      - discard the next k generates, where k is the split value */
X/*         of the RNG                                              */
X/*                                                                 */
X/* Special cases:                                                  */
X/*  Argument  Value  Meaning                                       */
X/*  --------  -----  -------                                       */
X/*  rng         0    Use most recently used or initialized rng.    */
X/*  n           0    Return  a single generate.                    */
X/*                   The value of v is ignored(it may be set to 0).*/
X/*                                                                 */
X/* All procedures return the first entry in the vector v.          */
X/*******************************************************************/
X
X/***********************************************/
X/* Buffered calling procedures.                */
X/***********************************************/
X/*******************/
X/* drandomrv       */
X/* Returns doubles */
X/*******************/
Xdouble drandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     double v[];
X{
X  long i;
X  double vdef[1],*vreturn;
X
X  if (RNGrange > MAXLONG + 1.0) {
X    randomrv(RET_DOUBLE, *RNGbuffer.dbufp++,
X         RET_LONG, longtodouble(*RNGbuffer.lbufp++)/RNGrange);
X  }
X  else {
X    randomrv(RET_DOUBLE, *RNGbuffer.dbufp++,
X         RET_LONG, *RNGbuffer.lbufp++/RNGrange);
X  }
X}
X
X/******************/
X/* frandomrv      */
X/* Returns floats */
X/******************/
Xfloat frandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     float v[];
X{
X  long i;
X  float vdef[1],*vreturn;
X
X  if (RNGrange > MAXLONG + 1.0) {
X    randomrv(RET_DOUBLE, *RNGbuffer.dbufp++,
X         RET_LONG, longtodouble(*RNGbuffer.lbufp++)/RNGrange);
X  }
X  else {
X    randomrv(RET_DOUBLE, *RNGbuffer.dbufp++,
X         RET_LONG, *RNGbuffer.lbufp++/RNGrange);
X  }
X}
X
X/*****************/
X/* lrandomrv     */
X/* Returns longs */
X/*****************/
Xlong lrandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     long v[];
X{
X  long vdef[1],*vreturn;
X
X  if (RNGrange > MAXLONG + 1.0) {
X    randomrv(RET_LONG, *RNGbuffer.lbufp++,
X         RET_DOUBLE, doubletolong(*RNGbuffer.dbufp++ * RNGrange));
X  }
X  else {
X    randomrv(RET_LONG, *RNGbuffer.lbufp++,
X         RET_DOUBLE, *RNGbuffer.dbufp++ * RNGrange);
X  }
X}
X
X/******************************************/
X/* brandomrv                              */
X/* Returns "high quality" bits,           */
X/*  i.e. each generate from the underlying*/
X/*  RNG is used to generate one bit.      */
X/******************************************/
Xint brandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     int v[];
X{
X  long i,j;
X  int vdef[1],*vreturn;
X
X  random_setup();
X
X  for (i=0;i<n;i++) {
X    if (RNGbuffer.nbleft-- == 0) { /* Bit buffer is empty, */
X      /* so fill it */
X      for (j=0;j<BBUF_SIZE;j++)
X    RNGbuffer.bbuf[j]=drandomr(rng)*2.0;
X      RNGbuffer.bbufp=RNGbuffer.bbuf;
X      RNGbuffer.nbleft=BBUF_SIZE-1; /* Bit buffer is full */
X    }
X    *v++ = *RNGbuffer.bbufp++;
X  }
X  return(*vreturn);
X}
X
X/******************************************/
X/* brandomrv_f                            */
X/* Returns "low quality" fast bits,       */
X/*  i.e. each generate from the underlying*/
X/*  RNG is used to generate 32 bits.      */
X/* Return -1 if the range of the RNG      */
X/*  is not 2^32.                          */
X/******************************************/
X/* Splitting of this fast bit stream is not currently supported. */
Xint brandomrv_f(rng,n,v)
X     RNGdata *rng;
X     long n;
X     int v[];
X{
X  long i,j,r;
X  int vdef[1],*vreturn;
X
X  random_setup();
X
X  for (i=0;i<n;i++) {
X    if (RNGbuffer.nbleft-- == 0) { /* Bit buffer is empty, */
X      /* so fill it */
X      r=lrandomr(rng);
X      for (j=0;j<BBUF_SIZE;j++) {
X    RNGbuffer.bbuf[j] = r & 1;
X    r = r >> 1;
X      }
X      RNGbuffer.bbufp=RNGbuffer.bbuf;
X      RNGbuffer.nbleft=BBUF_SIZE-1; /* Bit buffer is full */
X    }
X    *v++ = *RNGbuffer.bbufp++;
X  }
X  return(RNGrange == 4294967296.0 ? *vreturn: -1);
X}
X
X/***********************************************/
X/* Unbuffered calling procedures               */
X/***********************************************/
X/*******************/
X/* dxrandomrv      */
X/* Returns doubles */
X/*******************/
Xdouble dxrandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     double v[];
X{
X  long i,j;
X  double vdef[1],*vreturn;
X  double d;
X
X  xrandom_setup();
X  if (RNGreturns == RET_DOUBLE) {
X    DPP proc;
X    proc=dgen_a[RNGalg];
X    fill_vector(proc(rng));
X  }
X  else if (RNGreturns == RET_LONG) {
X    LPP proc;
X    proc=lgen_a[RNGalg];
X    if (RNGrange > MAXLONG + 1.0) /* 32-bit generator */
X      { fill_vector(longtodouble(lgen_a[RNGalg](rng))/RNGrange); }
X    else
X      { fill_vector((double)lgen_a[RNGalg](rng)/RNGrange); }
X  }
X  else
X    while (n--)
X      *v++ = 0;
X  return(*vreturn);
X}
X
X/******************/
X/* fxrandomrv     */
X/* Returns floats */
X/******************/
Xfloat fxrandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     float v[];
X{
X  long i,j;
X  float vdef[1],*vreturn;
X
X  xrandom_setup();
X  if (RNGreturns == RET_DOUBLE) {
X    DPP proc;
X    proc=dgen_a[RNGalg];
X    fill_vector(proc(rng));
X  }
X  else if (RNGreturns == RET_LONG) {
X    LPP proc;
X    proc=lgen_a[RNGalg];
X    if (RNGrange > MAXLONG + 1.0) /* 32-bit generator */
X      { fill_vector(longtodouble(lgen_a[RNGalg](rng))/RNGrange); }
X    else
X      { fill_vector((double)lgen_a[RNGalg](rng)/RNGrange); }
X  }
X  else
X    while (n--)
X      *v++ = 0;
X  return(*vreturn);
X}
X
X/*****************/
X/* lxrandomrv    */
X/* Returns longs */
X/*****************/
Xlong lxrandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     long v[];
X{
X  long i,j;
X  long vdef[1],*vreturn;
X
X  xrandom_setup();
X  if (RNGreturns == RET_LONG) {
X    LPP proc;
X    proc=lgen_a[RNGalg];
X    fill_vector(proc(rng));
X  }
X  else if (RNGreturns == RET_DOUBLE) {
X    DPP proc;
X    proc=dgen_a[RNGalg];
X    if (RNGrange > MAXLONG + 1.0) /* 32-bit generator */
X      { fill_vector(doubletolong(dgen_a[RNGalg](rng)*RNGrange)); }
X    else
X      { fill_vector(dgen_a[RNGalg](rng)*RNGrange); }
X  }
X  else
X    while (n--)
X      *v++ = 0;
X  return(*vreturn);
X}
X
X/*******************************/
X/* bxrandomrv                  */
X/* Return "high quality" bits. */
X/*******************************/
Xint bxrandomrv(rng,n,v)
X     RNGdata *rng;
X     long n;
X     int v[];
X{
X  long i,j;
X  int vdef[1],*vreturn;
X
X  xrandom_setup();
X  if (RNGreturns == RET_LONG) {
X    LPP proc;
X    proc=lgen_a[RNGalg];
X    fill_vector(((double)proc(rng)/RNGrange)*2.0);
X  }
X  else if (RNGreturns == RET_DOUBLE) {
X    DPP proc;
X    proc=dgen_a[RNGalg];
X    fill_vector(proc(rng)*2.0);
X    }
X  else
X    while (n--)
X      *v++ = 0;
X  return(*vreturn);
X}
X
X/*******************************/
X/* bxrandomrv_f                */
X/* Return "low quality" bits.  */
X/* Return -1 if range of RNG   */
X/*  is not 2^32                */
X/*******************************/
X/* This routine will "lose" random bits if they are not requested in
X *  multiples of (RNGsplit+1)*BBUF_SIZE, since bits are returned by
X *  pulling them off the random stream sequentially, using up exactly
X *  as many as are needed in order to generate the requested number
X *  of bits.
X * Splitting of this fast bit stream is not currently supported.
X */
Xint bxrandomrv_f(rng,n,v)
X     RNGdata *rng;
X     long n;
X     int v[];
X{
X  long i,j;
X  int vdef[1],*vreturn;
X  long value;
X  long index;
X
X  inc_counts(rng,n);
X  set_rng(rng);
X  random_setup();
X
X  index=0;
X
X  for (i=0;i<n;i++) {
X    if (index-- == 0) {
X      value=lxrandomr(rng);;
X      index=31;
X    }
X    v[i] = value & 1;
X    value = value >> 1;
X  }
X  return(RNGrange == 4294967296.0 ? *vreturn: -1);
X}
X
X/*************************************/
X/* seed_rng                          */
X/* Seeds the underlying RNG          */
X/*************************************/
Xvoid seed_rng(rng, seed)
X     RNGdata *rng;
X     long *seed;
X{
X  long i;
X
X  /* Keep a record of the seeds */
X  for (i=0;i<RNGseedsize;i++)
X    RNGseed[i]=seed[i];
X  RNGcount1 = 0;             /* Reset counter */
X  RNGcount2 = 0;
X  flush_rng(rng);            /* Flush the RNG's buffers */
X  seed_a[RNGalg](rng, seed); /* Seed the underlying RNG */
X  mru_rng = rng;             /* This RNG is now the most recently */
X                             /* initialized or restarted */
X}
X
X/*************************************/
X/* setRNGparams                      */
X/* Set RNGalg to alg, and set all    */
X/*  information which can be derived */
X/*  from RNGalg.                     */
X/* Return 0 and print message to     */
X/*  stderr if alg is out of range    */
X/* Return 1 otherwise                */
X/*************************************/
Xint setRNGparams(rng,alg)
X     RNGdata *rng;
X     long alg;
X{
X  if (alg<0 || alg>NUM_RNGS-1) /* alg is out of range */
X    return(0);
X  else { /* Set RNG parameters */
X    RNGalg = alg;
X    RNGstatesize = statesize_a[alg];
X    RNGseedsize = seedsize_a[alg];
X    RNGrange = range_a[alg];
X    strcpy(RNGname,RNGname_a[alg]);
X    RNGreturns = returns_a[alg];
X    RNGstatetype = statetype_a[alg];
X    return(1);
X  }
X}
X
X/*************************************/
X/* check_rng                         */
X/* Check the integrity of the RNG.   */
X/* Return 1 if ok, 0 if not.         */
X/*************************************/
Xint check_rng(rng)
X     RNGdata *rng;
X{
X  /* Temporary variables */
X  long statesize, seedsize;
X  double range;
X  signed int returns;
X  char name[RNGIDSTRLEN];
X  int statetype;
X
X  if (rng == 0)
X    return(0);
X
X  /* Have derived RNG values been overwritten? */
X  statesize = RNGstatesize;
X  seedsize = RNGseedsize;
X  range = RNGrange;
X  strcpy(name,RNGname);
X  returns = RNGreturns;
X  statetype = RNGstatetype;
X
X  if (!setRNGparams(rng,RNGalg))
X    return(0);
X
X  if (statesize != RNGstatesize || seedsize != RNGseedsize ||
X      range != RNGrange || strcmp(name,RNGname) || returns !=
X      RNGreturns || statetype != RNGstatetype)
X    return(0);
X
X  /* Now look at RNGstate (algorithm-specific tests) */
X  return(check_a[RNGalg](rng));
X}
X
X/*************************************/
X/* describe_rng                      */
X/* Write a short ASCII description   */
X/*  of the RNG to the user-supplied  */
X/*  string rngid, which must be of   */
X/*  length at least RNGIDSTRLEN.     */
X/* If the user has not initialized   */
X/*  the rng with init_rng() or       */
X/*  restart_rng(), abort with an     */
X/*  error message to stderr.         */
X/* Otherwise return rngid.           */
X/* Currently only supports RNGs with */
X/*  at most two seeds.               */
X/*************************************/
Xchar *describe_rng (rng,rngid)
X     RNGdata *rng;
X     char rngid[RNGIDSTRLEN];
X{
X  if (!check_rng(rng)) {
X    BADRNG_ERR;
X  }
X  sprintf(rngid, "RNG state identifier is (%ld, %ld: %ld, %ld; %ld,\
X%ld; %ld, %ld)\n",
X      RNGalg, RNGmrandom_alg, RNGseed[0], RNGseed[1],
X      RNGcount1, RNGcount2, RNGbuffer.size, RNGsplit);
X  return(rngid);
X}
X
X/*************************************/
X/* split_rng                         */
X/* Modify rng to "leapfrog" a number */
X/*  of generates determined by the   */
X/*  value of new_value               */
X/*  (new_value == 0 returns          */
X/*     every generate)               */
X/* Returns 0 if new_value < 0        */
X/* Returns 1 otherwise               */
X/* Exits on error if rng has not     */
X/*  been initialized                 */
X/*************************************/
Xint split_rng(rng,new_value)
X     RNGdata *rng;
X     long new_value;
X{
X  if (check_rng(rng)) {
X    if (new_value < 0)
X      return(0);
X    else {
X      RNGsplit=new_value;
X      return(1);
X    }
X  }
X  else {
X    BADRNG_ERR;
X  }
X}
X
X/*************************************/
X/* mralg_rng                         */
X/* Modify rng to use a different     */
X/*  algorithm for mrandom()          */
X/* Returns 0 if mralg is out of range*/
X/* Returns 1 otherwise               */
X/* Exits on error if rng has not     */
X/*  been initialized                 */
X/*************************************/
Xint mralg_rng(rng,new_value)
X     RNGdata *rng;
X     long new_value;
X{
X  if (check_rng(rng)) {
X    if (new_value == 0 || new_value == 1) {
X      RNGmrandom_alg=new_value;
X      return(1);
X    }
X    else
X      return(0);
X  }
X  else {
X    BADRNG_ERR;
X  }
X}
X
X/*************************************/
X/* range_rng                         */
X/* Return the range of the RNG       */
X/*  (i.e. rng is capable of producing*/
X/*  generates in the range           */
X/*  0...(range_rng(rng)-1)           */
X/* Exits on error if rng has not     */
X/*  been initialized                 */
X/*************************************/
Xdouble range_rng(rng)
X     RNGdata *rng;
X{
X  if (check_rng(rng))
X    return(RNGrange);
X  else {
X    BADRNG_ERR;
X  }
X}
X
X/* Auxiliary routine for allocating memory for an RNG */
XRNGdata *allocate_rng(rng, alg, bufsize)
X     RNGdata *rng;
X     long alg, bufsize;
X{
X  if (bufsize <= 0)
X    bufsize=1;
X  if ((rng = (RNGdata *)malloc(sizeof(RNGdata))) == 0)
X    return(0);
X  if (!setRNGparams(rng,alg))
X    return(0);
X  if ((RNGstate=(long *)calloc(RNGstatesize,sizeof(long))) == 0)
X    return(0);
X  /* Allocate and clear at least two elements for describe_rng */
X  if ((RNGseed=(long *)calloc(RNGseedsize > 2 ?
X                RNGseedsize:2,sizeof(long))) == 0)
X    return(0);
X  if (RNGreturns == RET_LONG) {
X    if ((RNGbuffer.lbuf=RNGbuffer.lbufp =
X     (long *) calloc(bufsize,sizeof(long))) == 0)
X      return(0);
X    RNGbuffer.dbuf=RNGbuffer.dbufp=0;
X  }
X  else if (RNGreturns == RET_DOUBLE) {
X    if ((RNGbuffer.dbuf=RNGbuffer.dbufp =
X     (double *) calloc(bufsize,sizeof(double))) == 0)
X      return(0);
X    RNGbuffer.lbuf=RNGbuffer.lbufp=0;
X  }
X  if ((RNGbuffer.bbuf=RNGbuffer.bbufp =
X       (int *) calloc(BBUF_SIZE,sizeof(int))) == 0)
X    return(0);
X  RNGbuffer.size=bufsize;
X  RNGbuffer.nleft=0;
X  RNGbuffer.nbleft=0;
X  return(rng);
X}
X
X/************/
X/* init_rng */
X/************/
X/* Initialize the general-purpose rng data area so that it holds the
X * state and other data required for the selected algorithm "alg", where
X *    alg = 0 is a trivial generator that returns state += seed2,
X *      where state is initially set to seed1.  Use for debugging only!
X *    alg = 1 is 4.3bsd random.c (non-linear additive feedback)
X *    alg = 2 is the Knuth/Bentley prand (lagged Fibonnacci)
X *    alg = 3 is the portable combined (multiplicative) RNG
X *    alg = 4 is 4.3bsd nrand48.c,
X *      alg = 5 is 4.3bsd rand
X *      alg = 6 is Press and Teukolsky's ran0
X *      alg = 7 is Press and Teukolsky's ran1
X *      alg = 8 is Press and Teukolsky's ran2
X *      alg = 9 is Marsaglia's Ultra RNG
X *
X * Note: Memory for rng is allocated by this routine.  Before calling
X *  this routine the user need only declare a pointer to the generator,
X *  for example
X *        RNGdata *myrng;
X *
X * The mrandom_alg parameter determines which algorithm will be used
X *  by mrandom when call with this RNG:
X *      0 = Thomborson's unbiased conversion
X *      1 = (int) (m * drandomr(rng))
X *
X * The bufsize parameter determines the number of entries in the
X *  buffer.  Buffer space is allocated dynamically during
X *  initialization.  Thirty-two entries are allocated for the bit
X *  buffer.
X *
X * The count1 parameter indicates how many times the selected RNG algorithm
X *  should be called prior to returning control to the calling routine.
X *  We recommend a high value, at least 10000, for this parameter, to minimize
X *  the ill effects of seeding.  The count2 parameter can be given a non-zero
X *  value if you wish to "cycle" the generator a huge number of times: it
X *  will be called (count2*1e9 + count1) times.  Probably count2 should always
X *  be 0 until computers become much, much faster than today.
X *
X * init_rng() returns a pointer to the initialized RNG unless an
X *  out-of-range argument is detected (rngalg >9 or < 0; count1 >1e9 or <0;
X *  or count2 <0), or if memory cannot be allocated for the RNG,
X *  in which case it returns a null pointer.
X *
X * Note: a single program may call init_rng() any number of times, to set up
X * multiple generators (possibly using more than one RNG algorithm), e.g with
X *        RNGdata *myrngs[8];
X *        long i,sum=0,seed[2];
X *              seed[1]=0;
X *        for (i=0; i<7; i++) {
X *          /* generator #i gets seed1 = i
X *                seed[0]=i;
X *          myrngs[i]=init_rng(2,0,seed,100000,0,1024);
X *        }
X *        /* our eighth RNG uses algorithm #3
X *              seed[7]=7;
X *        myrngs[7]=init_rng(3,0,seed,100000,0,1024);
X *        /* add 1-bit numbers, one from each generator
X *        for (i=0; i<7; i++) {
X *          sum += mrandom(&myrngs[i],2);
X *        }
X *
X * Warning: do not attempt to use multiple RNGdata areas for algorithm #1.
X *  The 4.3bsd random.c code has internal state that will not be modified
X *  correctly when you switch generators (this was a bug in the original
X *  implementation and it would be very difficult to fix here).
X *
X * Warning: Do NOT override previously-initialized RNGs with the results
X *  of this procedure.  If you have a pointer to a valid RNG and wish to
X *  initialize a new RNG using the same pointer, you should call
X *  kill_rng() before calling init_rng().  For example:
X *              ...
X *              i=lrandomr(rng); /* This RNG is in use
X *              kill_rng(rng);   /* Kill the RNG, and THEN
X *              rng=init_rng(3,1,seed,1000,0,256); /* init a new one
X *
X * We recommend that init_rng() be used very sparingly.  Except when
X * replicating experiments or debugging, it is better to restart an
X * existing generator (stored in a statefile) than to initialize a new one.
X */
XRNGdata *init_rng(alg, mrandom_alg, seed, count1, count2, bufsize)
X     long alg;         /* The RNG algorithm to use */
X     long mrandom_alg; /* Algorithm to use for mrandom:
X              0 = Thomborson's slow but unbasied conversion
X              1 = (int) (m*drandom()) */
X     long *seed;          /* Seed */
X     long count1, count2; /* Initial counts */
X     long bufsize;        /* Size of buffer */
X{
X  RNGdata *rng;
X  double x[64]; /* Temporary vector */
X  long i;       /* Loop variable    */
X
X  /* Are count1 and count2 in range? */
X  if (count1 >= BILLION || count1 < 0 || count2 < 0) {
X    return(0);
X  }
X  /* Is count2 nonzero? */
X  if (count2 != 0) {
X    fprintf(stderr, "Warning: this initialization will take a LONG time!\n");
X    fflush(stderr);
X  }
X  if ((rng=allocate_rng(rng,alg,bufsize))==0)
X    return(0); /* Allocate space for the RNG */
X  seed_rng(rng, seed);       /* Seed the RNG */
X  if (mralg_rng(rng, mrandom_alg) == 0)
X    return(0); /* Set default algorithm for mrandom */
X  if (split_rng(rng, SPLIT_DEF) == 0)
X    return(0); /* Set default split */
X  /* Cycle the generator, in blocks of 64 (for speed) */
X  for ( ; RNGcount2 < count2; ) { /* Get to the right BILLION count */
X    dxrandomrv(rng,64,x);
X  }
X  for ( ; RNGcount1+64 < count1; ) { /* Get close to the right count1 */
X    dxrandomrv(rng,64,x);
X  }
X  for ( ; RNGcount1 != count1; ) { /* and, finally, do the last few */
X    dxrandomrv(rng,1,x);
X  }
X  return(rng); /* normal exit */
X
X}
X
X/*******************************/
X/* kill_rng                    */
X/* Frees memory used by an RNG */
X/* Returns 0 if kill failed    */
X/* Returns 1 otherwise         */
X/*******************************/
Xint kill_rng(rng)
X     RNGdata *rng;
X{
X  if (check_rng(rng)) {
X    free(RNGstate);
X    free(RNGseed);
X    free(RNGbuffer.lbuf);
X    free(RNGbuffer.dbuf);
X    free(RNGbuffer.bbuf);
X    free(rng);
X    return(1);
X  }
X  else
X    return(0);
X}
X
X/***************************************/
X/* nextval                             */
X/* Return the next lxrandomrv() output */
X/*  without disturbing the RNG state.  */
X/***************************************/
Xlong nextval(rng)
X     RNGdata *rng;
X{
X  static long *state=0;
X  long i, r, tcount1, tcount2, retval;
X
X  if (state != 0)
X    free(state);
X  state = (long *)calloc(RNGstatesize,sizeof(long));
X  /* Check that this only allocates the first time */
X
X  /* Preserve old RNG state */
X  tcount1 = RNGcount1;
X  tcount2 = RNGcount2;
X  for (i=0; i<RNGstatesize; i++) {
X    state[i] = RNGstate[i];
X  }
X
X  /* Find the next value in this pseudorandom sequence */
X  retval=lxrandomr(rng);
X
X  /* Restore the RNG state */
X  RNGcount1 = tcount1;
X  RNGcount2 = tcount2;
X  /* Special fixup for 4.3bsd random()'s hidden state variables */
X  if (RNGalg == 1) for (i=1; i<31; i++) {
X    r = random();
X  }
X  for (i=0; i<RNGstatesize; i++) {
X    RNGstate[i] = state[i];
X  }
X
X  return(retval);
X}
X
X/* printvector */
X/* auxiliary routine for save_rng */
Xvoid printvector(fp,size,vector,type)
X     FILE *fp;
X     long size;
X     long *vector;
X     int type;
X{
X  long i,j,printval;
X  long windex,bindex;
X  long mask;
X
X  mask=0;
X  for(j=0;j < (type << 3); j++)
X    mask = (mask << 1) | 1;
X
X  windex=0;
X  bindex=31;
X  i=0;
X
X  while (windex<size) {
X    printval = 0;
X    for (j=0; j < (type << 3); j++) {
X      printval |= (vector[windex] & (1 << bindex));
X      if (bindex-- == 0) {
X    bindex = 31;
X    windex++;
X      }
X    }
X    fprintf(fp, RNGfileLINEnn[type],
X        (printval >> ((bindex+1)%32)) & mask);
X    if (++i == 4) {
X      fprintf(fp,"\n");
X      i = 0;
X    }
X  }
X  fprintf(fp,"\n");
X}
X
X/* scanvector */
X/* auxiliary routine for restart_rng */
Xvoid scanvector(fp,size,vector,type)
X     FILE *fp;
X     long size;
X     long *vector;
X     int type;
X{
X  long i,j,scanval;
X  long windex,bindex;
X
X  windex=-1;
X  bindex=0;
X  i=0;
X
X  while (windex<size) {
X    fscanf(fp,RNGfileLINEnn[type],&scanval);
X    for (j=0; j < (type << 3); j++) {
X      if (bindex-- == 0) {
X    bindex=31;
X    windex++;
X    vector[windex]=0;
X      }
X    }
X    vector[windex] |= (scanval << bindex);
X    if (++i == 4) {
X      fscanf(fp,"\n");
X      i = 0;
X    }
X  }
X  fscanf(fp,"\n");
X}
X
X/*********************************************/
X/* save_rng                                  */
X/* Save the RNG state to a statefile.        */
X/* Return 0 if RNG couldn't be saved.        */
X/* Returns 1 otherwise.                      */
X/*********************************************/
Xint save_rng(rng,filename)
X     RNGdata *rng;
X     char *filename;
X{
X  FILE *fp;
X  long i;
X  unsigned short *sarray;
X  RNGdata *temp_rng;
X
X  if (!check_rng(rng))
X    return(0);
X
X  /* Write the statefile */
X  fp = fopen(filename, "w");
X  if (!fp) /* Couldn't open file for writing */
X    return(0);
X
X  fprintf(fp, RNGfileLINE0, RNGalg);
X  fprintf(fp,RNGname_a[RNGalg]);
X  fprintf(fp, RNGfileLINE1,RNGbuffer.size);
X  fprintf(fp, RNGfileLINE2);
X  printvector(fp,RNGseedsize,RNGseed,RNGstatetype);
X  fprintf(fp, RNGfileLINE3, RNGcount2, RNGcount1);
X  fprintf(fp, RNGfileLINE4, nextval(rng));
X  fprintf(fp, RNGfileLINE5, RNGsplit+1);
X  fprintf(fp, RNGfileLINE6, RNGmrandom_alg);
X  fprintf(fp, RNGfileLINE7);
X  printvector(fp,RNGstatesize,RNGstate,RNGstatetype);
X  fclose(fp);
X
X  temp_rng=restart_rng(filename); /* Verify save */
X  if (temp_rng == 0)
X    return(0);
X  else {
X    kill_rng(temp_rng);
X    return(1);
X  }
X}
X
X/****************************************************************/
X/* restart_rng                                                  */
X/* Restart a generator from a statefile.                        */
X/* Return a null pointer if the restart failed due to a garbled */
X/*  or nonexistent statefile.                                   */
X/* Otherwise return a pointer to the restarted RNG.             */
X/* WARNING: An RNG which has been previously initialized using  */
X/*  init_rng() should NOT be overwritten with the return value  */
X/*  of this procedure.  In order to provide a "fresh" RNG for   */
X/*  this procedure, do one of the following:                    */
X/*     - Declare a new RNG                                      */
X/*     - Kill a previously initialized RNG using kill_rng()     */
X/****************************************************************/
XRNGdata *restart_rng(filename)
X     char *filename;
X{
X  FILE *fp;
X  unsigned short *sarray; /* for reading alg #4's statefile */
X  long i, m, r; /* temps */
X  long alg, bufsize;
X  int errflag; /* initially 0, becomes 1 if a problem is discovered */
X  RNGdata *rng;
X
X  /* If mru_rng == 0, we assume random() hasn't been called yet,
X   * so its internal state is at its startup values.
X   * If mru_rng != 0 and alg=1, then we've used random() already, so we
X   * must reset its internal state to its startup values.
X   */
X  if (mru_rng != 0 && mru_rng->rngalg == 1) {
X    m = (mru_rng->rngcount1%31) + mru_rng->rngcount2*(BILLION%31);
X    /* note: the hidden state variables are counters mod 31 */
X    for ( ; (m%31) != 0; m++) {
X      r = random();
X    }
X  }
X
X  /* Restore counter values, retrieve original RNG seeds and current state */
X  fp = fopen(filename, "r");
X  if (!fp)
X    return(0);
X  /* Read statefile */
X  fscanf(fp, RNGfileLINE0, &alg);
X  fscanf(fp,RNGname_a[alg]);
X  fscanf(fp, RNGfileLINE1, &bufsize);
X  if ((rng=allocate_rng(rng, alg, bufsize))==0)
X    return(0);
X  fscanf(fp, RNGfileLINE2);
X  scanvector(fp,RNGseedsize,RNGseed,RNGstatetype);
X  fscanf(fp, RNGfileLINE3, &RNGcount2, &RNGcount1);
X  fscanf(fp, RNGfileLINE4, &RNGnextval);
X  fscanf(fp, RNGfileLINE5, &RNGsplit);
X  RNGsplit--;
X  fscanf(fp, RNGfileLINE6, &RNGmrandom_alg);
X  fscanf(fp, RNGfileLINE7);
X  scanvector(fp,RNGstatesize,RNGstate,RNGstatetype);
X  fclose(fp);
X
X  errflag = 0; /* no errors detected yet */
X
X  /* If reconstruction will be rapid, do it as an error check. */
X  if (RNGcount1 < 1000 && RNGcount2 == 0) {
X    RNGdata *test_rng;
X    test_rng = init_rng(RNGalg, RNGmrandom_alg, RNGseed,
X            RNGcount1, RNGcount2, RNGbuffer.size);
X    /* See if we got to the same state */
X    for (i=0; i<RNGstatesize; i++) {
X      if ((test_rng->rngstate)[i] != RNGstate[i]) {
X        errflag = 1;
X      }
X    }
X    kill_rng(test_rng);
X  } else { /* Just update mru_rng */
X    /* first, some special hacks for alg 1 */
X    if (RNGalg == 1) {
X      /* tell random() we've got a 31-word statefile */
X      RNGstate[0] = MRNGSTATE0;
X      setstate(RNGstate);
X      /* and modify random()'s hidden state to correspond to the RNGcount */
X      m = RNGcount1%31 + RNGcount2*(BILLION%31);
X      for (i=0 ; i<(m%31); i++) {
X        r = random();
X      }
X    }
X    mru_rng = rng; /* remember this rng, for use by mrandom() and frandom() */
X  }
X
X  /* see if RNG state looks ok */
X  errflag = errflag || !check_rng(rng);
X
X  /* Check nextval() operation */
X  if (nextval(rng) != nextval(rng)) {
X    errflag = 1;
X  }
X
X  if (errflag) {
X    kill_rng(rng);
X    return(0);
X  }
X  return(rng);
X}
X
X/*********************************************/
X/* flush_rng                                 */
X/* Flush the contents of the RNG's buffers   */
X/* Returns 1 upon success, 0 upon failure    */
X/*********************************************/
Xint flush_rng(rng)
X     RNGdata *rng;
X{
X  if (check_rng(rng)) {
X    if (RNGreturns == RET_LONG) {
X      RNGbuffer.lbufp=RNGbuffer.lbuf;
X      RNGbuffer.dbuf=RNGbuffer.dbufp=0;
X    }
X    else if (RNGreturns == RET_DOUBLE) {
X      RNGbuffer.dbufp=RNGbuffer.dbuf;
X      RNGbuffer.lbuf=RNGbuffer.lbufp=0;
X    }
X    RNGbuffer.nleft=0;
X    RNGbuffer.bbufp=RNGbuffer.bbuf;
X    RNGbuffer.nbleft=0;
X    return(1);
X  }
X  else
X    return(0);
X}
X
X/*************/
X/* mrandomrv */
X/*************/
X/* Generate a length-n vector v of random longs, uniformly distributed
X * in the range 0..m-1, using the indicated rng.  Return a copy of the
X * first random variate, v[0].
X *
X * Special-case parameter values: if rng==0, use the RNG that was
X * the most-recently initialized or restarted; if n==0, return one
X * random variate and don't write into v[].
X *
X * Our code does not have a deterministic bias for any m, unlike the
X * typical "good" code
X *        (int) floor( drandom() * (double) m )
X * or the commonly-used, but hazardous (because it exposes the flaws
X * in many RNGs) code
X *        random()%m
X * We remove the bias by making multiple calls (on rare occasions)
X * to the underlying generator.  The expected number of RNG calls
X * is upper-bounded by n/(1 - (RNGrange%m)/RNGrange) < 2n.
X *
X * The program is aborted, with an error message, if
X *  m exceeds the range of the RNG.
X *
X * The program will also abort, again with an error message to stderr,
X * if the generator is behaving so non-randomly that our multiple-call
X * bias-removal algorithm makes an excessive number of calls to the
X * underlying generator.
X */
Xlong mrandomrv(rng,m,n,v)
X     RNGdata *rng;
X     long m,n,v[];
X{
X  long i,yield; /* temps */
X  long ncalls; /* counts number of RNG calls during this mrandomrv call */
X  long vdef[1],*vreturn; /* default value for v */
X  char rngdesc[RNGIDSTRLEN]; /* for a describe_rng() call */
X  double d; /* temporary variable */
X
X  /* we can avoid some calcs&tests if m and rng range are unchanged */
X  static long lastm=0,lastrangem1=0; /* m on last call */
X  static double lastrange=0;         /* range on last call */
X  static double maxv;  /* largest "useful" RNG output for this m */
X  static double ifreq; /* multiplicity of RNG -> mrandomr() mapping */
X  static double dm;    /* floated value of m */
X  static int is32bit;  /* is rng is 32-bit generator? */
X
X  set_rng(rng);
X
X  random_setup();
X
X  if (m != lastm) {
X    dm=longtodouble(m);
X    lastm=m;
X    if (dm > RNGrange) { /* Is m in range? */
X      fprintf(stderr,
X          "Error: mrandom() argument value, %ld, is out of range.\n", m);
X      fprintf(stderr,
X          "The range of this RNG is %.0f.\n", RNGrange);
X      fflush(stderr);
X      exit(1);
X    }
X  }
X  switch (RNGmrandom_alg) {
X  case 0: /* Thomborson's unbiased conversion */
X    /* see if RNGrange has changed recently */
X    if (RNGrange != lastrange) {
X      /* yes, save current RNGrange ... */
X      lastrange = RNGrange;
X
X      /* ... and recalculate ifreq and maxv */
X      /* Compute ifreq = multiplicity of RNG -> mrandom(m) mapping,
X       * and maxv = the largest "useful" RNG output for this m.
X       *
X       * You might want to rewrite this code if you don't have FP hardware.
X       */
X      ifreq = floor(RNGrange/dm);
X      maxv = RNGrange - (int)(RNGrange - ifreq*dm);
X      is32bit = (RNGrange > MAXLONG + 1.0);
X    }
X
X    /* the expected # of calls to underlying RNG is */
X    /* n/(1-xdiscard/range) < 2n */
X
X    ncalls = 0; /* number of RNG calls so far */
X    for (yield=0; yield<n;  ) {
X
X      /* fill (or re-fill) v[] */
X      lxrandomrv(rng,n-yield,&v[yield]);
X      /* make sure ncalls doesn't get ridiculously large */
X      ncalls += n-yield;
X      if (ncalls > 3*n+300) {
X    /* For yield ==n, mean(ncalls) < 2*n; std dev < sqrt(2*n);
X     * If ncalls > 3*n + 300, we are dozens of stddevs away from mean!
X     */
X    fprintf(stderr, "Trouble in mrandomrv, m = %ld\n",m);
X    fprintf(stderr, "Aborting program: this RNG is apparently stuck!\n");
X    fprintf(stderr, describe_rng(rng,rngdesc));
X    exit(1);
X      }
X
X      /* the following steps are coded for 32-bit and non-32-bit RNGs */
X      /* for reasons of efficiency (32-bit RNGs will be slower)       */
X      /* first, for all generators except 32-bit generators */
X      if (!is32bit) {
X    for (  ; yield<n; yield++) {
X      if (v[yield] > maxv) {
X        break;
X      }
X    }
X    for (i=yield+1; i<n; i++) {
X      if (v[i] <= maxv) {
X        v[yield] = v[i];
X        yield++;
X      }
X    }
X    for (i=0; i<n; i++) {
X      v[i] = (long) ((double)v[i] / ifreq);
X    }
X      }
X      else { /* for 32-bit generators */
X    /* find first out-of-range v[], if any, */
X    for (  ; yield<n; yield++) {
X      if (longtodouble(v[yield]) > maxv) {
X        break;
X      }
X    }
X    /* and move in-range values to front of v[] */
X    for (i=yield+1; i<n; i++) {
X      if (longtodouble(v[i]) <= maxv) {
X        v[yield] = v[i];
X        yield++;
X      }
X    }
X    /* map v[] values, with ifreq:1 multiplicity, */
X    /* into the range 0..m-1 */
X    for (i=0; i<n; i++) {
X      v[i] = (long) (longtodouble(v[i]) / ifreq);
X    }
X      }
X    }
X    return (*vreturn);
X    break;
X  case 1: /* (long) m*dxrandomr() */
X    for (i=0;i<n;i++)
X      v[i]= (long) (dm*dxrandomr(rng));
X    return(*vreturn);
X    break;
X  default:
X    break;
X  }
X}
END_OF_FILE
if test 42469 -ne `wc -c <'src/mrandom.c'`; then
    echo shar: \"'src/mrandom.c'\" unpacked with wrong size!
fi
# end of 'src/mrandom.c'
fi
echo shar: End of archive 3 \(of 6\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 6 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0