random numbers

[D. Jason Penney]

In article <9849 at alice.UUCP> andrew at alice.UUCP (Andrew Hume) writes:
>
>
>if you can tolerate a congruential style generator, try this:
>

If you can NOT tolerate a congruential style generator, try this:
(sorry, no test program or main -- shouldn't be a problem, though...)

/*========================================================================
*
* Name - %M% - random.hf
*
* Version:	%I%
*
* ccsid:	%W% - %G% %U%
* from: 	%F%
* date: 	%H% %T%
*
* %Q%
*  
* Description: 
*

*========================================================================*/
void RandomInit(void);

unsigned int Random(void);

unsigned int RandomBetween__and__(unsigned int lower, unsigned int upper);

/*========================================================================
*
* Name - %M% -- random.c
*
* Version:	%I%
*
* ccsid:	%W% - %G% %U%
* from: 	%F%
* date: 	%H% %T%
*
* %Q%
*  
* Description: Random Number Generator
*
*========================================================================*/

#include <random.hf>
#include <time.h>

/* This size must be power of two because of strange arithmetic we
 must go through: */
#define auxRandomTableSize 64

static unsigned int auxRandomTable[auxRandomTableSize];
static unsigned int
	randomSeed, /* for repeated operation */
	knuth_shuffler /* for Algorithm D */;

static unsigned int nextCyclicRandom(void)
/*
  Note: according to Knuth volume 2, this generator has some serious
  flaws, depending on how it is used.  We incorporate this generator with
  another randomization scheme to produce FUNCTION RANDOM.
*/
{
/* Note we take advantage of C arithmetic not to signal arithmetic overflow. */
 randomSeed = (unsigned int)((13849L + 27181L * randomSeed) & 65535L);
 return randomSeed;
}

void RandomInit(void)
/* random number initialization */
{
int i;

randomSeed = (unsigned int)(clock());
for (i = 0; i < auxRandomTableSize; i ++)
  auxRandomTable[i] = nextCyclicRandom();
knuth_shuffler = nextCyclicRandom();
}

/* Random Number Generation */

unsigned int Random(void)
/* See Knuth volume 2, page 32, "Algorithm B" */
{
int j;

/* B1 */
  j = (int)(knuth_shuffler * auxRandomTableSize / 65536);

/* B2 */
  knuth_shuffler = auxRandomTable[j];
  auxRandomTable[j] = nextCyclicRandom();
  return knuth_shuffler;
}

unsigned int RandomBetween__and__(unsigned int lower, unsigned int upper)
{
unsigned int modulus, threshold, probe;

if ((upper == 65535) && (lower == 0)) /* easy */
  return Random();

modulus = upper - lower + 1;
threshold = (int)(65535 / modulus * modulus);
do {
  probe = Random();
  } while (probe >= threshold); /* to make perfectly uniform */
return lower + probe % modulus;
}
-- 
D. Jason Penney                  Ph: (503) 629-8383
Beaverton, OR 97006              uucp: ...uunet!servio!penneyj
STANDARD DISCLAIMER:  Should I or my opinions be caught or killed, the
company will disavow any knowledge of my actions...