Message-Id:	<200004251246.IAA11271@delorie.com>
From:	"Dieter Buerssner" <buers AT gmx DOT de>
To:	djgpp-workers AT delorie DOT com
Date:	Tue, 25 Apr 2000 15:50:51 +0200
MIME-Version:	1.0
Subject:	rand() in libc
X-mailer:	Pegasus Mail for Win32 (v3.12b)
Reply-To:	djgpp-workers AT delorie DOT com

I have seen, that the rand() in libc improved very much in djgpp 2.03 
(or perhaps even earlier). It seems clear, that nobody will expect a 
high qualtity PRNG in the Standard C library. I nevertheless want to 
suggest one little improvement. The current code looks like this:

#include <stdlib.h>

static unsigned long long next = 1;

int
rand(void)
{
  /* This multiplier was obtained from Knuth, D.E., "The Art of
     Computer Programming," Vol 2, Seminumerical Algorithms, Third
     Edition, Addison-Wesley, 1998, p. 106 (line 26) & p. 108 */
  next = next * 6364136223846793005LL + 1;
  /* was: next = next * 0x5deece66dLL + 11; */
  return (int)((next >> 21) & RAND_MAX);
}

Substituting the return line with

  return (int)((next >> 32) & RAND_MAX);

will, If I understand correctly, make the period 2^11 times longer 
(2^63 vs. 2^52). With this, rand() will pass all diehard tests (a 
random number test suite by George Marsaglia), as well as many 
reparameterized diehard tests and some tests of my own. The original 
version will fail some tests. (Some of the tests probably will fail 
soon, when a bit more stressed, even with the larger shift.)
I think, no patch will be needed for this.

A preferrable alternative, that runs at about double speed here
(needs one MUL, one ADD, one ADC and a few MOVs, compared
to three MULs needed by the above code), that seems to pass all my 
tests very well is the following multiply with carry RNG. (Also 
suggested by Marsaglia)

#include <stdlib.h>

/* Multiply with carry RNG suggested by George Marsaglia.
  This just does a 32x32 -> 64 bit multiplication of the last
  random number with mul and adds the carry of the
  last multiplication to it. The new carry will be the high word,     
  the new random number the low word. This implementation is          
  somewhat microoptimized for Gcc. The period is ca. 2^62 */

static unsigned long long zseed = (unsigned long long)12345<<32;

int 
rand(void)
{
  static unsigned long mul=2051013963UL;
  unsigned long l1, l2;
  unsigned long long res;

  l1 = (unsigned long)(zseed&0xffffffffUL);
  l2 = zseed>>32;
  res = l2 + l1*(unsigned long long)mul;
  zseed = res;
  return (int)(res & RAND_MAX);
}

void 
srand(unsigned int s)
{
  zseed = (((unsigned long long)12345)<<32) | s;
}

When you want to consider this, please don't make the type of mul 
static const unsigned long.

I also want to comment the libc info of rand, where I read:

To get a random number in the range 0..N, use `rand()%(N+1)'.  Note
that the low bits of the `rand''s return value are not very random, so
`rand()%N' for small values of N could be not enough random.  The
alternative, but non-ANSI, function `random' is better if N is small.
*Note random::.

I think this is somewhat wrong advice. Eli told me that the libc info 
is not the place, to teach people about random numbers. This is 
certainly correct. However no advice might be better than wrong 
advice. 

1) I am not too certain, that random will be better than rand. It 
will fail some tests in a spectacular fashion. I think random should 
be depreciated for any other thing, than porting old BSD code. 

2) From the above paragraph, one would conclude, that rand()%N works 
well for large N. This is not the case, when N is not a power of two,
because different values in the range will have significantly 
different probabilities.

An alternative suggestion (for N much smaller RAND_MAX) could be

  (int)((rand()*(N+1.0))/(RAND_MAX+1.0))

or a similar Gcc specific formulation using unsigned long long 
instead of double.

Regards,

Dieter "today with signature" Buerssner

-- 
       A random number generator is like sex;
       When it's good, it's wonderful,
       And when it's bad, it's still pretty good.

       And when it's bad, try a twosome or threesome. (George Marsaglia)

- Raw text -