Mail Archives: djgpp-workers/1998/11/18/17:08:34
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--
- Raw text -