[comp.lang.c] random number generator

PEPRBV%CFAAMP.BITNET@husc6.harvard.EDU (Bob Babcock) (02/14/88)

>I need a random number generator with a Gaussian distribution. It
>should have a settable standard deviation. I would prefer it
>written in C but any language would be appreciated. Thanks

Check out Knuth's Seminumerical Algorithms book, section 3.4.1.
Several methods are given, in a form that could easily be
converted to code.

mcdonald@uxe.cso.uiuc.edu (02/16/88)

>Simplistic, but this works:
>Generate about 12 random numbers (assuming 0.0 <= i <1.0)

       ADD THEM TOGETHER

>and divide by 12.  This is sufficiently close to gaussian distribution
>for most purposes.  The more numbers you generate before you do the
>division, the closer you will get to purely gaussian distribution.

joseph@chromo.ucsc.edu (Joseph Reger) (02/18/88)

In article <225800004@uxe.cso.uiuc.edu> mcdonald@uxe.cso.uiuc.edu writes:
>
>>Simplistic, but this works:
>>Generate about 12 random numbers (assuming 0.0 <= i <1.0)
>
>       ADD THEM TOGETHER
>
>>and divide by 12.  This is sufficiently close to gaussian distribution
>>for most purposes.  The more numbers you generate before you do the
>>division, the closer you will get to purely gaussian distribution.


Devide by 12??

This method uses the central limit theorem and gives nearly
gaussian (normal) random numbers if the number of uniform
random numbers to be summed up is sufficiently large:
   Take N uniform random deviates on (0,1): u(1), u(2),...,u(N).
Then g = sqrt(12/N) * ( (u(1)+u(2)+...+u(N)) - N/2 )
yields a (nearly) normal variable with mean 0 and variance 1.
N=12 seems to be sufficiently large for most purposes and
avoids computing the sqrt() factor.
(So do not devide by 12, subtract 6!). 
g' = M + D * g gives a normal random number with mean M and
standard deviation D.

karl@haddock.ISC.COM (Karl Heuer) (02/19/88)

Followups to sci.math; this has nothing to do with C anymore.

>[A good approximation to a Gaussian distribution can be obtained by
>generating about 12 independent random numbers uniformly from [0,1), adding
>them together, and dividing by 12.]

Actually, the reason for using 12 of them (not "about 12") is that you don't
need to rescale.  The variance of the uniform distribution is 1/12, so adding
12 of them gives you a distribution with variance 1.  If you subtract 6, you
then have a good approximation of a standard Gaussian.

I happen to prefer the polar-coordinate version, because it's mathematically
exact, requires only one call to the uniform RNG per result, and is faster (or
so I've been told).

Karl W. Z. Heuer (ima!haddock!karl or karl@haddock.isc.com), The Walking Lint

gwc@root.co.uk (Geoff Clare) (03/01/88)

In article <225800004@uxe.cso.uiuc.edu> mcdonald@uxe.cso.uiuc.edu writes:
>
>>Simplistic, but this works:
>>Generate about 12 random numbers (assuming 0.0 <= i <1.0)
>
>       ADD THEM TOGETHER

	SUBTRACT 6   !!!!

>
>>and divide by 12.  This is sufficiently close to gaussian distribution
>>for most purposes.  The more numbers you generate before you do the
>>division, the closer you will get to purely gaussian distribution.

NO NO NO NO NO!

Only *exactly* 12 random numbers used in this way will give an
approximate Gaussian.  Any other number will not give the right variance.
(The variance of a uniform distribution over [0,1] is 1/12).

-- 

Geoff Clare    UniSoft Limited, Saunderson House, Hayne Street, London EC1A 9HH
gwc@root.co.uk   ...!mcvax!ukc!root44!gwc   +44-1-606-7799  FAX: +44-1-726-2750

ccel@chance.uucp (CCEL) (12/21/89)

Greetings 

I was trying to generate a random number bewteen 0 and maxlongint (i.e. 2^32)
on my machine (80386 with Interactive). The rand() function that I have
returns a psuedo-random number between 0 and MAXINT (which in this case
is 32768). I thought about generating two random numbers and multiplying
them together, but I realized that would produce a bell-curve-shaped
distrubution (the numbers in the middle would be more likely to occur).

Does anyone have any ideas or any code that generates random numbers,
or anyplace that I can check to find such code?

               Randy Tidd
rtidd@mwsun.mitre.org
Randy Tidd                   MITRE-McLean CCEL Lab
rtidd@mitre.arpa             ccel%community-chest@gateway.mitre.org
#define DISCLAIM TRUE

jonathan.forbes@canremote.uucp (JONATHAN FORBES) (12/22/89)

It depends how random you want your numbers to be.  If merely using the 
timer for randomising isn't enough, you can just add together the 
addresses of several of your variables and seed the generator with that.
 Poking around in memory and adding up random memory locations should 
produce a decent random number generator (although a very hacked up and 
non-portable one.)
---
 * Via ProDoor 3.1R 

gary@hpavla.HP.COM (Gary Jackoway) (12/22/89)

> I was trying to generate a random number bewteen 0 and maxlongint (i.e. 2^32)
> on my machine (80386 with Interactive). The rand() function that I have
> returns a psuedo-random number between 0 and MAXINT (which in this case
> is 32768). 
----------
Simply put two rand's together:

       longrand = (unsigned long) rand() | ((unsigned long) rand() << 32);

This should give you a uniform distribution.
I should warn you that rand() is often NOT a good random number generator,
so if the application is sensitive, be sure to run randomness tests on
your long random numbers.

No charge,

Gary Jackoway

condict@cs.vu.nl (Michael Condict) (12/22/89)

In article <83943@linus.UUCP> rtidd@mwsun.mitre.org writes:
|Greetings 
|
|I was trying to generate a random number bewteen 0 and maxlongint (i.e. 2^32)
|on my machine (80386 with Interactive). The rand() function that I have
|returns a psuedo-random number between 0 and MAXINT (which in this case
|is 32768). I thought about generating two random numbers and multiplying
|them together, but I realized that would produce a bell-curve-shaped
|distrubution (the numbers in the middle would be more likely to occur).
|
|Does anyone have any ideas or any code that generates random numbers,
|or anyplace that I can check to find such code?
|

Generate the two 16-bit numbers and concatenate them:

	long_random = (short_random_1 << 16) + short_random_2;

Or, if you prefer:

	long_random = short_random_1 * 65536 + short_random_2;

If the rand function generates uniformly distributed random shorts,
this will obviously produce uniformly distributed random longs.
-- 
Michael Condict		condict@cs.vu.nl
Vrije University
Amsterdam

rta@pixar.UUCP (Rick Ace) (12/23/89)

In article <83943@linus.UUCP>, ccel@chance.uucp (CCEL) writes:
> Does anyone have any ideas or any code that generates random numbers,
> or anyplace that I can check to find such code?
> 
>                Randy Tidd

Try Knuth, Seminumerical Algorithms.

Also, see June 1988 CACM for an article by Pierre L'Ecuyer, titled
   "Efficient and Portable Combined Random Number Generators"
There's theory and some practical examples you can put to use,
along with references at the end to other literature.

Rick Ace
Pixar
3240 Kerner Blvd, San Rafael CA 94901
...!{sun,ucbvax}!pixar!rta

bill@twwells.com (T. William Wells) (12/25/89)

In article <4943@condict.cs.vu.nl> condict@cs.vu.nl (Michael Condict) writes:
: Generate the two 16-bit numbers and concatenate them:
:
:       long_random = (short_random_1 << 16) + short_random_2;
:
: Or, if you prefer:
:
:       long_random = short_random_1 * 65536 + short_random_2;

In both cases, safety requires a cast:

	long_random = ((unsigned long)short_random_1 << 16) + short_random_2;
	long_random = (unsigned long)short_random_1 * 65536 + short_random_2;

Else the first time you compile on a machine with 16 bit ints,
you'll get a surprise.

---
Bill                    { uunet | novavax | ankh | sunvice } !twwells!bill
bill@twwells.com

bill@twwells.com (T. William Wells) (12/25/89)

In article <940001@hpavla.HP.COM> gary@hpavla.HP.COM (Gary Jackoway) writes:
:        longrand = (unsigned long) rand() | ((unsigned long) rand() << 32);
									^16?

Note that normally this code would be incorrect C code, in that
it twice evaluates a function that might produce different
effects on each invocation in such a way that the order of
evaluation is undefined.

But because rand() ought to be random, it should not matter which
order you get the results in. So, in this rare case, you can get
away with it. There is a cost, however: since the order of
evaluation is undefined, the compiler might choose to do it either
way in any given compilation, which could make reproducing the
results of your program difficult. (Say, for debugging.)

---
Bill                    { uunet | novavax | ankh | sunvice } !twwells!bill
bill@twwells.com

harish@csl36h.ncsu.edu (Harish Hiriyannaiah) (12/27/89)

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, if

	z = a1*X + a2*Y , 

then

       pdf(z) = (a1*pdf(X)) *** (a2*pdf(Y)).

Here *** denotes convolution, X and Y are independent random variables
and pdf() is the probability
density function.

If pdf(X) and pdf(Y) are uniform, then pdf(z) is a trapezoid,
degenerating to a triangular
pdf if pdf(X) = pdf(Y), with a1 = a2.

Now, if the convolution is between a very wide uniform distribution and
a very narrow
uniform distribution, then the resulting trapeziodal pdf MAY be
approximated to a uniform
distribution. So netters should be aware of this approximation and not
trust the
values of z near the ends of the distribution. i.e
   _
  | |
  | |
  | |
  | |                                               
_____________________________________
  | |    ***    __________________________    =    /                    
                  \
  | |          |                          |       /                     
                   \
   w1                      w2                     ^^^^             w3   
                 ^^^^
                                                  Areas of
                                                  approximation not
                                                  to be trusted.


In our case, if pdf(rand1) = pdf(rand2) ( which is the case), and if a1
>> a2, (which is also
true), then the above scenario results. Since I am lazy to delve into a
book on probability
to analytically determine the triangular portions of the trapezoidal
distribution in terms
of a1 and a2, I suggest that the netters carefully look into it before
using the r.v's like
rand. (Actually, if w1 is the width of the narrow unif pdf, and w2 that
of the wide pdf, then
w3 = w1 + w2, and each triagular region is of width w1, with the flat
region having a height
of 1/(w1*w2). )
      
harish pu. hi.						harish@ecebucolix.ncsu.edu
							harish@ecelet.ncsu.edu

harish@csl36h.ncsu.edu (Harish Hiriyannaiah) (12/27/89)

Damn XRN. It messed up the previous posting.

harish pu. hi			harish@ecelet.ncsu.edu

wilson@uhccux.uhcc.hawaii.edu (Tom Wilson) (12/27/89)

In article <1989Dec26.164103.17945@ncsuvx.ncsu.edu> harish@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, if
>
>	z = a1*X + a2*Y , 
>
>then
>
>       pdf(z) = (a1*pdf(X)) *** (a2*pdf(Y)).
>
>Here *** denotes convolution, X and Y are independent random variables
>and pdf() is the probability
>density function.
>
>If pdf(X) and pdf(Y) are uniform, then pdf(z) is a trapezoid,
>degenerating to a triangular
>pdf if pdf(X) = pdf(Y), with a1 = a2.

>harish pu. hi.					harish@ecebucolix.ncsu.edu
>						harish@ecelet.ncsu.edu

I am wiling to accept your analysis for two continuous pdf's.  However, in the
case being considered here, rand1 and rand2 are uniform discrete distributions
on [0 .. 65535] (note 65535 = 2^16-1, a bit pattern of all 1s in a 16-bit
integer).  Then rand1*65536 + rand2 is the same as shifting rand1 left by
16 bits and or'ing rand2 into the lower 16 bits.  For each integer in 
[0 .. 65536^2 -1 ], there is a unique way that it can arise from this sum.
So it seems to me that P(k) = (1/65536)^2 for k in [0..65536^2-1],
because the two underlying terms each occur with P = 1/65536; therefore, it
is a uniform random distribution on the 32-bit integers.  

In general, 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.

Another thought:  this can be viewed not as a sum, but as picking
a coordinate in [0..65535] x [0..65535] using the two independent uniform
random variables.  The mapping from this rectangular region to [0..65536^2-1]
is 1-1 and onto, so the resulting distribution is uniform.


-- 
Tom Wilson                        wilson@uhccux.uhcc.Hawaii.Edu (Internet)
                                  wilson@uhccux.UUCP || wilson@uhccux (Bitnet)

john@frog.UUCP (John Woods) (12/28/89)

In article <940001@hpavla.HP.COM>, gary@hpavla.HP.COM (Gary Jackoway) writes:
> > I was trying to generate a random number bewteen 0 and maxlongint (i.e. 2^32)
> > on my machine (80386 with Interactive). The rand() function that I have
> > returns a psuedo-random number between 0 and MAXINT (which in this case
> > is 32768). 

I assume you mean 32767.

>        longrand = (unsigned long) rand() | ((unsigned long) rand() << 32);

This is the second wrong answer I've seen.  A correct answer is

	typedef unsigned long ulong;

	longrand = (ulong)rand() | ((ulong)rand() << 15)
				 | (((ulong)rand()&3) << 30);

The first wrong answer shifted by 16, but rand() as described returns a 15-bit
integer.

A better solution (IMHO) is to grab Knuth Vol. 2 and code up a good 32 bit
RNG of your own.
-- 
John Woods, Charles River Data Systems, Framingham MA, (508) 626-1101
...!decvax!frog!john, john@frog.UUCP, ...!mit-eddie!jfw, jfw@eddie.mit.edu

Happiness is Planet Earth in your rear-view mirror.	- Sam Hurt

chris@mimsy.umd.edu (Chris Torek) (12/28/89)

>In article <1989Dec26.164103.17945@ncsuvx.ncsu.edu>
>harish@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@uhccux.uhcc.hawaii.edu> wilson@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@cs.umd.edu	Path:	uunet!mimsy!chris

staceyc@sco.COM (Stacey Campbell) (12/29/89)

In article <940001@hpavla.HP.COM>, gary@hpavla.HP.COM (Gary Jackoway) writes:
>I was trying to generate a random number bewteen 0 and maxlongint (i.e. 2^32)

Stop using rand().  The Unix version you are using has an excellent
pseudo-random number generator in the C library.  Check the man page
for drand48().  Long random numbers are returned by lrand48().
-- 
Stacey Campbell                                             _--_|\
{uunet,ucscc,decwrl,att,microsoft,wyse}!sco!staceyc        /      \
staceyc@sco.com                                            \_.--._/
                                                                 v

tswenson@dgis.dtic.dla.mil (Timothy Swenson) (07/24/90)

	Does anyone have a random number generator
written in C?  Or better yet, Small-C.  I really need
one for a project I'm working on, and my compiler does
not have one in it's library.  I'll take anything that 
I can get (even helpfull hints).
	Thanks in advance,

	Tim Swenson
	tswenson@dgis.dtic.dla.mil

wozniak@utkux1.utk.edu (Bryon Lape) (07/24/90)

In article <941@dgis.dtic.dla.mil> tswenson@dgis.dtic.dla.mil (Timothy Swenson) writes:
>
>	Does anyone have a random number generator
>written in C?  Or better yet, Small-C.  I really need
>one for a project I'm working on, and my compiler does
>not have one in it's library.  I'll take anything that 

	I thought that rand() and srand() are standard library
functions!  These will generate random numbers.  srand() is the seed
function while rand() returns a pseudo-random number.


-bryon lape-

paul@u02.svl.cdc.com (paul kohlmiller) (07/25/90)

wozniak@utkux1.utk.edu (Bryon Lape) writes:

>In article <941@dgis.dtic.dla.mil> tswenson@dgis.dtic.dla.mil (Timothy Swenson) writes:
>>
>>	Does anyone have a random number generator
>>written in C?  Or better yet, Small-C.  I really need
>>one for a project I'm working on, and my compiler does
>>not have one in it's library.  I'll take anything that 

>	I thought that rand() and srand() are standard library
>functions!  These will generate random numbers.  srand() is the seed
>function while rand() returns a pseudo-random number.

Indeed the standard does define rand and srand and even puts the code in the
standard itself.
static unsigned long int next = 1;
int rand(void) /* RAND-MAX assumed to be 32767 */
{
next=next*1103515245+12345;
return (unsigned int) (next/65536) % 32768;
}
The srand function simply resets the value of next.
This code is an example only and it is not a requirement that an ANSI 
conforming implementation use this version of rand.
Paul Kohlmiller
standard disclaimers