Random number generator

From: James Turner (turnerjh@XTN.NET)
Date: 07/10/98 

I've watched this thread (as well as read the original slashdot
story).  Interesting stuff.  George is right, though, the GPL code
can't be incorporated.  Recently I redid the random number generator
in my code because I wanted it to be faster (well, I also wanted to
play with RNGs in general).  Below is a program to compare my code and
the circle code.  I base my algorithm pretty much verbatim from
Knuth's Art of Computer Programming volume 2 (Semi-numerical
algorithms).  It uses no multiplication and is quite fast, with a
tradeoff of maybe 220 bytes of memory.

It has a period of (at least) 2^55 - 1.  It can be adapted to generate
floating point random numbers as well... I've considered doing this
but have refrained.  In my code I needed numbers generated on a normal
distribution (ie the bell curve thing) and use floats for that
purpose.  Anyway, here's the code.  It's faster, has a longer period,
and has been proven in extensive tests (by others) to be quite
reliable.  My own experience indicates this as well.  Knuth writes
that he would recommend it completely except that though it has
withstood a great deal of empirical testing, little theoretical work
has made headway on it.

The output of the program without optimizations is:
Old: 51, new: 19

With -mpentium and -O9 on egcs 1.03a, the ouput is
Old: 21, new: 6

So this one is about 3 times faster with a longer period than
stock's.  Use it anywhere you like, I don't care, and you don't have
to put me in any credits or whatnot :)

One thing.  The code does use circle's old random number generator
when it seeds itself, so you'll need to keep it around.

--cut here--
#include <time.h>
#include <stdio.h>

#define m  (unsigned long)2147483647
#define q  (unsigned long)127773

#define a (unsigned int)16807
#define r (unsigned int)2836

/*
** F(z) = (az)%m
**      = az-m(az/m)
**
** F(z)  = G(z)+mT(z)
** G(z)  = a(z%q)- r(z/q)
** T(z)  = (z/q) - (az/m)
**
** F(z)  = a(z%q)- rz/q+ m((z/q) - a(z/m))
**       = a(z%q)- rz/q+ m(z/q) - az
*/

static unsigned long seed;

void circle_srandom_old(unsigned long initial_seed)
{
    seed = initial_seed;
}

unsigned long circle_random_old(void)
{
   register int lo, hi, test;

    hi   = seed/q;
    lo   = seed%q;

    test = a*lo - r*hi;

    if (test > 0)
        seed = test;
    else
        seed = test+ m;

    return seed;
}

/* new code below */

static unsigned long rand_list[56];
static unsigned long rand_j, rand_k;

void
circle_srandom_new(unsigned long seed)
{
  int i;

  circle_srandom_old(seed);

  for (i = 0; i <= 55; i++)
    rand_list[i] = circle_random_old();

  rand_list[24] |= 1;  /* make sure all aren't even */

  rand_j = 24;
  rand_k = 55;
}

unsigned long
circle_random_new()
{
  unsigned long retval;

  retval = rand_list[rand_k] = (rand_list[rand_k] + rand_list[rand_j]);

  rand_k--;
  rand_j--;

  if (rand_j == 0)
    rand_j = 55;
  else if (rand_k == 0)
    rand_k = 55;

  return retval;

}

#define NUM_TESTS     100000000

int
main()
{
  int i, delta1, delta2;

  circle_srandom_old(time(NULL));
  circle_srandom_new(time(NULL));

  delta1 = time(NULL);
  for (i = 0; i < NUM_TESTS; i++)
    circle_random_old();

  delta1 = time(NULL) - delta1;

  delta2 = time(NULL);
  for (i = 0; i < NUM_TESTS; i++)
    circle_random_new();

  delta2 = time(NULL) - delta2;

  printf("Old: %d, new: %d\n", delta1, delta2);
  return 0;
}

--cut here--


--
James Turner                turnerjh@pattern.net         UIN: 1102038
                            http://www.vuse.vanderbilt.edu/~turnerjh/


     +------------------------------------------------------------+
     | Ensure that you have read the CircleMUD Mailing List FAQ:  |
     | http://democracy.queensu.ca/~fletcher/Circle/list-faq.html |
     +------------------------------------------------------------+

This archive was generated by hypermail 2b30 : 12/15/00 PST