## uniform floating point random number generator - summary

### uniform floating point random number generator - summary

A few days ago I requested help on generating uniform random numbers
in the range (0..1] on a BSD 4.2 system.  Here is a summary of the help
I received and the method I finally chose.

Note:  the information described here is also directly applicable to
BSD 4.3 .  This method could also be adapted for use on SYSV, but I
a specific implementation.

Note:  in the description below, I will not quote any contributors by
name, since I have neither asked for, nor received their permission to
do so.

To begin with, it was pointed out to me that I had misstated the desired
range of the random number generator.  It should be [0..1) .  Oops :-).

Several people sent me references for random number generation.

1. P.  L'Ecuyer, "Efficient  and  portable  combined  random  number
generators", Communications of the ACM 31 (1988) 742-749.

2. P. Bratley, B.L. Fox and L.E. Schrage, "A  guide to simulation",
(Springer-Verlag, New York, NY, 1987).

3. "Efficient and Portable Combined Random Number Generators",
Pierre L'Ecuyer, Communications of the ACM, Vol. 31, No. 6 (June 1989).

I suspect the 1st and 3rd references are the same; it's probable that
one of the cited dates is incorrect.

I can also add the reference for the algorithm that I mentioned
when I posted my query (which overflows when implemented with unsigned
arithmetic in C):

4. J. Banks and J. S. Carson, "Discrete-Event System Simulation",
Prentice-Hall, 1984, ISBN 0-13-215582-6, Chapter 7, Example 7.12,
pages 266-267.  Chapter 7 of this book describes random number
generation in detail and the methods used to test generators for
uniformity and independence.

Continuing, one individual told me to check out Knuth's books and
one individual advised me to consult the Numerical Recipes In
C/Pascal/Fortran books.  Another individual sent me his implementation
of an algorithm cited in reference #1 above.  If you'd like a copy of
his implementation, drop me a note.

Finally, two individuals told me to use the generator provided by the
"random()" library function.  At first, I discounted these suggestions
because random() returns integers.  But I carefully re-read reference
4 above (specifically example 7.12) and immediately noticed that
if random() returns uniform integers in a _known_ range, then it
can be used to generate uniform floating point values in the range
[0..1) with one arithmetic operation (division).  So here's the
method I chose to use:

/*
* Return a "uniformly" distributed double precision random
* number in the range [0..1).
*
* CAVEAT:  only works on a 32-bit architecture.
*/
double
drandom()
{
return (((double) random()) / ((unsigned) 0x80000000));

Quote:}

The key to note here is that random() returns integers in the range
0 to 2^31-1 on MY MACHINE (your mileage WILL vary).  To obtain the
desired floating point values, it was only necessary to divide
random()'s result by 2^31 (0x80000000).

But does it work (is the result uniform)?  That, of course, depends
upon whether or not random() is uniform.  I asked a statistician how
to test the results.  She said, "There are many methods.  Three of
the simplest are:

- generate a large sample size, bin the results, and look at the
distribution.  If the bins look skewed, you probably have problems.
If the bins look uniformly distributed, verify the results with one
of the following two tests.

- Chi-squared goodness of fit test.

- Kolmogrov-Smirnov test."

So, I generated 10,000 samples and binned the results:

Bin  Total
0: 504
1: 504
2: 514
3: 495
4: 515
5: 489
6: 478
7: 511
8: 469
9: 521
10: 493
11: 525
12: 498
13: 477
14: 532
15: 479
16: 473
17: 520
18: 491
19: 512

As far as I was concerned, this was plenty uniform for my purposes.
I did not pursue any other tests.

[ I'm sorry I couldn't thread this to my original posting, but it
expired on our news server. ]

I would like to thank everyone who replied.  The nearly unanimous

gave the best summary:

The xdr library seems to be available on every machine that supports
nfs and rpc -- just about everybody these days.  xdr uses the IEEE-754
standard for floating point, so the conversion from double is a no-op
for most machines.

Given the near universal usage of IEEE-754 representation, a fact I
was previously unaware of, I have decided to use it directly in my
protocol.  To keep things simple, I have decided not to use the xdr
library until such time as a conversion to IEEE-754 is actually
required.

Peter Johansson

6. Test