One issue that comes up fairly often around here is the poor
quality of the pseudo-random number generators supplied with many
C implementations. As a result, we have to recommend things like
using the high-order bits returned by rand() instead of the
low-order bits, avoiding using rand() for anything that wants
decently random numbers, not using rand() if you want more than
approx. UINT_MAX total different sequences, and so on.
So I wrote some PRNG code that uses RC4, with a seed of up to 256
bytes. Here it is. I believe it is clc-compliant. Comments on
this and anything else are welcome.
----------------------------------------------------------------------
/*
* prng.c - Portable, ISO C90 and C99 compliant high-quality
* pseudo-random number generator based on the alleged RC4
* cipher. This PRNG should be suitable for most general-purpose
* uses. Not recommended for cryptographic or financial
* purposes. Not thread-safe.
*/
/*
* Copyright (c) 2004 Ben Pfaff <bl*@cs.stanford.edu>.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or
* without modification, are permitted provided that the
* following conditions are met:
*
* 1. Redistributions of source code must retain the above
* copyright notice, this list of conditions and the following
* disclaimer.
*
* 2. Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials
* provided with the distribution.
*
* 3. The author's and contributors' names may not be used to
* endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS
* IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
* FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
* SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
* THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
* OF SUCH DAMAGE.
*
*/
#include "prng.h"
#include <assert.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <time.h>
/* RC4-based pseudo-random state. */
static unsigned char s[256];
static int s_i, s_j;
/* Nonzero if PRNG has been seeded. */
static int seeded;
/* Take advantage of inlining if possible. */
#if __STDC_VERSION__ >= 199901L
#define INLINE inline
#else
#define INLINE
#endif
/* Swap bytes. */
static INLINE void
swap_byte (unsigned char *a, unsigned char *b)
{
unsigned char t = *a;
*a = *b;
*b = t;
}
/* If KEY is non-null and SIZE is nonzero, seeds the
pseudo-random number generator based on the SIZE bytes in BUF.
At most the first 256 bytes of BUF are used. If KEY is a null
pointer and SIZE is zero, seeds the pseudo-random number
generator based on the current time.
If this function is not called by the user before any prng_get
function, it is called automatically to obtain a time-based
seed. */
void
prng_seed (const void *key_, size_t size)
{
const unsigned char *key = key_;
size_t key_idx;
int i, j;
assert ((key != NULL && size > 0) || (key == NULL && size == 0));
if (key == NULL)
{
static time_t t;
t = (t == 0) ? time (NULL) : t + 1;
key = (void *) &t;
size = sizeof t;
}
for (i = 0; i < 256; i++)
s[i] = i;
for (key_idx = 0, i = j = 0; i < 256; i++)
{
j = (j + s[i] + key[key_idx]) & 255;
swap_byte (s + i, s + j);
if (++key_idx >= size)
key_idx = 0;
}
s_i = s_j = 0;
seeded = 1;
}
/* Returns a pseudo-random integer in the range [0,255]. */
unsigned char
prng_get_octet (void)
{
if (!seeded)
prng_seed (NULL, 0);
s_i = (s_i + 1) & 255;
s_j = (s_j + s[s_i]) & 255;
swap_byte (s + s_i, s + s_j);
return s[(s[s_i] + s[s_j]) & 255];
}
/* Returns a pseudo-random integer in the range [0,UCHAR_MAX]. */
unsigned char
prng_get_byte (void)
{
if (UCHAR_MAX == 255)
{
/* Normal case for 8-bit bytes. */
return prng_get_octet ();
}
else
{
/* Handle oversize bytes. */
unsigned byte = prng_get_octet ();
unsigned done = 255;
while ((unsigned char) done != UCHAR_MAX)
{
byte = (byte << 8) | prng_get_octet ();
done = (done << 8) | 255;
}
return byte;
}
}
/* Fills BUF with SIZE pseudo-random bytes. */
void
prng_get_bytes (void *buf_, size_t size)
{
unsigned char *buf;
for (buf = buf_; size-- > 0; buf++)
*buf = prng_get_byte ();
}
/* Returns a pseudo-random unsigned long in the range [0,
ULONG_MAX]. */
unsigned long
prng_get_ulong (void)
{
unsigned long ulng = prng_get_octet ();
unsigned long done = 255;
while (done != ULONG_MAX)
{
ulng = (ulng << 8) | prng_get_octet ();
done = (done << 8) | 255;
}
return ulng;
}
/* Returns a pseudo-random long in the range [0, LONG_MAX]. */
long
prng_get_long (void)
{
for (;;)
{
unsigned ulng = prng_get_ulong ();
if (ulng <= LONG_MAX)
return ulng;
}
}
/* Returns a pseudo-random unsigned int in the range [0,
UINT_MAX]. */
unsigned
prng_get_uint (void)
{
unsigned uint = prng_get_octet ();
unsigned done = 255;
while (done != UINT_MAX)
{
uint = (uint << 8) | prng_get_octet ();
done = (done << 8) | 255;
}
return uint;
}
/* Returns a pseudo-random int in the range [0, INT_MAX]. */
int
prng_get_int (void)
{
for (;;)
{
unsigned uint = prng_get_uint ();
if (uint <= INT_MAX)
return uint;
}
}
/* Returns a pseudo-random floating-point number from the uniform
distribution with range [0,1). */
double
prng_get_double (void)
{
for (;;)
{
unsigned long ulng = prng_get_ulong ();
if (ulng != ULONG_MAX)
return (double) ulng / ULONG_MAX;
}
}
/* Returns a pseudo-random floating-point number from the
distribution with mean 0 and standard deviation 1. (Multiply
the result by the desired standard deviation, then add the
desired mean.) */
double
prng_get_double_normal (void)
{
/* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
Algorithm P. */
static double next_normal = DBL_MAX;
double this_normal;
if (next_normal != DBL_MAX)
{
this_normal = next_normal;
next_normal = DBL_MAX;
}
else
{
double v1, v2, s;
do
{
double u1 = prng_get_double ();
double u2 = prng_get_double ();
v1 = 2.0 * u1 - 1.0;
v2 = 2.0 * u2 - 1.0;
s = v1 * v1 + v2 * v2;
}
while (s >= 1);
this_normal = v1 * sqrt (-2. * log (s) / s);
next_normal = v2 * sqrt (-2. * log (s) / s);
}
return this_normal;
}
----------------------------------------------------------------------
#ifndef PRNG_H_INCLUDED
#define PRNG_H_INCLUDED
#include <stddef.h>
void prng_seed (const void *, size_t);
unsigned char prng_get_octet (void);
unsigned char prng_get_byte (void);
void prng_get_bytes (void *, size_t);
unsigned long prng_get_ulong (void);
long prng_get_long (void);
unsigned prng_get_uint (void);
int prng_get_int (void);
double prng_get_double (void);
double prng_get_double_normal (void);
#endif /* prng.h */
----------------------------------------------------------------------
--
"There's only one thing that will make them stop hating you.
And that's being so good at what you do that they can't ignore you.
I told them you were the best. Now you damn well better be."
--Orson Scott Card, _Ender's Game_
