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

Text File | 1994-05-06 | 87.7 KB | 2,533 lines

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

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

#! /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 2 (of 6)."
# Contents:  man/mrandom.3 man/mrandom.3.txt src/mrandom.h src/mrtest.c
# Wrapped by [e-mail] on Fri May  6 21:42:55 1994
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'man/mrandom.3' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'man/mrandom.3'\"
else
echo shar: Extracting \"'man/mrandom.3'\" \(21315 characters\)
sed "s/^X//" >'man/mrandom.3' <<'END_OF_FILE'
X.\" mrandom.3 3.1 5/28/93
X.TH MRANDOM 3 "5/28/93"
X.SH NAME
Xinit_rng, kill_rng, \
Xdxrandomrv, dxrandomr, dxrandomv, dxrandom, \
Xfxrandomrv, fxrandomr, fxrandomv, fxrandom, \
Xlxrandomrv, lxrandomr, lxrandomv, lxrandom, \
Xbxrandomrv, bxrandomr, bxrandomv, bxrandom, \
Xbxrandomrv_f, bxrandomr_f, bxrandomv_f, bxrandom_f, \
Xdrandomrv, drandomr, drandomv, drandom, \
Xfrandomrv, frandomr, frandomv, frandom, \
Xlrandomrv, lrandomr, lrandomv, lrandom, \
Xbrandomrv, brandomr, brandomv, brandom, \
Xbrandomrv_f, brandomr_f, brandomv_f, brandom_f, \
Xmrandomrv, mrandomr, mrandomv, mrandom, \
Xflush_rng, save_rng, restart_rng, seed_rng, check_rng, \
Xdescribe_rng, mralg_rng, split_rng, range_rng \
X\- a uniform interface to several random number generators
X.SH SYNOPSIS
X.nf
X.B RNGdata *init_rng(alg, mrandom_alg, seed, count1, count2, bufsize)
X.B long alg, mrandom_alg, *seed, count1, count2, bufsize;
X.LP
X.B int kill_rng(rng)
X.B RNGdata *rng;
X.LP
X.B double dxrandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B double v[];
X.LP
X.B float fxrandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B float v[];
X.LP
X.B long lxrandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B long v[];
X.LP
X.B long lxrandomv(n, v)
X.B long n;
X.B long v[];
X.LP
X.B int bxrandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B int v[];
X.LP
X.B int bxrandomrv_f(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B int v[];
X.LP
X.B double drandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B double v[];
X.LP
X.B float frandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B float v[];
X.LP
X.B long lrandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B long v[];
X.LP
X.B int brandomrv(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B int v[];
X.LP
X.B int brandomrv_f(rng, n, v)
X.B RNGdata *rng;
X.B long n;
X.B int v[];
X.LP
X.B long mrandomrv(rng, m, n, v)
X.B RNGdata *rng;
X.B long m;
X.B long n;
X.B long v[];
X.LP
X.B int save_rng(rng, filename)
X.B RNGdata *rng;
X.B char *filename;
X.LP
X.B RNGdata *restart_rng(rng, filename)
X.B char *filename;
X.LP
X.B void seed_rng(rng, seed)
X.B RNGdata *rng;
X.B long *seed;
X.LP
X.B int check_rng(rng);
X.B RNGdata *rng;
X.LP
X.B char *describe_rng(rng, rngid)
X.B RNGdata *rng;
X.B char rngid[RNGIDSTRLEN];
X.LP
X.B int mralg_rng(rng,new_value)
X.B RNGdata *rng;
X.B long new_value;
X.LP
X.B int split_rng(rng, new_value)
X.B RNGdata *rng;
X.B long new_value;
X.LP
X.B double range_rng(rng)
X.B RNGdata *rng;
X.fi
X.IX  "init_rng()"  ""  "\fLinit_rng\fP \(em initialize random generator"
X.IX  "kill_rng()"  ""  "\fLkill_rng\fP \(em kill random generator"
X.IX  "dxrandomrv()"  ""  "\fldxrandomrv\fP \(em generate vector of \
X    random doubles"
X.IX  "fxrandomrv()"  ""  "\flfxrandomrv\fP \(em generate vector of \
X    random floats"
X.IX  "lxrandomrv()"  ""  "\fllxrandomrv\fP \(em generate vector of \
X    random longs"
X.IX  "bxrandomrv()"  ""  "\flbxrandomrv\fP \(em generate vector of \
X    random bits, one bit from each generate"
X.IX  "bxrandomrv_f()"  ""  "\flbxrandomrv_f\fP \(em generate vector of \
X    random bits, 32 bits from each generate"
X.IX  "dxrandomrv()"  ""  "\fldxrandomrv\fP \(em generate buffered vector of \
X    random doubles"
X.IX  "fxrandomrv()"  ""  "\flfxrandomrv\fP \(em generate buffered vector of \
X    random floats"
X.IX  "lxrandomrv()"  ""  "\fllxrandomrv\fP \(em generate buffered vector of \
X    random longs"
X.IX  "bxrandomrv()"  ""  "\flbxrandomrv\fP \(em generate buffered vector of \
X    random bits, one bit from each generate"
X.IX  "bxrandomrv_f()"  ""  "\flbxrandomrv_f\fP \(em generate buffered \
X    vector of random bits, 32 bits from each generate"
X.IX  "mrandomrv()"  ""  "\fLmrandomrv\fP \(em generate vector of \
X    random integers mod m"
X.IX  "save_rng()"  "" "\fLsave_rng\fP \(em save random generator to file"
X.IX  "restart_rng()"  "" \
X    "\fLrestart_rng\fP \(em restart random generator from file"
X.IX  "seed_rng()"  ""  "\fLseed_rng\fP \(em seed random generator"
X.IX  "check_rng()"  ""  "\fLcheck_rng\fP \(em check integrity of \
X    random generator"
X.IX  "describe_rng()"  "" \
X    "\fLdescribe_rng\fP \(em construct short string describing \
X     generator state"
X.IX  "mralg_rng()"  ""  "\fLmralg_rng\fP \(em set mrandomrv algorithm \
X    number"
X.IX  "split_rng()"  ""  "\fLsplit_rng\fP \(em set generator split value"
X.IX  "range_rng()"  ""  "\fLrange_rng\fP \(em return generator range"
X.IX  "random number generator"  "\fLinit_rng()\fP"
X.IX  "random number generator"  "\fLkill_rng()\fP"
X.IX  "random number generator"  "\fLdxrandomrv()\fP"
X.IX  "random number generator"  "\fLfxrandomrv()\fP"
X.IX  "random number generator"  "\fLlxrandomrv()\fP"
X.IX  "random number generator"  "\fLbxrandomrv()\fP"
X.IX  "random number generator"  "\fLbxrandomrv_f()\fP"
X.IX  "random number generator"  "\fLdrandomrv()\fP"
X.IX  "random number generator"  "\fLfrandomrv()\fP"
X.IX  "random number generator"  "\fLlrandomrv()\fP"
X.IX  "random number generator"  "\fLbrandomrv()\fP"
X.IX  "random number generator"  "\fLbrandomrv_f()\fP"
X.IX  "random number generator"  "\fLmrandomrv()\fP"
X.IX  "random number generator"  "\fLdxrandomrv()\fP"
X.IX  "random number generator"  "\fLfxrandomrv()\fP"
X.IX  "random number generator"  "\fLlxrandomrv()\fP"
X.IX  "random number generator"  "\fLbxrandomrv()\fP"
X.IX  "random number generator"  "\fLflush_rng()\fP"
X.IX  "random number generator"  "\fLsave_rng()\fP"
X.IX  "random number generator"  "\fLrestart_rng()\fP"
X.IX  "random number generator"  "\fLseed_rng()\fP"
X.IX  "random number generator"  "\fLcheck_rng()\fP"
X.IX  "random number generator"  "\fLdescribe_rng()\fP"
X.IX  "random number generator"  "\fLmralg_rng()\fP"
X.IX  "random number generator"  "\fLsplit_rng()\fP"
X.IX  "random number generator"  "\fLrange_rng()\fP"
X.IX  "generate random numbers"  "\fLinit_rng()\fP"
X.IX  "generate random numbers"  "\fLkill_rng()\fP"
X.IX  "generate random numbers"  "\fLdxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLfxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLlxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLbxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLbxrandomrv_f()\fP"
X.IX  "generate random numbers"  "\fLdrandomrv()\fP"
X.IX  "generate random numbers"  "\fLfrandomrv()\fP"
X.IX  "generate random numbers"  "\fLlrandomrv()\fP"
X.IX  "generate random numbers"  "\fLbrandomrv()\fP"
X.IX  "generate random numbers"  "\fLbrandomrv_f()\fP"
X.IX  "generate random numbers"  "\fLmrandomrv()\fP"
X.IX  "generate random numbers"  "\fLdxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLfxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLlxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLbxrandomrv()\fP"
X.IX  "generate random numbers"  "\fLflush_rng()\fP"
X.IX  "generate random numbers"  "\fLsave_rng()\fP"
X.IX  "generate random numbers"  "\fLrestart_rng()\fP"
X.IX  "generate random numbers"  "\fLseed_rng()\fP"
X.IX  "generate random numbers"  "\fLcheck_rng()\fP"
X.IX  "generate random numbers"  "\fLdescribe_rng()\fP"
X.IX  "generate random numbers"  "\fLmralg_rng()\fP"
X.IX  "generate random numbers"  "\fLsplit_rng()\fP"
X.IX  "generate random numbers"  "\fLrange_rng()\fP"
X.SH DESCRIPTION
X.LP
X.I mrandomrv(rng,m,n,v)
Xgenerates n random integers in the range 0 to m-1, using the generator
Xrng, and storing them in vector v.
XAn initialization routine,
X.I init_rng(),
Xdescribed below, allows you to select an appropriate generator.
X
XUsing this package, instead of direct calls to random() or some
Xother RNG code, offers the following advantages.
X.IP 1.
XYou can switch to another RNG, and re-run a
Xportion of your experiment to check its validity,
Xwithout changing any of your code other
Xthan a single parameter in your
X.I
Xinit_rng()
Xcall.
X.IP 2.
XYou can use my unbiased method for generating random integers
Xin the range 0..m-1.  By contrast, the typical integer-generating codes
X"random()%m" or "(int) m * ((double) random()/MAXLONG)"
Xhave a slight bias for large m.
X.IP 3.
XFloating point numbers, uniformly distributed in [0.0, 1.0),
Xare available (as on most RNG interfaces).
X.IP 4.
XAs with most RNG interfaces, you can have many simultaneously-active RNGs.
X(Unfortunately, I was unable to find an efficient work-around for 4.3bsd
Xrandom()'s state-saving bug; you'll have to use one of the other generators
Xif you want to use multiple RNG streams in a single program invocation.)
X.IP 5.
XTime-efficient vectorized calls, returning multiple uniform variates,
Xare available.
X.IP 6
XThe ability to "split" RNGs to produce parallet output streams using the
X"leapfrog" method.
X.IP 7.
XThe package offers a shorthand notation for completely specifying the
Xalgorithm and current state of your RNG(s),
Xin an 80-character human-readable ASCII string.
XThis is very useful in experimental documentation and replication.
X.IP 8.
XA complete RNG state can be (serially) reconstructed from its
Xshorthand notation, and
X.IP 9.
XA file-I/O interface allows fast saves and restarts of complete
XRNG state vectors, without the time overhead of serial reconstruction.
X
X.LP
XAlso included in this software release is a (unsupported) driver routine
X.I mrtest.c.
XThis routine implements some of the simpler tests of
Xrandomness, e.g. equidistribution, pairwise (both short-
Xand long-range) correlation, and 3-tuple correlation [Marsaglia, 1985].
XThese tests were chosen to illustrate the use of the components
Xof the
X.I mrandom()
Xpackage, as well as to exhibit the known
Xdefects of the 4.3bsd generators rand() and nrand48().
XSee the man page on
X.I mrtest
Xfor more details.
X
XAnother use for
X.I mrtest.c
Xis as a template for your own calls to the mrandom() package.
X(I, for one, like to program by example...)
X
XBelow are detailed descriptions of the interfaces to the
Xvarious C-language functions in the mrandom() package.
XFor technical details, consult the file
X.I mrandom.tex,
Xalso included in this distribution.
X
XIn order to use an RNG, it must first be initialized.  This
Xis accomplished by first declaring a pointer to an RNGdata
Xstructure, and then calling
X.I init_rng()
X, which allocates memory for
Xthe RNG and readies it for use by the other routines in the
Xpackage.
X
X.I init_rng()
Xreturns a pointer to an initialized RNG.  A pointer
Xreturned by
X.I init_rng() is valid for use by all other mrandom routines.
X
XIn an init_rng() call, the
X.I alg
Xparameter should have a value between 0 and 9, to indicate which RNG
Xalgorithm is desired:
X.IP 1.
X4.3bsd random(), a non-linear additive feedback RNG,
X.IP 2.
XThe Knuth/Bentley prand(), a lagged-Fibonnacci generator
Xwith a state table of size 55 [Bentley, 1992],
X.IP 3.
X4.3bsd nrand48(),
Xa 48-bit multiplicative congruential RNG,
X.IP 4.
XL'Ecuyer's ``portable combined'' 32-bit multiplicative congruential
Xgenerator [L'Ecuyer, 1988],
X.IP 5.
X4.3bsd rand(),
Xthe discredited 32-bit multiplicative congruential RNG, and
X.IP 6-8.
XPress & Teukolsky's ran0, ran1, and ran2, respectively [Press &
XTeukolsky, 1992],
X.IP 9.
XMarsaglia's subtract-with-borrow Ultra generator [Marsaglia and Zaman, 1991],
X.IP 0.
Xfor testing purposes,
Xa linear additive generator, "long state=seed1; state += seed2",
Xcapable of generating any 32-bit constant or any arithmetic sequence.
X
X.I mrandom_alg
Xis the algorithm to be used by
X.I mrandomrv() when
Xcalled with the RNG.
X
X.LP
XThe seed and count parameters are used for initialization and
X``cycling'' of the generator before control is returned to the calling
Xprogram.  The generator will be cycled count1+(1.0e9)*count2 times, so
Xbeware: if count2 is non-zero, the underlying RNG will be called
Xbillions of times before control is returned to the calling program!
XMost generators take just one 32-bit seed, but the
X.I seed
Xparameters points to a vector of seeds, whose length is determined by
Xthe needs of the underlying generator.
X
X.I init_rng()
Xreturns a pointer to the allocated and initialized RNG.
X
X.I bufsize
Xis the size of the RNG's main buffer.  A non-positive value of
X.I bufsize will be interpreted as a value of 1.
X
X.I kill_rng()
Xdestroys 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
X.I not
Xuse an
X.I RNGdata
X
Xpointer which points to an active RNG to store the return value of .I
Xinit_rng().  In order to initialize an RNG, you should either declare a
Xnew
X.I RNGdata
Xpointer, and then use it to store the return value of
X.I init_rng(),
Xor, use
X.I kill_rng()
Xto de-initialize an
X.I RNGdata
Xpointer which points to an active RNG, and then use that pointer to
Xstore the return value of
X.I init_rng()
X
XIn general, I believe that users should extend the sequence of an existing RNG,
Xwhenever possible, instead of seeding a new one.
XI suggest this methodology because it is so difficult to properly seed
Xan RNG when performing multiple program runs during a single experiment.
XWhere, after all, can you get truly random seeds?
XTo use the time of day, or a process ID, is to invite disaster in the form
Xof subtle experimental correlations or the (often not-so-subtle) effects of
Xinadvertent reuse of a seed.
XNote that if you use a ``perfectly random'' RNG to generate seeds
Xuniformly distributed in 0..2^{31}-1,
Xyou will are very likely to reuse a seed -- and thus risk a duplicated
Xprogram run -- after running, say, thirty thousand experiments.
X
XAccordingly, my package makes it easy to save and restart RNGs,
Xto minimize the attraction of reseeding with init_rng() calls.
XThe call
X.I save_rng(rng, filename)
Xwill write a complete state table to the specified file.
XThe files are in human-readable ASCII,
Xand are never more than about 1000 characters.
XThe return value of save_rng is 1 if the file is successfully created;
X0 otherwise.
X
XThe call
X.I restart_rng(filename)
Xreads the specified file
Xinto and returns a pointer to the RNG constructed from the file.
XA null pointer is returned, and a message is printed to stderr,
Xif the restart fails due to a garbled or unreadable statefile.
XOtherwise restart_rng returns the value 1.
X
XA major advantage of using restart_rng
Xinstead of init_rng
Xis that the time required to initialize the RNG is independent
Xof the value of the count1, count2 parameters.
X
X.LP
XSeveral routines are available for generating pseudorandom numbers.
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.  For detailed information about
Xbuffering, seed mrandom.tex, included in this distribution.
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
Xwith an invalid RNGdata pointer.
X
XThe two routines for generating bits deserve some extra attention.
X.I bxrandomrv()
Xand
X.I brandomrv()
Xeach use one generate from the RNG to generate each bit.
X.I bxrandomrv_f()
Xand
X.I brandomrv_f()
Xuse each generate to produce 32 bits.  These two routines can only be
Xused with 32-bit generators; they return -1 otherwise.
X
X.I mrandomrv()
Xfills the vector v with n generates in the range 0..m-1.  If
Xrange_rng(rng) < m, the program aborts with an error.
X
XThe algorithm used by
X.I mrandomrv()
Xto fill v is set by
X.I init_rng
Xor by
X.I mralg_rng.
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
X.I mrandomrv()
X, however, this code is slower than algorithm 1, due to the
Xnecessity of recomputing 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
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
X.SH NOTES
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
X.LP
X.I seed_rng()
Xseeds rng with the seed table pointed to by seed.  The RNG's counter is
Xreset to 0.
X
X.LP
X.I check_rng()
Xchecks the integrity of the RNG, in order to determine whether it can be
Xused by the other mrandom library routines.
X
X.LP
X.I describe_rng()
Xplaces a human-readable description of rng in the string rngid.
X
X.LP
X.I mralg_rng()
Xsets the mrandom algorithm number of rng.
X
X.LP
X.I split_rng()
Xsets the split value of rng.
X
X.LP range_rng()
Xreturns the range of rng.
X
X.SH AUTHOR
XRobert Plotkin, [e-mail] and Clark Thomborson,
[e-mail]
X.SH DIAGNOSTICS
XIf error-checking code in any of the routines discovers a problem, an
Xerror message is printed on the stderr stream.
X
X.SH "SEE ALSO"
Xrandom(3), rand(3C), drand48(3), mrtest(1)
X
X.SH BUGS
XA little slower than
X.I random().
X
XA source-code rewrite of random() would allow efficient, bug-free,
Xstate-saving and restarting.  As things stand, you can get
X"non-random" and "non-restartable" RNG streams by calling
Xinit_rng() several times with alg=1.  Perhaps I should add code to
Xgenerate an error message in this case.
X
X.SH THEORY
XAside from the bug with multiple simultaneous generators,
XI know of no reason to choose
X4.3bsd random() over the Knuth/Bentley RNG.
XBoth are likely to fail Marsaglia's Birthday-Spacings test,
Xalthough I don't know that this has been tested (and I don't see
Xthat this is likely to pose a problem in any ``real'' application).
XThe portable-combined RNG is noticeably slower, but arguably superior
Xto both the above.
XThe defects of 4.3bsd rand() and nrandom() are fairly well-known,
Xand they are only included in this package for reasons of
Xbackward compatibility and testing.
X
XFor more information on this package, see the LaTex files mrandom.tex
Xand soda.tex included with this distribution.
X
X.SH REFERENCES
X
XJon Louis Bentley,
X``The software exploratorium: Some random thoughts.''
X.I UNIX Review 10,
XNumber 6, June 1992.
X
XPierre L'Ecuyer,
X``Efficient and portable combined random number generators.''
X.I Communications of the ACM, 31(6):
X742--774, June 1988.
X
XGeorge Marsaglia,
X``A current view of random number generators.''
XIn L. Billard, editor,
X.I Computer Science and Statistics: The Interface,
Xpages 3--10. Elsevier Science Publishers, 1985.
X
XGeorge Marsaglia and Arif Zaman,
X``A New Class of Random Number Generators.''
X.I The Annals of Applied Probability, 1(3):
X1991.
X
XOra E. Percus and Malvin H. Kalos,
X``Random number generators for MIMD parallel processors.''
X.I Journal of Parallel and Distributed Computing, 6:
X477--497, 1989.
X
XWilliam H. Press and Saul A. Teukolsky,
X``Portable Random Number Generators.''
X.I Computers in Physics, 6(5):
XSept/Oct. 1992.
X
XRobert Sedgewick,
X.I Algorithms in C,
XAddison-Wesley, 1990.
X
XClark Thomborson, "Tools for Randomized Experimentation," to appear in the
X.I Proceedings of the 25th Symposium on the Interface:
X.I Computing Science and Statistics,
X1993.
X
XClark Thomborson.
X``Mrandom (version 1).''
X.I Comp.sources.unix 25(23),
XDecember 1991.
X
END_OF_FILE
if test 21315 -ne `wc -c <'man/mrandom.3'`; then
    echo shar: \"'man/mrandom.3'\" unpacked with wrong size!
fi
# end of 'man/mrandom.3'
fi
if test -f 'man/mrandom.3.txt' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'man/mrandom.3.txt'\"
else
echo shar: Extracting \"'man/mrandom.3.txt'\" \(21353 characters\)
sed "s/^X//" >'man/mrandom.3.txt' <<'END_OF_FILE'
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
XNAME
X     init_rng, kill_rng, dxrandomrv, dxrandomr, dxrandomv, dxran-
X     dom, fxrandomrv, fxrandomr, fxrandomv, fxrandom, lxrandomrv,
X     lxrandomr,  lxrandomv,  lxrandom,   bxrandomrv,   bxrandomr,
X     bxrandomv, bxrandom, bxrandomrv_f, bxrandomr_f, bxrandomv_f,
X     bxrandom_f, drandomrv, drandomr,  drandomv,  drandom,  fran-
X     domrv,  frandomr,  frandomv,  frandom,  lrandomrv, lrandomr,
X     lrandomv, lrandom, brandomrv, brandomr,  brandomv,  brandom,
X     brandomrv_f,  brandomr_f,  brandomv_f, brandom_f, mrandomrv,
X     mrandomr,   mrandomv,    mrandom,    flush_rng,    save_rng,
X     restart_rng,  seed_rng,  check_rng, describe_rng, mralg_rng,
X     split_rng, range_rng - a uniform interface to several random
X     number generators
X
XSYNOPSIS
X     RNGdata *init_rng(alg, mrandom_alg, seed, count1, count2,
X     long alg, mrandom_alg, *seed, count1, count2,
X
X     int kill_rng(rng)
X     RNGdata *rng;
X
X     double dxrandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     double v[];
X
X     float fxrandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     float v[];
X
X     long lxrandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     long v[];
X
X     long lxrandomv(n, v)
X     long n;
X     long v[];
X
X     int bxrandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     int v[];
X
X     int bxrandomrv_f(rng, n, v)
X     RNGdata *rng;
X     long n;
X     int v[];
X
X     double drandomrv(rng, n, v)
X     RNGdata *rng;
X
X
X
XSun Release 4.1       Last change: 5/28/93                      1
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     long n;
X     double v[];
X
X     float frandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     float v[];
X
X     long lrandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     long v[];
X
X     int brandomrv(rng, n, v)
X     RNGdata *rng;
X     long n;
X     int v[];
X
X     int brandomrv_f(rng, n, v)
X     RNGdata *rng;
X     long n;
X     int v[];
X
X     long mrandomrv(rng, m, n, v)
X     RNGdata *rng;
X     long m;
X     long n;
X     long v[];
X
X     int save_rng(rng, filename)
X     RNGdata *rng;
X     char *filename;
X
X     RNGdata *restart_rng(rng, filename)
X     char *filename;
X
X     void seed_rng(rng, seed)
X     RNGdata *rng;
X     long *seed;
X
X     int check_rng(rng);
X     RNGdata *rng;
X
X     char *describe_rng(rng, rngid)
X     RNGdata *rng;
X     char rngid[RNGIDSTRLEN];
X
X     int mralg_rng(rng,new_value)
X     RNGdata *rng;
X     long new_value;
X
X
X
X
X
XSun Release 4.1       Last change: 5/28/93                      2
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
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
XDESCRIPTION
X     _m_r_a_n_d_o_m_r_v(_r_n_g,_m,_n,_v) generates  n  random  integers  in  the
X     range 0 to m-1, using the generator rng, and storing them in
X     vector v.  An initialization routine, _i_n_i_t__r_n_g(),  described
X     below, allows you to select an appropriate generator.
X
X     Using this package, instead of direct calls to  random()  or
X     some other RNG code, offers the following advantages.
X
X     1.   You can switch to another RNG, and re-run a portion  of
X          your experiment to check its validity, without changing
X          any of your code other than a single parameter in  your
X          _i_n_i_t__r_n_g() call.
X
X     2.   You can use my unbiased method  for  generating  random
X          integers in the range 0..m-1.  By contrast, the typical
X          integer-generating codes "random()%m"  or  "(int)  m  *
X          ((double)  random()/MAXLONG)"  have  a  slight bias for
X          large m.
X
X     3.   Floating point numbers, uniformly distributed in  [0.0,
X          1.0), are available (as on most RNG interfaces).
X
X     4.   As  with  most  RNG  interfaces,  you  can  have   many
X          simultaneously-active   RNGs.   (Unfortunately,  I  was
X          unable to find  an  efficient  work-around  for  4.3bsd
X          random()'s  state-saving bug; you'll have to use one of
X          the other generators if you want to  use  multiple  RNG
X          streams in a single program invocation.)
X
X     5.   Time-efficient  vectorized  calls,  returning  multiple
X          uniform variates, are available.
X
X     6    The ability to "split" RNGs to produce parallet  output
X          streams using the "leapfrog" method.
X
X     7.   The package offers a shorthand notation for  completely
X          specifying  the  algorithm  and  current  state of your
X          RNG(s), in an 80-character human-readable ASCII string.
X          This  is  very useful in experimental documentation and
X          replication.
X
X     8.   A complete RNG state can  be  (serially)  reconstructed
X          from its shorthand notation, and
X
X
X
X
XSun Release 4.1       Last change: 5/28/93                      3
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     9.   A file-I/O interface allows fast saves and restarts  of
X          complete  RNG  state vectors, without the time overhead
X          of serial reconstruction.
X
X
X     Also included in this software release  is  a  (unsupported)
X     driver routine _m_r_t_e_s_t._c. This routine implements some of the
X     simpler tests of randomness, e.g. equidistribution, pairwise
X     (both short- and long-range) correlation, and 3-tuple corre-
X     lation [Marsaglia, 1985].  These tests were chosen to illus-
X     trate the use of the components of the _m_r_a_n_d_o_m() package, as
X     well as to exhibit the known defects of the  4.3bsd  genera-
X     tors  rand()  and nrand48().  See the man page on _m_r_t_e_s_t for
X     more details.
X
X     Another use for _m_r_t_e_s_t._c is as a template for your own calls
X     to  the  mrandom() package.  (I, for one, like to program by
X     example...)
X
X     Below are detailed descriptions of  the  interfaces  to  the
X     various  C-language functions in the mrandom() package.  For
X     technical  details,  consult  the  file  _m_r_a_n_d_o_m._t_e_x,   also
X     included in this distribution.
X
X     In order to use an RNG, it must first be initialized.   This
X     is  accomplished  by first declaring a pointer to an RNGdata
X     structure, and then calling  _i_n_i_t__r_n_g()  ,  which  allocates
X     memory  for the RNG and readies it for use by the other rou-
X     tines in the package.
X
X     _i_n_i_t__r_n_g() returns a  pointer  to  an  initialized  RNG.   A
X     pointer returned by _i_n_i_t__r_n_g() _i_s _v_a_l_i_d _f_o_r _u_s_e _b_y
X
X     In an init_rng() call, the _a_l_g parameter should have a value
X     between 0 and 9, to indicate which RNG algorithm is desired:
X
X     1.   4.3bsd random(), a non-linear additive feedback RNG,
X
X     2.   The Knuth/Bentley prand(), a lagged-Fibonnacci  genera-
X          tor with a state table of size 55 [Bentley, 1992],
X
X     3.   4.3bsd nrand48(), a 48-bit multiplicative  congruential
X          RNG,
X
X     4.   L'Ecuyer's ``portable combined'' 32-bit  multiplicative
X          congruential generator [L'Ecuyer, 1988],
X
X     5.   4.3bsd rand(), the  discredited  32-bit  multiplicative
X          congruential RNG, and
X
X     6-8. Press & Teukolsky's ran0, ran1, and ran2,  respectively
X          [Press & Teukolsky, 1992],
X
X
X
XSun Release 4.1       Last change: 5/28/93                      4
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     9.   Marsaglia's subtract-with-borrow Ultra generator  [Mar-
X          saglia and Zaman, 1991],
X
X     0.   for testing  purposes,  a  linear  additive  generator,
X          "long state=seed1; state += seed2", capable of generat-
X          ing any 32-bit constant or any arithmetic sequence.
X
X          _m_r_a_n_d_o_m__a_l_g is the algorithm to be used by  _m_r_a_n_d_o_m_r_v()
X          _w_h_e_n called with the RNG.
X
X
X     The seed and count parameters are  used  for  initialization
X     and  ``cycling'' of the generator before control is returned
X     to the  calling  program.   The  generator  will  be  cycled
X     count1+(1.0e9)*count2  times,  so  beware: if count2 is non-
X     zero, the underlying RNG will be called  billions  of  times
X     before  control  is  returned  to the calling program!  Most
X     generators take just one 32-bit seed, but the  _s_e_e_d  parame-
X     ters points to a vector of seeds, whose length is determined
X     by the needs of the underlying generator.
X
X     _i_n_i_t__r_n_g() returns a pointer to the allocated  and  initial-
X     ized RNG.
X
X     _b_u_f_s_i_z_e is the size  of  the  RNG's  main  buffer.   A  non-
X     positive value of _b_u_f_s_i_z_e _w_i_l_l _b_e _i_n_t_e_r_p_r_e_t_e_d _a_s _a
X
X     _k_i_l_l__r_n_g() destroys the RNG,  making  it  invalid  for  use.
X     This  procedure  de-allocates the space used by the RNG, and
X     should therefore be used to kill RNGs which will  no  longer
X     be used.
X
X     Do _n_o_t use an _R_N_G_d_a_t_a
X
X     pointer which points to an active RNG to  store  the  return
X     value  of .I init_rng().  In order to initialize an RNG, you
X     should either declare a new _R_N_G_d_a_t_a pointer, and then use it
X     to  store the return value of _i_n_i_t__r_n_g(), or, use _k_i_l_l__r_n_g()
X     to de-initialize an  _R_N_G_d_a_t_a  pointer  which  points  to  an
X     active  RNG,  and  then use that pointer to store the return
X     value of _i_n_i_t__r_n_g()
X
X     In general, I believe that users should extend the  sequence
X     of  an existing RNG, whenever possible, instead of seeding a
X     new one.  I suggest this methodology because it is so diffi-
X     cult  to  properly seed an RNG when performing multiple pro-
X     gram runs during a single experiment.  Where, after all, can
X     you  get  truly  random seeds?  To use the time of day, or a
X     process ID, is to invite disaster  in  the  form  of  subtle
X     experimental   correlations  or  the  (often  not-so-subtle)
X     effects of inadvertent reuse of a seed.  Note  that  if  you
X     use  a  ``perfectly random'' RNG to generate seeds uniformly
X
X
X
XSun Release 4.1       Last change: 5/28/93                      5
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     distributed in 0..2^{31}-1, you  will  are  very  likely  to
X     reuse  a  seed  -- and thus risk a duplicated program run --
X     after running, say, thirty thousand experiments.
X
X     Accordingly, my package makes it easy to  save  and  restart
X     RNGs,   to   minimize   the  attraction  of  reseeding  with
X     init_rng() calls.  The  call  _s_a_v_e__r_n_g(_r_n_g,  _f_i_l_e_n_a_m_e)  will
X     write  a  complete  state  table to the specified file.  The
X     files are in human-readable ASCII, and are never  more  than
X     about 1000 characters.  The return value of save_rng is 1 if
X     the file is successfully created; 0 otherwise.
X
X     The call _r_e_s_t_a_r_t__r_n_g(_f_i_l_e_n_a_m_e) reads the specified file into
X     and  returns a pointer to the RNG constructed from the file.
X     A null pointer is returned, and  a  message  is  printed  to
X     stderr,  if the restart fails due to a garbled or unreadable
X     statefile.  Otherwise restart_rng returns the value 1.
X
X     A major advantage of using restart_rng instead  of  init_rng
X     is  that the time required to initialize the RNG is indepen-
X     dent of the value of the count1, count2 parameters.
X
X
X     Several routines are available for  generating  pseudorandom
X     numbers.   Both  buffered  and  unbuffered routines are pro-
X     vided.  Unbuffered routines call the underlying RNG only  as
X     many  times as are needed to produce the requested number of
X     generates, while buffered routines maintain buffers of  gen-
X     erates,  so  that generates may be produced efficiently even
X     when requested in small quantities.  Roughly, buffered  rou-
X     tines  are  preferable when generates are requested one at a
X     time or in small quantities, while unbuffered  routines  are
X     preferable when generates are requested in large quantities.
X     For detailed information about buffering, seed  mrandom.tex,
X     included in this distribution.
X
X     The name of a routine denotes the type of  the  value  which
X     the  routine  returns and whether the routine is buffered or
X     unbuffered. The first letter of a routine denotes  the  type
X     of  value  which  it returns: ``d'' for double precision and
X     ``f'' for single  precision  floating  point  in  the  range
X     [0,1);    ``l''    for    long    integer   in   the   range
X     0..(range_rng(rng)-1), and ``b'' for bit (either a  0  or  a
X     1).  If the second letter of the routine's name is an ``x'',
X     then the routine is unbuffered.  Otherwise, the  routine  is
X     buffered.
X
X     For convenience in  user  programming,  we  also  provide  a
X     number  of macros that supply default parameter values.  The
X     last two letters of all our fundamental routines is  ``rv''.
X     This means that they must be provided with both a pointer to
X     an RNGdata structure and a vector  to  fill  with  generates
X
X
X
XSun Release 4.1       Last change: 5/28/93                      6
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     from  the  RNG.   Macros whose names do not contain an ``r''
X     have the RNGdata pointer omitted from their parameter  list;
X     they  use  the most-recently initialized or restarted RNG to
X     produce generates.  Macros whose  names  do  not  contain  a
X     ``v''  have  the vector and number of generates omitted from
X     their parameter list; they produce and return a single  gen-
X     erate.
X
X     All generating routines abort with a message  to  stderr  if
X     called with an invalid RNGdata pointer.
X
X     The two routines for  generating  bits  deserve  some  extra
X     attention.   _b_x_r_a_n_d_o_m_r_v()  and _b_r_a_n_d_o_m_r_v() each use one gen-
X     erate from the RNG to generate each bit.  _b_x_r_a_n_d_o_m_r_v__f() and
X     _b_r_a_n_d_o_m_r_v__f()  use  each generate to produce 32 bits.  These
X     two routines can only be used with 32-bit  generators;  they
X     return -1 otherwise.
X
X     _m_r_a_n_d_o_m_r_v() fills the vector v with n generates in the range
X     0..m-1.   If  range_rng(rng) < m, the program aborts with an
X     error.
X
X     The algorithm used by  _m_r_a_n_d_o_m_r_v()  to  fill  v  is  set  by
X     _i_n_i_t__r_n_g or by _m_r_a_l_g__r_n_g.
X
X     Algorithm 0 is Thomborson's unbiased method, which  produces
X     unbiased  long  integers in the range [0..m).  The algorithm
X     discards any outputs from rng which are larger than r-(r mod
X     m), where r is equal to range_rng(rng).  At worst, this code
X     will discard (on long-term average) at most one value  of  r
X     for  every  one  that  is  useful.   This worst case is only
X     encountered for extremely large m; for fixed and moderate m,
X     this  code  will  rarely  discard a value, and thus will run
X     essentially as fast as algorithm 1.  When  the  value  of  m
X     changes  on each call to _m_r_a_n_d_o_m_r_v() , however, this code is
X     slower than algorithm 1, due to the necessity of recomputing
X     r-(r mod m).
X
X     The program aborts with an error message to stderr if rng is
X     behaving  so  non-randomly  that  Algorithm  0  must make an
X     excessive number of calls to rng in  order  to  produce  the
X     requested number of generates.
X
X     Algorithm 1 is the standard (long)(m*dxrandomr(rng)).   This
X     algorithm may be biased: for large m, some results may be be
X     more likely than others.  The bias is (r mod m)/m, which  is
X     upper-bounded  by  0.1%  if m is less than a million and the
X     range r of rng is at least a billion.
X
X     We do not support, and indeed we actively  discourage,  gen-
X     erating   restricted-range  integers  with  lrandomr(rng)%m.
X     Many RNGs have poor behavior under this transformation, most
X
X
X
XSun Release 4.1       Last change: 5/28/93                      7
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     noticeably when m is a power of 2.  When m is not a power of
X     2, fixed-point division required by an  ``%''  operation  is
X     time-consuming on many workstations.
X
X
XNOTES
X     The  mrandomrv  procedure  is  capable  of  generating  long
X     integers  in  the  full  range  of  any  RNG  for which 1 <=
X     range_rng(rng) <= 2^32.  In order to accomplish  this,  with
X     the parameter m a signed long integer, the following mapping
X     is used:
X
X     Range(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
X
X     _s_e_e_d__r_n_g() seeds rng with the seed table pointed to by seed.
X     The RNG's counter is reset to 0.
X
X
X     _c_h_e_c_k__r_n_g() checks the integrity of the  RNG,  in  order  to
X     determine  whether  it  can  be  used  by  the other mrandom
X     library routines.
X
X
X     _d_e_s_c_r_i_b_e__r_n_g() places a human-readable description of rng in
X     the string rngid.
X
X
X     _m_r_a_l_g__r_n_g() sets the mrandom algorithm number of rng.
X
X
X     _s_p_l_i_t__r_n_g() sets the split value of rng.
X
X
X     returns the range of rng.
X
X
XAUTHOR
X     Robert Plotkin, [e-mail] and  Clark  Thombor-
X     son, [e-mail]
X
XDIAGNOSTICS
X     If error-checking code in any of the  routines  discovers  a
X     problem, an error message is printed on the stderr stream.
X
X
XSEE ALSO
X     random(3), rand(3C), drand48(3), mrtest(1)
X
X
X
X
X
XSun Release 4.1       Last change: 5/28/93                      8
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
XBUGS
X     A little slower than _r_a_n_d_o_m().
X
X     A source-code rewrite of  random()  would  allow  efficient,
X     bug-free, state-saving and restarting.  As things stand, you
X     can get "non-random" and "non-restartable"  RNG  streams  by
X     calling  init_rng()  several  times  with  alg=1.  Perhaps I
X     should add code to generate an error message in this case.
X
X
XTHEORY
X     Aside from the bug with multiple simultaneous generators,  I
X     know  of  no  reason  to  choose  4.3bsd  random()  over the
X     Knuth/Bentley RNG.  Both  are  likely  to  fail  Marsaglia's
X     Birthday-Spacings  test, although I don't know that this has
X     been tested (and I don't see that this is likely to  pose  a
X     problem in any ``real'' application).  The portable-combined
X     RNG is noticeably slower, but arguably superior to both  the
X     above.   The  defects  of  4.3bsd  rand()  and nrandom() are
X     fairly well-known, and they are only included in this  pack-
X     age for reasons of backward compatibility and testing.
X
X     For more information on this package, see  the  LaTex  files
X     mrandom.tex and soda.tex included with this distribution.
X
X
XREFERENCES
X     Jon Louis Bentley, ``The software exploratorium: Some random
X     thoughts.'' _U_N_I_X _R_e_v_i_e_w _1_0, Number 6, June 1992.
X
X     Pierre L'Ecuyer, ``Efficient and  portable  combined  random
X     number  generators.'' _C_o_m_m_u_n_i_c_a_t_i_o_n_s _o_f _t_h_e _A_C_M, _3_1(_6): 742-
X     -774, June 1988.
X
X     George Marsaglia, ``A current view of random number  genera-
X     tors.''  In L. Billard, editor, _C_o_m_p_u_t_e_r _S_c_i_e_n_c_e _a_n_d _S_t_a_t_i_s_-
X     _t_i_c_s: _T_h_e _I_n_t_e_r_f_a_c_e, pages 3--10. Elsevier Science  Publish-
X     ers, 1985.
X
X     George Marsaglia and Arif Zaman, ``A  New  Class  of  Random
X     Number  Generators.''  _T_h_e  _A_n_n_a_l_s  _o_f  _A_p_p_l_i_e_d _P_r_o_b_a_b_i_l_i_t_y,
X     _1(_3): 1991.
X
X     Ora E. Percus and Malvin H. Kalos, ``Random  number  genera-
X     tors for MIMD parallel processors.'' _J_o_u_r_n_a_l _o_f _P_a_r_a_l_l_e_l _a_n_d
X     _D_i_s_t_r_i_b_u_t_e_d _C_o_m_p_u_t_i_n_g, 477--497, 1989.
X
X     William H. Press and Saul A.  Teukolsky,  ``Portable  Random
X     Number  Generators.''  _C_o_m_p_u_t_e_r_s _i_n _P_h_y_s_i_c_s, _6(_5): Sept/Oct.
X     1992.
X
X     Robert Sedgewick, _A_l_g_o_r_i_t_h_m_s _i_n _C, Addison-Wesley, 1990.
X
X
X
XSun Release 4.1       Last change: 5/28/93                      9
X
X
X
X
X
X
XMRANDOM(3)             C LIBRARY FUNCTIONS             MRANDOM(3)
X
X
X
X     Clark Thomborson, "Tools for Randomized Experimentation," to
X     appear in the _P_r_o_c_e_e_d_i_n_g_s _o_f _t_h_e _2_5_t_h _S_y_m_p_o_s_i_u_m _o_n _C_o_m_p_u_t_i_n_g
X     _S_c_i_e_n_c_e _a_n_d _S_t_a_t_i_s_t_i_c_s, 1993.
X
X     Clark     Thomborson.      ``Mrandom     (version      1).''
X     _C_o_m_p._s_o_u_r_c_e_s._u_n_i_x _2_5(_2_3), December 1991.
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XSun Release 4.1       Last change: 5/28/93                     10
X
X
X
END_OF_FILE
echo shar: 652 control characters may be missing from \"'man/mrandom.3.txt'\"
if test 21353 -ne `wc -c <'man/mrandom.3.txt'`; then
    echo shar: \"'man/mrandom.3.txt'\" unpacked with wrong size!
fi
# end of 'man/mrandom.3.txt'
fi
if test -f 'src/mrandom.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/mrandom.h'\"
else
echo shar: Extracting \"'src/mrandom.h'\" \(20309 characters\)
sed "s/^X//" >'src/mrandom.h' <<'END_OF_FILE'
X/* mrandom.h 3.1 5/28/93 */
X/*
X *        mrandom.h
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#define MRANDOM
X#endif
X
X/*****************************/
X/* Miscellaneous definitions */
X/*****************************/
X#define RNGIDSTRLEN 80 /* Maximum length of string written by */
X               /* describe_rng() */
X#define BILLION    1000000000 /* The modulus for 32-bit RNGcount1 in our */
X               /* statefile, i.e. */
X                           /* (Number of calls to RNG since seeding) = */
X               /* RNGcount1 + BILLION * RNGcount2 */
X#define BBUF_SIZE    32 /* Size of the bit buffer */
X#define NUM_RNGS     10 /* Number of RNGs currently installed */
X#define SPLIT_DEF    0  /* Default split value for all RNGs */
X
X#define RET_LONG     4  /* Types of values returned by */
X#define RET_DOUBLE   -8 /* RNG generating routines */
X
X#define STATE_CHAR   1  /* Types of values stored in */
X#define STATE_INT    2  /* RNG state and seed vectors */
X#define STATE_LONG   4
X#define STATE_FLOAT  -4 /* STATE_FLOAT and STATE_DOUBLE are */
X#define STATE_DOUBLE -8 /* currently unsupported */
X
X/**************************/
X/* RNGdata                */
X/* Data describing an RNG */
X/**************************/
Xstruct rngdata {
X  /* The following data is taken directly from the statefile */
X  long    rngalg;        /* algorithm used to generate pseudorandoms:
X             * 0 = trivial rng (long state=seed1; state += seed2)
X             * 1 = 4.3bsd random.c,
X             * 2 = bentley.c,
X             * 3 = pcrand.c,
X             * 4 = 4.3bsd nrand48.c,
X             * 5 = 4.3bsd rand.c,
X             * 6 = Press & Teukolsky's ran0 (ran0.c)
X             * 7 = Press & Teukolsky's ran1 (ran1.c)
X             * 8 = Press & Teukolsky's ran2 (ran2.c)
X             * 9 = Marsaglia's Ultra (ultra.c)
X             */
X  long mrandom_alg; /* Algorithm used by mrandomrv */
X  long *rngstate;   /* The RNG's state vector */
X  long *rngseed;    /* The seed originally used to initialize rngstate */
X  long rngcount1;   /* mod-BILLION counter of calls to the RNG */
X  long rngcount2;   /* div-BILLION counter */
X  struct {
X    long size;              /* Number of entries in buffer */
X    long nleft;             /* Number of values left in buffer */
X    int nbleft;             /* Number of values in bit buffer */
X    double *dbuf,*dbufp;    /* RNG double buffer [0,1) */
X    long *lbuf,*lbufp;      /* RNG long buffer 0...RNGrange-1 */
X    int *bbuf,*bbufp;       /* RNG bit buffer 0,1 */
X  } buffer;
X  long rngnextval;          /* the RNG's next output */
X  long rngsplit;            /* Return every rngsplit-th generate */
X
X  /* The next six values are inferred from rngalg */
X  char rngname[RNGIDSTRLEN]; /* Name of the RNG */
X  long rngstatesize;
X  long rngseedsize; /* It is permissible to have an RNG which can use */
X            /* different numbers of seeds.  In that case, use */
X            /* the maximum permissible number of seeds for */
X            /* this value. */
X  double rngrange;  /* The RNG's range, expressed as a d.p. float */
X  signed int rngreturns; /* The kind of value the rng returns
X                - Positive numbers indicate integers in the
X                  range 0...RNGrange-1
X                - Negative numbers indicate floats in the
X                  range [0,1)
X                - Absolute value indicates the number of
X                  bytes per RNG output
X                  e.g. 4 means a 32-bit integer,
X                  -16 means a 128 bit float
X                  Currently only 4 and -8 are supported.
X                  These are defined as RET_LONG and
X                  RET_DOUBLE, respectively. */
X  int rngstatetype; /* The type of the value stored in the RNG's state */
X            /* and seed vectors.  Currently, three types are */
X            /* supported:
X               Type            Value
X               -----------     ----------
X               8-bit char      STATE_CHAR
X               16-bit int      STATE_INT
X               32-bit long     STATE_LONG
X               */
X};
Xtypedef struct rngdata RNGdata;
X
X/* Convenient names for data in our RNGdata struct */
X#define RNGalg        (rng->rngalg)
X#define RNGname         (rng->rngname)
X#define RNGseed            (rng->rngseed)
X#define RNGcount1    (rng->rngcount1)
X#define RNGcount2    (rng->rngcount2)
X#define RNGnextval    (rng->rngnextval)
X#define RNGstatesize    (rng->rngstatesize)
X#define RNGseedsize     (rng->rngseedsize)
X#define RNGrangem1    (rng->rngrangem1)
X#define RNGrange    (rng->rngrange)
X#define RNGstate    (rng->rngstate)
X#define RNGmrandom_alg  (rng->mrandom_alg)
X#define RNGbuffer       (rng->buffer)
X#define RNGreturns      (rng->rngreturns)
X#define RNGstatetype    (rng->rngstatetype)
X#define RNGsplit        (rng->rngsplit)
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(/*RNGdata *rng, long n, double v[]*/);
X#define drandomr(rng)  drandomrv(rng,0,0)
X#define drandomv(n,v)  drandomrv(0,n,v)
X#define drandom()      drandomrv(0,0,0)
X
X/******************/
X/* frandomrv      */
X/* Returns floats */
X/******************/
Xfloat frandomrv(/*RNGdata *rng, long n, float v[]*/);
X#define frandomr(rng)  frandomrv(rng,0,0)
X#define frandomv(n,v)  frandomrv(0,n,v)
X#define frandom()      frandomrv(0,0,0)
X
X/*****************/
X/* lrandomrv     */
X/* Returns longs */
X/*****************/
Xlong lrandomrv(/*RNGdata *rng, long n, long v[]*/);
X#define lrandomr(rng)  lrandomrv(rng,0,0)
X#define lrandomv(n,v)  lrandomrv(0,n,v)
X#define lrandom()      lrandomrv(0,0,0)
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(/*RNGdata *rng, long n, int v[]*/);
X#define brandomr(rng)  brandomrv(rng,0,0)
X#define brandomv(n,v)  brandomrv(0,n,v)
X#define brandom()      brandomrv(0,0,0)
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(/*RNGdata *rng, long n, int v[]*/);
X#define brandomr_f(rng)  brandomrv_f(rng,0,0)
X#define brandomv_f(n,v)  brandomrv_f(0,n,v)
X#define brandom_f()      brandomrv_f(0,0,0)
X
X/*********************************/
X/* Unbuffered calling procedures */
X/*********************************/
X/*******************/
X/* dxrandomrv      */
X/* Returns doubles */
X/*******************/
Xdouble dxrandomrv(/*RNGdata *rng, long n, double v[]*/);
X#define dxrandomr(rng) dxrandomrv(rng,0,0)
X#define dxrandomv(n,v) dxrandomrv(0,n,v)
X#define dxrandom()     dxrandomrv(0,0,0)
X
X/******************/
X/* fxrandomrv     */
X/* Returns floats */
X/******************/
Xfloat fxrandomrv(/*RNGdata *rng, long n, float v[]*/);
X#define fxrandomr(rng) fxrandomrv(rng,0,0)
X#define fxrandomv(n,v) fxrandomrv(0,n,v)
X#define fxrandom()     fxrandomrv(0,0,0)
X
X/*****************/
X/* lxrandomrv    */
X/* Returns longs */
X/*****************/
Xlong lxrandomrv(/*RNGdata *rng, long n, long v[]*/);
X#define lxrandomr(rng) lxrandomrv(rng,0,0)
X#define lxrandomv(n,v) lxrandomrv(0,n,v)
X#define lxrandom()     lxrandomrv(0,0,0)
X
X/*******************************/
X/* bxrandomrv                  */
X/* Return "high quality" bits. */
X/*******************************/
Xint bxrandomrv(/*RNGdata *rng, long n, int v[]*/);
X#define bxrandomr(rng) bxrandomrv(rng,0,0)
X#define bxrandomv(n,v) bxrandomrv(0,n,v)
X#define bxrandom()     bxrandomrv(0,0,0)
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(/*RNGdata *rng, long n, int v[]*/);
X#define bxrandomr_f(rng) bxrandomrv_f(rng,0,0)
X#define bxrandomv_f(n,v) bxrandomrv_f(0,n,v)
X#define bxrandom_f()     bxrandomrv_f(0,0,0)
X
X/****************************/
X/* seed_rng                 */
X/* Seeds the underlying RNG */
X/****************************/
Xvoid seed_rng(/*RNGdata *rng, long *seed*/);
X
X/*************************************/
X/* check_rng                         */
X/* Check the integrity of the RNG.   */
X/* Return 1 if ok, 0 if not.         */
X/*************************************/
Xint check_rng(/*RNGdata *rng*/);
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 (/*RNGdata *rng, char rngid[RNGIDSTRLEN]*/);
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(/* RNGdata *rng, long new_value*/);
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(/* RNGdata *rng, long new_value*/);
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(/*RNGdata *rng*/);
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 L'Ecuyer's 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(/* long alg, long mrandom_alg, long *seed,
X             long count1, long count2, long bufsize*/);
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(/*RNGdata *rng*/);
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(/* RNGdata *rng, char *filename*/);
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(/* char *filename*/);
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(/*RNGdata *rng*/);
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 (/* RNGdata *rng, long m, long n, long v*/);
X#define mrandomr(rng,m)  mrandomrv(rng,m,0,0)
X#define mrandomv(m,n,v)  mrandomrv(0,m,n,v)
X#define mrandom(m)       mrandomrv(0,m,0,0)
END_OF_FILE
if test 20309 -ne `wc -c <'src/mrandom.h'`; then
    echo shar: \"'src/mrandom.h'\" unpacked with wrong size!
fi
# end of 'src/mrandom.h'
fi
if test -f 'src/mrtest.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/mrtest.c'\"
else
echo shar: Extracting \"'src/mrtest.c'\" \(21361 characters\)
sed "s/^X//" >'src/mrtest.c' <<'END_OF_FILE'
X/* mrtest.c 3.1 5/28/93 */
X/*
X *        mrtest.c
X *
X *    Test routine for mrandom.c
X *
X *    Original Implementation:
X *        Clark Thomborson, September 1991
X *
X *    Modifications:
X *
X *    Clark Thomborson, May 1992 -- used new mrandom.h features,
X *        -added Chi-square test (based on Sedgewick's {it Algorithms}),
X *        -added -dnnn option, to demonstrate a problem with nrand48,
X *        -added Marsaglia's k-tuple test (from {\it Computer Science and
X *         Statistics: The Interface}, L. Billard (ed.), Elsevier 1985,
X *         pp. 3--10,
X *        -added -v, -f, and -p options; added -DVECTOR cc option
X *        -added median and mod 3 tests
X *        -added -t,-e options; added lots of argerr() calls
X *        -used xsq instead of Chi-square.  Statistic is not chi-squared,
X *         except in large-n limit.
X *        -improved accuracy of xsq calcs for very large n
X *        -adjusted code, so "mrtest 0" prints current rngstate
X *        -added large-m equidistribution test
X *        -debugged default RNGfile creation code (Aug 1992)
X *
X *      Robert Plotkin, May 1993
X *              -update for use with mrandom 3.0
X *
X *    Possible future additions:
X *        -more RNGs
X *        -faster (hash-based) large-m equidistribution test
X *        -add Marsaglia's "Birthday-Spacings" test, to exhibit a
X *         shortcoming of Bentley's prand() (as well as 4.3bsd random()):
X *         sort n results from mrandom(m), compute the list of "spacings"
X *         (differences) between adjacent output values in the sorted
X *         list, then let Y be the number of values that appear
X *         more than once among the spacings.  Y is asymptotically
X *         Poisson with parameter $\lambda = n^3/(4m)$ [Marsaglia,
X *         ``A current view of random number generators'', Computer
X *         Science and Statistics: The Interface, Elsevier, 1985;
X *         Marsaglia says the ``proof, due to Janos Komlos, and
X *         detailed discussion of the test will appear elsewhere.''].
X *         According to Marsaglia, lagged-Fibbonacci generators
X *         based on binary addition, e.g. Bentley's RNG, fail this
X *         test.  Quite probably random() would fail it as well.
X *        -add a Poisson interpretation routine, to support the
X *         Birthday Spacings test.  This could also be useful in
X *         a possible improvement to the interpret_xsq routine,
X *         since Poisson statistics could be used to analyze frequency
X *         counts in (our usual) case of very small n/k.
X *
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.
X */
X
X# include <sys/file.h> /* we use access() */
X# include "mrandom.h"
X# include "xsq.h" /* for interpretation of X-squared values */
X# include <math.h>
X# include <values.h> /* we need MAXLONG */
X
X# ifndef RNGFILENAME
X# define RNGFILENAME "RNGstatefile" /* where the RNG state is stored */
X# endif
X
X/* max(m,n) for which equidistribution test will be performed */
X# define EQUIDISTMAX (1<<18)
X/* max m for which 2-tuple correlation will be tested */
X# define PAIRMAXM 256
X
Xvoid argerr(progname)
Xchar *progname;
X{
X  printf(
X    "Usage: %s [ nn -dnnnn -S[n[,n[,n[,n[,n]]]]] -mnn -Mnn -q -t -e]",
X     progname);
X#ifndef VECTOR
X  printf(" -f -v -p");
X#endif
X  printf(" ]\n");
X  printf("  nn sets number of random generates to be tested.  Default 10.\n");
X  printf(
X   "  -dnnnn discards nnnn generates between tested values.  Default 0.\n");
X  printf("  -S[n[,n[,n[,n[,n[,n]]]]]] initializes an RNG, as follows:\n");
X  printf("    parameter #1 sets the RNG algorithm:\n");
X  printf("      0 is an additive linear scheme (for testing only);\n");
X  printf("      1 is 4.3bsd random,\n");
X  printf("      2 is the Knuth/Bentley prand,\n");
X  printf("      3 is L'Ecuyer's portable combined multiplicative RNG,\n");
X  printf("      4 is 4.3bsd nrandom48, and\n");
X  printf("      5 is 4.3bsd rand;\n");
X  printf("      6 is Press & Teukolsky's ran0\n");
X  printf("      7 is Press & Teukolsky's ran1\n");
X  printf("      8 is Press & Teukolsky's ran2\n");
X  printf("      9 is Marsaglia's Ultra\n");
X  printf("    parameter #2 sets the first RNG seed\n");
X  printf("    parameter #3 sets the second RNG seed (if any)\n");
X  printf("    parameter #4 sets the number of times to cycle the RNG\n");
X  printf("      before starting the tests, mod 1 billion.\n");
X  printf("    parameter #5 sets the number of times to cycle the RNG\n");
X  printf("      before starting the tests, div 1 billion.\n");
X  printf("    defaults are 1,1,1,0,0 respectively.\n");
X  printf("  -mnnnn sets the range of the RNG to be nnnn.  Default 100.\n");
X  printf("  -Mnn sets the range of the RNG to be 2^nn, 0<=nn<=31.\n");
X  printf("  -q or -quiet doesn't print any random generates.\n");
X  printf("  -t eliminates most RNG tests, for timing measurements.\n");
X  printf("  -e echos the command line (useful in scripts).\n");
X# ifdef VECTOR
X  printf("  This version of the code is optimized for speed.\n");
X  printf("  Recompile without -DVECTOR if you want any of the following:\n");
X# endif
X  printf("  -f uses (int(dxrandom()*m): faster, but slightly biased.\n");
X  printf("  -p uses random()%%m, a poor method.\n");
X  exit(1);
X} /* end argerr */
X
X#define NORMAL      0
X#define VECTORIZED  1
X#define FAST        2
X#define POOR        3
X
X/* declare storage for a vector of random generates */
X# define VECLENGTH 64
Xstatic long rpt=VECLENGTH;
Xstatic long rvals[VECLENGTH];
X
X# ifndef VECTOR
X/* slower code with runtime options */
Xlong random_value(rng,m,method)
XRNGdata *rng;
Xlong m,method;
X{
X
X  switch (method) {
X    case NORMAL:
X     return( mrandom(m) );
X    case VECTORIZED:
X     if (rpt==VECLENGTH) {
X       rpt = 1;
X       return( mrandomrv(rng,m,VECLENGTH,rvals) );
X     } else {
X       return( rvals[rpt++] );
X     }
X    case FAST:
X     return( (int) (dxrandom()*m) );
X    case POOR:
X     return( lxrandom()%m );
X    default:
X     return( 0 );
X  }
X}
X# else
X/* fast code, an optimized version of the -v option */
X# define random_value(rng,m,method) \
X    (rpt==VECLENGTH ? \
X       ( rpt = 1, mrandomrv(rng,m,VECLENGTH,rvals) ) : \
X       ( rvals[rpt++] ) \
X    )
X# endif /* VECTOR */
X
X/* Comparison routine, for qsort() */
Xint comparelongs(a,b)
Xlong *a,*b;
X{
X  return( ((*a)<(*b))? -1 : (((*a)==(*b))? 0 : 1) );
X}
X
X/* A simple test-driver */
Xint main(argc,argv)
Xint argc; char *argv[];
X{
X  /* command-line arguments, with default values */
X  int reseeding=0; /* nonzero if RNG state will be re-initialized */
X  int quiet=0; /* nonzero if we don't want to print random generates */
X  int echo=0; /* nonzero if we want to echo the mrtest invocation line */
X  int timing=0; /* nonzero if we don't want to make a lot of RNG tests */
X# ifndef VECTOR
X  int method=NORMAL; /* how we should call the mrandom package;
X    * other defined values are VECTORIZED (the -v option),
X    * FAST (-f), and POOR (-p).
X    */
X# else
X  int method = VECTORIZED; /* default (and only) option in this version */
X# endif
X  long alg=1;
X  long mrandom_alg=0;
X  long bufsize=1;
X  long seed1=1,seed2=1, seed[2]; /* seed for RNG */
X  long count1=0,count2=0; /* init call count for RNG */
X  long n=10; /* number of randoms to generate */
X  long m=100; /* desired range of random outputs: 0..m-1 */
X  long discards=0; /* how many generates to discard between outputs */
X
X  long i,j,k; /* temp counters */
X  char rngdesc[RNGIDSTRLEN]; /* for a describe_rng() call */
X
X  RNGdata *myrng; /* area to store this RNG's state & other data */
X
X  double dm,dn; /* = (double) m, (double) n */
X
X  /* temps used to calculate x-squared values */
X  double expf,f,sumfsq; /* we use doubles, not longs, to avoid overflow */
X  long sumf;  /* sum of all freqs, for error check */
X  double zerofreqs; /* number of freqs = 0 */
X  /* flags to indicate which tests will be run */
X  int testequidist; /* 0: max(m,n) too big, don't test;
X             * 1: count freqs online,
X             * 3: m is huge, count freqs offline.
X             */
X  int test2tuples, test3tuples;
X  /* various xsq stats */
X  double xsq,xsq2,xsq3, xsq2l, xsq3l;
X  /* random generates: current, previous, previous^2, first, second */
X  long x, prevx, pprevx, firstx, secondx;
X
X  /* var for max test */
X  long maxx;
X  /* var for median test */
X  long lowcount,median;
X  /* vars for mod 3 test */
X  long mod0, mod1;
X  /* data for equi-distribution test:
X   * freq[i] = sum_{1 \leq j \leq n} (x[j]==i), if testequidist == 1
X   * freq[i] = x[i], if testequidist == 3
X   * freq[i] = undefined, otherwise
X   */
X  long freq[EQUIDISTMAX];
X  /* data for 2-tuple test */
X  long pairs[PAIRMAXM][PAIRMAXM];
X  /* vars for 3-tuple test on low-order 3 bits */
X  long lowpairs[8][8];
X  long lowtrips[8][8][8];
X
X  if(argc > 1) {
X    for(i=1;i<argc;i++) {
X      if(argv[i][0] >= '0' && argv[i][0] <= '9') {
X    n = atol(&(argv[i][0]));
X    if (n < 0) {
X      printf("Illegal value, %ld, for number of random generates.\n",n);
X      argerr((char*) argv[0]);
X        }
X      } else if(argv[i][0] == '-') {
X    switch(argv[i][1]) {
X    case 'S': /* new seed(s) for rng */
X      seed1 = 1; /* defaults */
X      seed2 = 1;
X      count1 = 0;
X      count2 = 0;
X      sscanf (&(argv[i][2]), "%ld,%ld,%ld,%ld,%ld,&ld,&ld",
X         &alg, &seed1, &seed2, &count1, &count2);
X      seed[0]=seed1; seed[1]=seed2;
X      reseeding = 1;
X      if (seed1<0) {
X        printf("Illegal value, %ld, for RNG seed.\n",seed1);
X        argerr((char*) argv[0]);
X          }
X      if (seed2<0) {
X        printf("Illegal value, %ld, for second RNG seed.\n",seed2);
X        argerr((char*) argv[0]);
X          }
X      if (count1<0) {
X        printf("Illegal value, %ld, for number of times ");
X        printf(" to cycle rng before starting tests.\n",count1);
X        argerr((char*) argv[0]);
X          }
X      if (count2<0) {
X        printf("Illegal value, %ld, for number of billions of times ");
X        printf(" to cycle rng before starting tests.\n",count1);
X        argerr((char*) argv[0]);
X          }
X      break;
X    case 'd': /* adjust number of discards */
X      discards = atol(&(argv[i][2]));
X      if (discards < 0) {
X        printf("Illegal value, %ld, for number of discards.\n",discards);
X        argerr((char*) argv[0]);
X          }
X      break;
X    case 'm': /* adjust range of rng */
X      m = atol(&(argv[i][2]));
X      break;
X    case 'M': /* adjust log range of rng */
X      j = atol(&(argv[i][2]));
X      if (j<0 || j>31) {
X        printf("Illegal value, %ld, for log range of rng.\n",j);
X        argerr((char*) argv[0]);
X          }
X      if (j == 31) {
X        m = ~MAXLONG; /* note: m is an unsigned long */
X      }
X      else if (j==0) {
X        m = 0;
X      }
X      else {
X        m = (long)pow( 2.0, (double) j );
X      }
X      break;
X    case 'q': /* quiet! */
X      quiet = 1;
X      break;
X    case 't': /* strip testing code from inner loop, for timing */
X      timing = 1;
X      quiet = 1; /* -t implies -q */
X      break;
X    case 'e':
X      echo = 1; /* echo mrtest invocation */
X      break;
X# ifndef VECTOR
X    case 'v':
X      method = VECTORIZED;
X      break;
X    case 'f':
X      method = FAST;
X      break;
X    case 'p':
X      method = POOR;
X      break;
X# endif VECTOR
X    default:
X      argerr((char*) argv[0]);
X    }
X      } else {
X    argerr((char*) argv[0]);
X      }
X    }
X  }
X
X  if (echo) for (i=0; i<argc; i++) {
X    printf("%s ",argv[i]);
X  }
X  printf("\n");
X
X  dn = (double) n;
X
X  if (m > 0) {
X    dm = (double) m;
X  }
X  else if (m < 0) {
X    dm = (double)(m & 0x7fffffff) + 2147483648.0; /* 2^31 */
X  }
X  else {
X    dm = 4294967296.0; /* 2^32 */
X  }
X
X  if (!reseeding ) {
X    if (access(RNGFILENAME, R_OK)) {
X      printf("There is no RNG statefile in this directory, so ");
X      printf("I'll make one for you.\n");
X      reseeding = 1;
X    }
X  }
X
X  if (reseeding) { /* create a new statefile from scratch */
X    printf("Initializing RNG.\n");
X    myrng=init_rng(alg,mrandom_alg,seed,count1,count2,bufsize);
X  } else { /* use an existing statefile */
X    if (n != 0) printf("Restarting RNG.\n");
X    if (!(myrng=restart_rng(RNGFILENAME))) {
X      exit(1);
X    }
X  }
X  printf(describe_rng(myrng,rngdesc));
X
X  if (n == 0) {
X    if (reseeding) {
X      exit( save_rng(myrng,RNGFILENAME) ); /* save new rng, exit */
X    } else {
X      exit(0); /* immediate exit if n == 0 */
X    }
X  }
X
X  /* set flags: will we run the various tests? */
X  if (timing) { /* avoid time-consuming tests */
X    testequidist = 0;
X    test2tuples = 0;
X    test3tuples = 0;
X  } else {
X    if (m == ~MAXLONG) { /* special case for M==31 (m==2^31) */
X      testequidist = (n <= EQUIDISTMAX ? 3 : 0);
X      test2tuples = 0;
X      test3tuples=1;
X    } else {
X      testequidist = ( m <= EQUIDISTMAX ? 1 : (n <= EQUIDISTMAX ? 3 : 0));
X      test2tuples = ( m <= PAIRMAXM );
X      test3tuples = ( (m%8) == 0 );
X    }
X  }
X
X  /* initialize the various frequency-counting arrays, if needed */
X  if (testequidist==1) for (i=0; i<m; i++) freq[i] = 0;
X  if (test2tuples) for (i=0; i<m; i++) for (j=0; j<m; j++) pairs[i][j] = 0;
X  if (test3tuples) {
X    for (i=0; i<8; i++) {
X      for (j=0; j<8; j++) {
X    lowpairs[i][j] = 0;
X        for (k=0; k<8; k++) {
X      lowtrips[i][j][k] = 0;
X    }
X      }
X    }
X  }
X
X  /* initialize counter, threshhold for median test */
X  lowcount = 0;
X  median = ( (m == ~MAXLONG) ? 1<<30 : m/2 );
X
X  /* initialize counters for mod 3 test */
X  mod0 = 0;
X  mod1 = 0;
X
X  /* initialize storage for max test */
X  maxx = -1;
X
X  if (timing) discards=0; /* -d is meaningless with -t */
X
X  if (!quiet) {
X    printf("Here %s %ld random value", ((n == 1)? "is" : "are"), n);
X  } else {
X    printf("Generating %ld random value", n);
X  }
X  printf("%s", ((n == 1)? " " : "s "));
X  printf("in the range 0 to %ld\n", ((m == ~MAXLONG)? MAXLONG : m-1));
X
X  if (discards > 0) {
X    printf("Note: discards = %ld, so mrandom() will be called a",discards);
X    printf(" total of %.0f times.\n", (1.+(float)discards)*(float)n);
X  }
X
X
X  if (timing) {
X
X    /* vectorized inner loop */
X    for (i=0; i+VECLENGTH<n; i+=VECLENGTH) {
X      mrandomrv(myrng,m,VECLENGTH,rvals);
X      for (j=0; j<VECLENGTH; j++) {
X    /* Note: there are no subroutine calls in this loop! */
X    x = rvals[j];
X    if (x > maxx) maxx = x;
X      }
X    }
X    /* the last (partial) block, if any */
X    for ( ; i<n; i++) {
X      x = mrandomr(myrng,m);
X      if (x > maxx) maxx = x;
X    }
X
X  } else for (i=0; i<n; i++) { /* do lots of RNG tests */
X
X    x = random_value(myrng,m,method);
X    if ( (x >= m && m != ~MAXLONG) || x < 0) {
X      printf("Illegal output, %ld, from mrandom(%ld).\n",x,m);
X      printf("Please contact [e-mail].\n");
X      printf(describe_rng(myrng,rngdesc));
X      exit(1);
X    }
X
X    /* max test */
X    if (x > maxx) maxx = x;
X
X    /* count number of occurrences of each distinct output, if range is small,
X     * otherwise just keep track of outputs (we'll count freqs later)
X     */
X    if (testequidist==1) {
X      freq[x]++;
X    } else if (testequidist==3) {
X      freq[i] = x;
X    }
X
X    /* count pairs */
X    if (i == 0) {
X      firstx = x; /* keep track of first generate for "circular" pairs count */
X    } else {
X      if (test2tuples) {
X    pairs[x][prevx]++;
X    if (i == n-1) pairs[firstx][x]++; /* handle the "wraparound pair" */
X      }
X      if (test3tuples) {
X    lowpairs[x&7][prevx&7]++;
X    if (i == n-1) pairs[firstx&7][x&7]++; /* handle the "wraparound pair" */
X      }
X    }
X
X    /* count triples */
X    if (i == 1) {
X      secondx = x; /* needed for the last "wraparound triplet" */
X    }
X    if (i > 1) {
X      if (test3tuples) {
X    lowtrips[x&7][prevx&7][pprevx&7]++;
X    if (i == n-1) { /* time to do the wraparound triplets? */
X      lowtrips[firstx&7][x&7][prevx&7]++;
X      lowtrips[secondx&7][firstx&7][x&7]++;
X    }
X      }
X    }
X
X    /* count number below median */
X    if (x < median) lowcount++;
X
X    /* count number ==0 and ==1, mod 3 */
X    j = x%3;
X    if (j==0) {
X      mod0++;
X    } else if (j==1) {
X      mod1++;
X    }
X
X    /* keep track of previous two generates, for pair and triple counts */
X    pprevx = prevx;
X    prevx = x;
X
X    /* use fixed-width format, to make it easy to sort the output */
X    if (!quiet) {
X      if (m < 0 && x < 0) { /* print an unsigned long */
X        printf("  %.0f\n", (double)(x&MAXLONG)+(double)MAXLONG+1.0);
X      } else if (m <= 10) {
X        printf("  %d\n", x);
X      } else if (m <= 100) {
X        printf("  %2d\n", x);
X      } else if (m <= 1000) {
X        printf("  %3d\n", x);
X      } else {
X        printf("  %4d\n", x);
X      }
X    }
X    /* now discard some generates: d=2^k, k>9, gives nrand48() trouble */
X    for (j=0; j<discards; j++) {
X      /* note: we don't even check whether x is in range */
X      x = random_value(myrng,m,method);
X    }
X  }
X
X  printf("Max value from %ld call%s", n, (n==1? "" : "s"));
X  if (discards) printf(" (discards = %ld)",discards);
X  printf(" to mrandom(%ld) = %ld\n",m,maxx);
X  if ( (maxx >= m && m != ~MAXLONG) || maxx < 0) {
X    printf("Illegal value for max!\n");
X    printf("Please contact [e-mail].\n");
X    printf(describe_rng(myrng,rngdesc));
X    exit(1);
X  }
X
X  if (!timing && n>2) {
X    printf("Equi-distribution test results:\n");
X    if (testequidist == 0) {
X      printf("  Range m and/or number of random generates n is too large.\n");
X      printf("  Max(m,n) for this test is %d.\n", EQUIDISTMAX);
X    } else if (m == 1) {
X      printf("  You need m > 1 for this test.\n");
X    } else {
X      sumfsq = 0.0;
X      expf = dn/dm; /* expected freq[] */
X      if (testequidist == 1) { /* freq[] array has frequencies */
X        sumf = 0; /* error check: is ((sumf = \sum_i{freq[i]}) == n)? */
X        for (i=0; i<m; i++) {
X      sumf += freq[i];
X          f = ((double)freq[i]) - expf;
X          sumfsq += f*f;
X        }
X        if (sumf != n) {
X      printf("Warning: a frequency counter has overflowed.");
X      printf("  Test invalid!\n");
X      goto testinvalid;
X        }
X      } else if (testequidist == 3) {
X    /* freq[] array has the random generates: sort, count dups */
X        qsort(freq,n,sizeof(long),comparelongs);
X    f = 1.0;
X    zerofreqs = dm; /* number of integers NOT seen in RNG output */
X    for (i=0; i<n; i++) {
X      if (i == (n-1) || freq[i] != freq[i+1]) {
X        sumfsq += (f-expf)*(f-expf); /* end of one run of duplicates... */
X        zerofreqs -= 1.0;
X        f = 1.0; /* ...and the beginning of another run */
X      } else {
X        f += 1.0; /* a "match" extends a run */
X      }
X    }
X    /* add contributions from zero-frequency outputs */
X    sumfsq += zerofreqs*(-expf)*(-expf);
X      }
X      /* Note: many texts suggest summing t += freq[i]^2, then calculating
X       * xsq by (dm*t/dn)-dn.  This is slightly faster, but MUCH less
X       * accurate when dn is much larger (e.g. 2^20 times larger) than dm.
X       */
X      xsq = sumfsq/expf;
X      interpret_xsq(xsq, dm, dn);
X      /* End of equi-distribution test code */
X
X      printf("Pairwise correlation test results:\n");
X      if (!test2tuples) {
X        printf("  Range is too large.  Max m for this test is %d.\n", PAIRMAXM);
X      } else {
X        sumfsq = 0.0;
X        expf = dn/(dm*dm); /* expected freq[] */
X        for (i=0; i<m; i++) for (j=0; j<m; j++) {
X      f = (double)pairs[i][j] - expf;
X          sumfsq += f*f;
X        }
X        xsq2 = sumfsq/expf;
X    /* note that we use the xsq from the equi-distribution test */
X        interpret_xsq(xsq2-xsq, dm*dm-dm+1.0, dn);
X      } /* end of pairwise correlation test code */
X    }
X
X    printf("3-tuple, low-order 3 bits, correlation test results:\n");
X    if (!test3tuples) {
X      printf("  You need (m mod 8) == 0 for this test.\n");
X    } else {
X      sumfsq = 0.0;
X      expf = dn/64.0; /* expected freq in each bin */
X      for (i=0; i<8; i++) {
X        for (j=0; j<8; j++) {
X          f = (double)lowpairs[i][j] - expf;
X          sumfsq += f*f;
X        }
X      }
X      xsq2l = sumfsq/expf;
X      sumfsq = 0.0;
X      expf = dn/512.0; /* expected freq in each bin */
X      for (i=0; i<8; i++) {
X        for (j=0; j<8; j++) {
X      for (k=0; k<8; k++) {
X            f = (double)lowtrips[i][j][k] - expf;
X            sumfsq += f*f;
X      }
X        }
X      }
X      xsq3l = sumfsq/expf;
X      interpret_xsq(xsq3l-xsq2l, 8.0*8.0*8.0-8.0*8.0+1.0, dn);
X    } /* end of 3-tuple test code */
X
X    printf("Most-significant bit test results:\n");
X    f = (double)lowcount; /* freq below m/2 */
X    expf = dn*(floor(dm/2.))/dm; /* expected val of f */
X    sumfsq = ((f-expf)*(f-expf))/expf; /* normalized dev^2 from expectation */
X    f = dn-f; /* freq at or above m/2 */
X    expf = dn-expf;
X    sumfsq += (f-expf)*(f-expf)/expf;
X    interpret_xsq(sumfsq,2.0,dn);
X    /* end of MSB test code */
X
X    printf("Mod-3 test results:\n");
X    /* We need the freq[] to be reasonably equal for xsq analysis */
X    if (!(m == 3 || m > 5 || m == ~MAXLONG)) {
X      printf("  You need m==3 or m>5 to run this test.\n");
X    } else {
X      f = (double)mod0; /* freq of x==0 mod 3 */
X      expf = dn*(ceil(dm/3.))/dm;
X      sumfsq = ((f-expf)*(f-expf))/expf;
X      f = (double)mod1; /* freq of x==1 mod 3 */
X      expf = dn*(ceil(2.*dm/3.))/dm - expf;
X      sumfsq += (f-expf)*(f-expf)/expf;
X      f = dn-f-(double)mod0; /* freq of x==2 mod 3 */
X      expf = dn-dn*(ceil(2.*dm/3.))/dm;
X      sumfsq += (f-expf)*(f-expf)/expf;
X      interpret_xsq(sumfsq,3.0,dn);
X    } /* end of mod-3 test code */
X  } else if (n < 3) {
X    printf("At least 3 random generates are required for the testing module.\n"
X      );
X  }
Xtestinvalid:
X  /* terminate job */
X  printf("Final ");
X  printf(describe_rng(myrng,rngdesc));
X
X  return( save_rng(myrng,RNGFILENAME) );
X
X} /* end main */
END_OF_FILE
if test 21361 -ne `wc -c <'src/mrtest.c'`; then
    echo shar: \"'src/mrtest.c'\" unpacked with wrong size!
fi
# end of 'src/mrtest.c'
fi
echo shar: End of archive 2 \(of 6\).
cp /dev/null ark2isdone
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