(fwd) Random Numbers - Results of testing BSD random()

[Jim choate <ravage@bga.com>]

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis
Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd
From: deleyd@netcom.com
Subject: Random Numbers - Results of testing BSD random()
Message-ID: <deleydCsIAyv.K6n@netcom.com>
Organization: NETCOM On-line Communication Services (408 261-4700 guest)
Date: Wed, 6 Jul 1994 06:49:42 GMT
Lines: 119
Xref: bga.com sci.stat.math:1314 sci.math:15352 sci.math.num-analysis:3353

BSD random()

Here are the partial results.  Further tests were not performed due to
lack of time.  So far the generator appears to be comparable to a
shuffled linear congruential generator.

          DEFINITION:
            Generating polynomial: x^31 + x^3 + 1 (primitive polynomial)
            Initialize circular queue of 31 elements using ANSI C linear
             congruential generator.
            Recursion formula: a[i] = a[i] + a[i-3]

          RATING:
             1-D FAILS above 800,000 bpd      (bins per dimension)
             2-D FAILS above 3000 bpd
             3-D FAILS above 210 bpd
             4-D PASSES at 50 bpd (highest tested so far)
             5-D not tested
             6-D not tested
             7-D not tested
             8-D not tested

This is an additive congruential type random number generator.  An array
table[31] is initially filled with random numbers using the ANSI C
linear congruential random number generator.  Random numbers are then
generated using the recursion formula:

     table[k] = (table[k-31] + table[k-3]) mod 32

(Note that x**31 + x**3 + 1 is a primitive polynomial mod 2, which is
being used here as a generator.)  Since we are using the array table[]
as a circular queue with 31 elements then table[k-31] is just table[k]
before it gets replaced with the new value.  The recursion formula
becomes:

     table[k] = table[k] + table[k-3]

The generator works well in practice.  Knuth claims the sequence will
have period 2**31 - 1.  Knuth also claims there is very little theory to
prove that this generator does or does not have desirable random
properties.  I would be interested if anyone knows of any recent
developments in this area.

-David Deley
deleyd@netcom.com

(So sorry, I lost the name of the original person who posted this code below
 which was used in the tests. -D.D.)

/***
   Code to implement random() & srandom() of BSD Unix. It was taken
   (though coded somewhat differently) from the Gnu BSD implementation.
 ***/

#include <stdio.h>
#include <stdlib.h>
#define LONG31

#ifdef LONG31  /* x^31 + x^3 + 1 */
#define SIZE  31
#define SIZE1 30
#define P1 3
#define P2 0
#else  /* LONG63: x^63 + x + 1 */
#define SIZE  63
#define SIZE1 62
#define P1 1
#define P2 0
#endif

#define LONG_MAX  0x7fffffff


int p1=P1, p2=P2;  
long table[SIZE];

/*** return a "random" number in range [0, LONG_MAX] */

long xrand ()
{
  int r;

  table[p1] = table[p1] + table[p2]; /* add two table elements */
  r = (table[p1] >> 1) & LONG_MAX;   /* throw least significant bit away */

  if (p1 == SIZE1) { /* increment the table indexes */
    p1 = 0; 
    p2 = p2 + 1; 
  }
  else if (p2 == SIZE1) {
    p1 = p1 + 1;    
    p2 = 0;
  }
  else {
    p1 = p1 + 1;    
    p2 = p2 + 1;
  }

  return (r);
}


/*** use a linear congruential type generator to seed the 
     state table & cycle the entire table 10 times */

void sxrand (seed)
long seed;
{
  int i;

  table[0] = seed;
  for (i=1; i<SIZE; ++i)
    table[i] = (table[i-1] * 1103515145) + 12345;  /* lousy */

  for (i=0; i<10*SIZE; ++i)
    (void) xrand();
}