From: buers AT gmx DOT de (Dieter Buerssner)
Newsgroups: comp.os.msdos.djgpp
Subject: Re: Random numbers
Date: 30 May 2000 08:48:59 GMT
Lines: 101
Message-ID: <8h06hp.3vvrb0l.0@buerssner-17104.user.cis.dfn.de>
References: <392CB951 DOT 34B432E6 AT chemistry DOT uq DOT edu DOT au> <Pine DOT SUN DOT 3 DOT 91 DOT 1000525084705 DOT 21360M-100000 AT is> <cs95jscumbafalh5jj0f8l6dm1u0uiq5o5 AT 4ax DOT com>
NNTP-Posting-Host: pec-45-170.tnt3.s2.uunet.de (149.225.45.170)
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 8bit
X-Trace: fu-berlin.de 959676539 2292376 149.225.45.170 (16 [17104])
X-Posting-Agent: Hamster/1.3.13.0
User-Agent: Xnews/03.02.04
To: djgpp AT delorie DOT com
DJ-Gateway: from newsgroup comp.os.msdos.djgpp
Reply-To: djgpp AT delorie DOT com

Damian Yerrick wrote:

>The FAQ entry describes the "rescaling" technique:
>
>#include <stdlib.h>
>
>int random_number =
>    low + (double)rand () * (high - low + 1) / RAND_MAX;

This should be

int random_number =
    low + (double)rand () * (high - low + 1) / (RAND_MAX+1.0);

>When I do this on my target platform (386 PC), the cast to double
>produces SIGNOFPE. [...] Does this mean I have to bundle emu387.dxe 
>with all programs that use random numbers?  

No. There are alternatives.

>Or can I do this to get rid of the unrandom low-order
>bits?
>  nextPiece = (rand() >> 4) % 5;

I suggest, to not use this technique. While this may yield acceptable
results with DJGPP rand(), it has some severe drawbacks. I.e., for 
portability, you have to assume, that RAND_MAX can be (and often is)
0x7fff. Then you'll get some numbers in the range with a probability
of 819/4096 and the other numbers with a probability of 820/4096.
(4096 = (RAND_MAX+1)>>4; 819 = floor(4096/5); 820 = ceil(4096/5).)
Also, when you change 5 to 4 or to 8, you can try the small program
I sent in this thread, and you'll see, that DJGPP rand() shows
failure with this method, in cases, it would do well with the
alternatives.

If you are not interested in portability (or only C99 and/or gcc)
the obvious way would be:

  nextPiece = (rand() * 5ULL)/(RAND_MAX+1U);

Unfortunately, this will produce rather inefficient code. For DJGPP
only, you can hope, that the value of RAND_MAX will never change.
Then you can use (with unsigned range=5U)

  nextPiece = ((unsigned long long)(unsigned)rand()*(range*2U))>>32;

This will be translated to one MUL and one ADD, and probably is the 
fastest you can get.

I once sent the KISS PRNG in a reply to your question. Here you
know the "RAND_MAX" for sure and for every platform (0xffffffffUL).
Then write

  nextPiece = ((unsigned long long)kiss_rand()*range)>>32;

One alternative, that should work with any C Implementation:

  nextPiece = rand()/(RAND_MAX/5 + 1);

All those methods are most probably sufficient for a range of 5,
but they don't generalize to a large range.

For a large range, the following should work with any C implementation.
It "developed" from an email correspondence with Hans-Bernhard Broeker.

#include <stdlib.h>

/* Call with 1 < range <= RAND_MAX+1, returns 0 <= r < range. */
unsigned rand_range(unsigned range)
{
  unsigned rmax, r, d;
  /* find the largest number rmax <= RAND_MAX, for which 
     (rmax+1) % range == 0.
     All returns from rand() > rmax will be skipped, to guarantee  
     equal probability for all return values. */
  /* The following depends on RAND_MAX + 1U <= UINT_MAX. */
  /* d = (RAND_MAX+1U) / range;  */
  d = (RAND_MAX+1U-range) / range + 1; 
  rmax = d * range - 1; /* -1 to avoid "overflow to zero" */
  do
    r = rand();
  while (r > rmax);
  return r/d;
}

You can change unsigned to unsigned long, RAND_MAX to
0xffffffffUL and rand() to kiss_rand(), if you like.

When you have to call rand_range often with the same range,
you can "cache" the values of rmax and d (or even calculate
them at compile time). Then the function will need on average 
one call to random, one compare and one division, when the range 
is small, and so it may even perform better, than 
(double)rand()/(RAND_MAX+1.0)*range. The worst case is 
range = (RAND_MAX+1)/2+1, where the function will need
two calls to rand() on average. (In this case all the other
methods would fail, because you get one number with only half the
probability of the other numbers.)

-- 
Regards, Dieter Buerssner