Random number generator

[Chris Torek]

>In article <1989Dec26.164103.17945 at ncsuvx.ncsu.edu>
>harish at ecebucolix.ncsu.edu writes:
>>A lot of netters have suggested using  "rand = rand1*655536 + rand2" to
>>create a uniform random number, where rand1 and rand2 are two uniform
>>r.v's.  It should be noted that rand will NOT be uniform, since ...
[analysis deleted]

In article <5825 at uhccux.uhcc.hawaii.edu> wilson at uhccux.uhcc.hawaii.edu
(Tom Wilson) writes:
[it is uniform because]
>... we are looking at the very special case of two discrete 
>uniform distributions on [0..n-1], and forming the new distribution
>rand = n*rand1 + rand2.

Wilson is right for the reasons he gives.  Unfortunately, both analyses
are ultimately wrong, because the `random number generators' used in
computer programming are not truly random: although the two random
variables are discrete and (probably) uniformly distributed, they are
likely to be codependent.  In particular, look at the common (bad)
random number generator found on many Unix systems:

	static long randx = 1;

	void
	srand(x)
		unsigned x;
	{

		randx = x;
	}

	rand()
	{

		return ((randx = randx * 1103515245 + 12345) & RAND_MAX);
	}

This may return a uniform distribution---I do not know that it does---but
it is not very random at all.  In fact, the low $n$ bits repeat with a
period of $n$, at least for n in 0..7.  (See program below.)  Hence
we know, for instance, that (if this rand() is used) the two values
returned from two successive calls will be alternately odd and even.
This will skew the distribution of the combined numbers.

In short, all this reasoning about random distributions is a nice example
of how to apply mathematics to computers to get a brilliantly wrong proof.
:-)  The premise is wrong: rand() does not return random values.

Appendix A: a program to see whether the low bits cycle in the first 512
calls.  On most Unix systems, it produces no output.

#define n_ones(n) ((1 << (n)) - 1)

main()
{
	register int i, t;
	int r[256];

	for (i = 0; i < 256; i++)
		r[i] = rand();
	for (i = 0; i < 256; i++) {
		t = rand();
#define test(k) \
	if ((t & n_ones(k)) != (r[i] & n_ones(k))) \
		printf("low %d bits do not match (%2.2x vs %2.2x)\n", k, \
			t & n_ones(k), r[i] & n_ones(k));
		test(1);
		test(2);
		test(3);
		test(4);
		test(5);
		test(6);
		test(7);
		test(8);
	}
	exit(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