Random long-number generator

[Chris Torek]

In article <6482 at uhccux.uhcc.hawaii.edu> caprio at uhccux.uhcc.hawaii.edu
(Mike Caprio) quotes an algorithm from
>Dudewicz, E.J., Z. A. Karian and R. J. Marshall III.  1985.
>Random number generation on microcomputers. pp. 9-14. In R. G.
>Lavery [ed.], Modeling and simulation on microcomputers:1985.
>Society for Computer Simulation, La Jolla.

which goes as follows:

>long  seed;
> 
>long lrand ()
>{
> 
>  seed = seed * 1220703125L;
>  if (seed < 0)
>   seed += 2147483648L;
>}

This code is not portable, because it depends on integer overflow being
ignored.  (It also fails to return a value from lrand---no doubt simply
an oversight.)  It can be made portable quite simply: replace the
`long's with `unsigned long's, add a mask (in case ULONG_MAX >
0xffffffff, e.g., if longs are 64 or 256 bits or something [yes, 256
bit longs are legal, if unlikely]), and use a cast or two:

	unsigned long seed;

	long lrand()
	{

		seed *= 1220703125;	/* the L suffix is unnecessary */
		seed &= 0xffffffff;	/* mask to 32 bits */
		if ((long)seed < 0)
			seed += 0x80000000;
		return ((long)seed);
	}

Next, the test for ((long)seed < 0) can be eliminated by observing the
following:

	1. (long)seed will be < 0 in both 1s and 2s complement arithmetic
	   only when (a) long is 32 bits and (b) bit 31 is set.
	2. In both cases, adding 0x80000000 will (since long is 32 bits
	   and unsigned arithmetic is modular) simply clear bit 31.

Hence, to obtain the same behaviour as before, we can simply mask bit 31
while masking other possible bits:

	long lrand()
	{
		seed *= 1220703125;
		seed &= 0x7fffffff;
		return ((long)seed);
	}

Finally, although I am no numerical analyst nor statistician and get my
random numbers from Knuth vol. 2 (or more simply, from the 4BSD `random()'
library), one observation:

>[lrand()] produces integers in the range 0-2^31 ....

As written, lrand() must never produce the value 0, because if it does
it will continue to produce 0 indefinitely.  Most multiplicative random
number generators involve an addition step as well, so that the algorith
is

	seed = (A * seed) + B;

for two (fixed) nonzero values of A and B.  A zero value for either one
results in a generator with at least one fixed point (when seed is 0).
-- 
In-Real-Life: Chris Torek, Univ of MD Comp Sci Dept (+1 301 454 7163)
Domain:	chris at cs.umd.edu	Path:	uunet!mimsy!chris