From: Kbwms AT aol DOT com
Message-ID: <88406dd7.3653446f@aol.com>
Date: Wed, 18 Nov 1998 17:04:31 EST
To: djgpp-workers AT delorie DOT com
Mime-Version: 1.0
Subject: Fwd: Changes to Functions rand.c and random.c
Content-type: multipart/mixed;
	boundary="part0_911426672_boundary"
X-Mailer: AOL 3.0 16-bit for Windows sub 38
Reply-To: djgpp-workers AT delorie DOT com

This is a multi-part message in MIME format.

--part0_911426672_boundary
Content-ID: <0_911426672 AT inet_out DOT mail DOT aol DOT com DOT 1>
Content-type: text/plain; charset=US-ASCII


--part0_911426672_boundary
Content-ID: <0_911426672 AT inet_out DOT mail DOT aol DOT com DOT 2>
Content-type: message/rfc822
Content-transfer-encoding: 7bit
Content-disposition: inline

From: Kbwms AT aol DOT com
Return-path: <Kbwms AT aol DOT com>
To: dj AT delorie DOT com
Subject: Changes to Functions rand.c and random.c
Date: Wed, 18 Nov 1998 16:55:06 EST
Mime-Version: 1.0
Content-type: text/plain; charset=US-ASCII
Content-transfer-encoding: 7bit

Dear DJ Delorie,

In some recent tests, I discovered that function random() failed
the birthday spacings test due to George Marsaglia [1].  I have
modified the function accordingly using a method by Marsaglia [2].
After making the changes, I subjected the function to the scrutiny
of several other tests.  The full list of tests is:

    Name of test                Page in Knuth [1]
    ------------                -----------------
    birthday spacings test      71-72
    collision test              70-71
    frequency test              61
    maximum-of-t test           70
    permutation test            65-66
    run test                    66-69, 77 & 564 (ex. 14 & answer)
    serial test                 62
    serial correlation test     72-73

While I was in the area, I cranked a new version of rand() through
them, too.  As mentioned in a earlier e-mail, the original rand()
failed the spectral test and the collision test.  The change in
multiplier satisfied the spectral test but the generator still failed
the collision test in some dimensions.  I changed the right shift
count to 21.  Both functions rand() and random() now pass all tests.

Patches to update rand.c and random.c are appended.  The patches to
rand.c assume that the function will be restored to its simple form.


K.B. Williams

1.  Donald E. Knuth, *The Art of Computer*, Volume 2, Seminumerical
    Algorithms, Third Edition, Reading Massachusetts, Addison-Wesley
    Publishing Company.

2.  George Marsaglia and Arif Zaman, "Toward a universal random number
    generator," Statistics and Probability Letters 8(1990), pp 35-35,
    Elsevier Science Publishers B.V (North Holland)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*** rand.x      Sun Nov 15 14:15:08 1998
--- rand.c      Wed Nov 18 16:20:44 1998
***************
*** 2,38 ****
  /* Copyright (C) 1996 DJ Delorie, see COPYING.DJ for details */
  /* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */
  #include <stdlib.h>
- #include <pc.h>
- #include <go32.h>
- #include <libc/farptrgs.h>

! static unsigned long long next = 0;
! static int srand_called = 0;

  int
  rand(void)
  {
-   if (!srand_called)
-   {
-     unsigned int lsb, msb, tics;
-     outportb(0x43, 0x00);
-     lsb = inportb(0x40);
-     msb = inportb(0x40);
-     tics = _farpeekl(_dos_ds, 0x46c);
-     srand(lsb | (msb<<8) || (((long long)tics)<<16));
-   }
-
    next = next * 6364136223846793005LL + 1;
    /* was: next = next * 0x5deece66dLL + 11; */
!   return (int)((next >> 16) & RAND_MAX);
  }

  void
  srand(unsigned seed)
  {
    next = seed;
-   srand_called = 1;
-   rand();
-   rand();
-   rand();
  }
--- 2,20 ----
  /* Copyright (C) 1996 DJ Delorie, see COPYING.DJ for details */
  /* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */
  #include <stdlib.h>

! static unsigned long long next = 1;

  int
  rand(void)
  {
    next = next * 6364136223846793005LL + 1;
    /* was: next = next * 0x5deece66dLL + 11; */
!   return (int)((next >> 21) & RAND_MAX);
  }

  void
  srand(unsigned seed)
  {
    next = seed;
  }
*** random.x    Sun Oct  1 01:41:00 1995
--- random.c    Wed Nov 18 16:45:38 1998
***************
*** 3,8 ****
--- 3,9 ----
  /* This file may have been modified by DJ Delorie (Jan 1995).  If so,
  ** these modifications are Coyright (C) 1993 DJ Delorie, 24 Kirsten Ave,
  ** Rochester NH, 03867-2954, USA.
+ ** Also, modified by K.B. Williams (Nov 1998), kbwms AT aol DOT com
  */

  /*
***************
*** 141,148 ****
   * to point to randtbl[1] (as explained below).
   */

! static  long          *fptr                   = &randtbl[ SEP_3 + 1 ];
! static  long          *rptr                   = &randtbl[ 1 ];

  /*
   * The following things are the pointer to the state information table,
--- 142,149 ----
   * to point to randtbl[1] (as explained below).
   */

! static  long  *fptr           = (long *)&randtbl[ SEP_3 + 1 ];
! static  long  *rptr           = (long *)&randtbl[ 1 ];

  /*
   * The following things are the pointer to the state information table,
***************
*** 156,166 ****
   * the front and rear pointers have wrapped.
   */

! static  long          *state                  = &randtbl[ 1 ];
  static  int           rand_type               = TYPE_3;
  static  int           rand_deg                = DEG_3;
  static  int           rand_sep                = SEP_3;
! static  long          *end_ptr                = &randtbl[ DEG_3 + 1 ];

  /*
   * srandom:
--- 157,172 ----
   * the front and rear pointers have wrapped.
   */

! static  long  *state          = (long *)&randtbl[ 1 ];
  static  int   rand_type       = TYPE_3;
  static  int   rand_deg        = DEG_3;
  static  int   rand_sep        = SEP_3;
! static  long  *end_ptr        = (long *)&randtbl[ DEG_3 + 1 ];
!
! #define       Start_K1 735218L
! static        long    K1 = Start_K1;          /* Parameters for Arithmetic
*/
! static        long    K2 = 987654321;         /* Sequence Generator
*/
! static        long    M  = 2147483647;        /* Modulus = 2^31-1
*/

  /*
   * srandom:
***************
*** 178,184 ****
  int
  srandom(int x)
  {
!   int i, j;

    if (rand_type == TYPE_0)
    {
--- 184,190 ----
  int
  srandom(int x)
  {
!   int i;

    if (rand_type == TYPE_0)
    {
***************
*** 186,192 ****
    }
    else
    {
!     j = 1;
      state[ 0 ] = x;
      for (i = 1; i < rand_deg; i++)
      {
--- 192,200 ----
    }
    else
    {
!     /* initialize Arithmetic Sequence Generator */
!     K1 = Start_K1;
!
      state[ 0 ] = x;
      for (i = 1; i < rand_deg; i++)
      {
***************
*** 194,199 ****
--- 202,208 ----
      }
      fptr = &state[rand_sep];
      rptr = &state[0];
+     /* Warm up the generator */
      for( i = 0; i < 10*rand_deg; i++ )
        random();
    }
***************
*** 292,299 ****
  setstate(char *arg_state)
  {
    long *new_state = (long *)arg_state;
!   int type = new_state[0]%MAX_TYPES;
!   int rear = new_state[0]/MAX_TYPES;
    char *ostate = (char *)( &state[ -1 ] );

    if (rand_type == TYPE_0)
--- 301,308 ----
  setstate(char *arg_state)
  {
    long *new_state = (long *)arg_state;
!   int type = (int)new_state[0]%MAX_TYPES;
!   int rear = (int)new_state[0]/MAX_TYPES;
    char *ostate = (char *)( &state[ -1 ] );

    if (rand_type == TYPE_0)
***************
*** 326,336 ****
   * random:
   * If we are using the trivial TYPE_0 R.N.G., just do the old linear
   * congruential bit.  Otherwise, we do our fancy trinomial stuff, which is
the
!  * same in all ther other cases due to all the global variables that have
been
   * set up.  The basic operation is to add the number at the rear pointer
into
   * the one at the front pointer.  Then both pointers are advanced to the
next
!  * location cyclically in the table.  The value returned is the sum
generated,
!  * reduced to 31 bits by throwing away the "least random" low bit.
   * Note: the code takes advantage of the fact that both the front and
   * rear pointers can't wrap on the same call by not testing the rear
   * pointer if the front one has wrapped.
--- 335,350 ----
   * random:
   * If we are using the trivial TYPE_0 R.N.G., just do the old linear
   * congruential bit.  Otherwise, we do our fancy trinomial stuff, which is
the
!  * same in all the other cases due to all the global variables that have
been
   * set up.  The basic operation is to add the number at the rear pointer
into
   * the one at the front pointer.  Then both pointers are advanced to the
next
!  * location cyclically in the table.  The final value returned to the caller
!  * is the sum reduced to 31 bits by throwing away the "least random" low
bit.
!  * This value is combined with the output of a simple arithmetic sequence
!  * modulo 2^31-1 to create the final output.  This last step is required
!  * because all unvarnished lagged-Fibonacci generators fail the birthday
!  * spacings test.
!  *
   * Note: the code takes advantage of the fact that both the front and
   * rear pointers can't wrap on the same call by not testing the rear
   * pointer if the front one has wrapped.
***************
*** 340,355 ****
  long
  random(void)
  {
!   long i;

    if (rand_type == TYPE_0)
    {
!     i = state[0] = ( state[0]*1103515245 + 12345 )&0x7fffffff;
    }
    else
    {
      *fptr += *rptr;
!     i = (*fptr >> 1)&0x7fffffff; /* chucking least random bit */
      if (++fptr >= end_ptr )
      {
        fptr = state;
--- 354,373 ----
  long
  random(void)
  {
!   long Vout;

    if (rand_type == TYPE_0)
    {
!     Vout = state[0] = ( state[0]*1103515245 + 12345 )&0x7fffffff;
    }
    else
    {
      *fptr += *rptr;
!
!     Vout = (*fptr >> 1)&0x7fffffff; /* chucking least random bit */
!     K1 = (K1 - K2) > 0 ? (K1 - K2) : (K1 - K2) + M;
!     Vout = (Vout - K1) > 0 ? (Vout - K1) : (Vout - K1) + M;
!
      if (++fptr >= end_ptr )
      {
        fptr = state;
***************
*** 361,365 ****
        rptr = state;
      }
    }
!   return i;
  }
--- 379,383 ----
        rptr = state;
      }
    }
!   return Vout;
  }

--part0_911426672_boundary--