Article View: comp.os.msdos.programmer
Article #5070Re: Random Numbers
From: Jesse Chisholm
Date: Thu, 07 Mar 1991 03:33
Date: Thu, 07 Mar 1991 03:33
599 lines
15355 bytes
15355 bytes
In article <r7bGl?wz@cs.psu.edu> melling@cs.psu.edu (Michael D Mellinger) writes:
>
>How can I get a random number on on DOS machines? Seeding the random
># generator with the system clock seemed like a good idea, but not all
>machines have clocks. Here is what I am doing now.
>
>srand( (unsigned int) time(theTime))
>
>Then I successively call rand() to get my random #.
>
>Anyone have a better idea?
>
>
>-Mike
Better only in the sense that the built in random number
generator in most compilers isn't very random.
The attached shar file is my implementation of a randum number
generator that is ver random and has a very long period before
the sequence of numbers repeat. There is also a test main
routine to make sure it is working correctly.
This random number generator does use the time of day if it must
as the seed. It also saves the current value in a file so it
can continue a sequence across program calls.
It is written in MicroSoft C 5.1 and should be portable to any
other language you are using.
===== jesse@Altos86.Altos.COM tel:1-408-432-6200x4810 fax:-434-0278
//// |||| |||| |||| || /
------------++++-----------++++- --------+++++--------++-------
|| //// | |||| | | /
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# mrn.c
# This archive created: Wed Mar 6 19:29:31 1991
export PATH; PATH=/bin:/usr/bin:$PATH
if test -f 'mrn.c'
then
echo shar: "will not over-write existing file 'mrn.c'"
else
cat << \SHAR_EOF > 'mrn.c'
/* ***************************************************
Purpose: This file has routines to implement a
Combined Prime Multiplicative Congruential
Psuedo-Random Number Generator.
Author: Jesse Chisholm
Auditor:
Creation Date: 91.03.04
Modifications:
*************************************************** */
/* include directives */
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
/* defines and macros */
#if !defined(DOS) && !defined(UNIX)
#error Must define one of: DOS, UNIX
#endif
#define NORMAL 4 /* this many uniforms at a time */
/*
** takes bits 30..15 of each long and glues them together
*/
#define V(a,b) ((((a) << 1) & 0xFFFF0000) | (((b) >> 14) & 0x0000FFFF))
/* global variables */
/*
** these two are visible to the outside world if they want
*/
long _i_seed1_ = 1L; /* current value of first generator */
long _i_seed2_ = 1L; /* current value of second generator */
/* static global variables */
/*
** these are all internal values
*/
static long _i_prim1_ = 2147483563L;
static long _i_prim2_ = 2147483399L;
static long _i_root1_ = 40014L;
static long _i_root2_ = 40692L;
static long _i_quo1_ = 53668L; /* _i_prim1_ / _i_root1_ */
static long _i_quo2_ = 52774L; /* _i_prim2_ / _i_root2_ */
static long _i_rem1_ = 12211L; /* _i_prim1_ % _i_root1_ */
static long _i_rem2_ = 3791L; /* _i_prim2_ % _i_root2_ */
/*
** The 10000th seed from 1L for arithmatic checking
*/
#define PRCHK1 1919456777L
#define PRCHK2 2006618587L
/* ***************************************************
Name: unsetran
Purpose: This routine saves the current values of the
two seeds on disk, so that the next program
that uses this generator can get the continued
sequence.
Sample Call: unsetran(); called automatically at program exit
Inputs: none
Outputs: none
may create a file on HOME directory
may change contents of file on HOME directory
Returns: none
Algorithm: attempt to open $HOME/.random UNIX
attempt to open $HOME\RANDOM.NUM DOS
write out the current values
close the file
*************************************************** */
static void unsetran()
{
FILE *fp;
char buf[128];
/*
** Attempt to locate the HOME random number file
*/
strcpy(buf, getenv("HOME"));
#ifdef DOS
if ((buf[strlen(buf)-1] != '\\') && (buf[strlen(buf)-1] != '/'))
strcat(buf, "\\");
strcat(buf, "RANDOM.NUM");
#endif /* DOS */
#ifdef UNIX
if (buf[strlen(buf)-1] != '/')
strcat(buf, "/");
strcat(buf, ".random");
#endif /* UNIX */
/*
** if we can find it, then write the new seed numbers out
*/
if ((fp = fopen(buf, "w")) != NULL) {
fprintf(fp, "%ld\n%ld\n", _i_seed1_, _i_seed2_);
fclose(fp);
}
}
/* ***************************************************
Name: setran
Purpose: This routine sets the random number
sequence at a known place. Or it can
continue a sequence saved by unsetran().
Sample Call: setran(1L,1L);
Inputs: long seed1; starting seed for 1st generator
long seed2; starting seed for 2nd generator
if seeds==0 try to use file
if seeds still ==0 then use time()
Outputs: none
Returns: none
Algorithm: ensure unsetran gets eventually called
if either seed is ==0L
try to open and read seed values
from the storage file
while either seed is ==0L
use time() value for seed
return combined seed value
*************************************************** */
void
setran(long seed1, long seed2)
{
long time();
/*
** ensure the values get saved when program is done
** this allows progression to continue each invokation
*/
atexit(unsetran);
/*
** if seeds==0L, attempt to read seed from file
*/
if ((seed1 == 0L) || (seed2 == 0L)) {
FILE *fp;
char buf[128];
strcpy(buf, getenv("HOME"));
#ifdef DOS
if ((buf[strlen(buf)-1] != '\\') && (buf[strlen(buf)-1] != '/'))
strcat(buf, "\\");
strcat(buf, "RANDOM.NUM");
#endif /* DOS */
#ifdef UNIX
if (buf[strlen(buf)-1] != '/')
strcat(buf, "/");
strcat(buf, ".random");
#endif /* UNIX */
if ((fp = fopen(buf, "r")) != NULL) {
fscanf(fp, "%ld\n%ld\n", &seed1, &seed2);
fclose(fp);
}
}
/*
** make sure neither seed is zero
*/
for(_i_seed1_=seed1, _i_seed2_=seed2;
(_i_seed1_==0L)||(_i_seed2_==0L);
_i_seed1_=time((long*)0)%_i_prim1_,_i_seed2_=time((long*)0)%_i_prim2_)
continue;
}
/* ***************************************************
Name: _rand
Purpose: This is the basic routine for generating
random numbers. The range is 00000000-FFFFFFFF.
The period is humongous.
Sample Call: l = _rand();
Inputs: none
Outputs: none
Returns: long l; return value is 32 bits wide
Algorithm: if either seed is ==0 call setran()
seed1 = seed1 * root1 % prime1
seed2 = seed2 * root2 % prime2
return combined most random bits
This implementation used the following numbers:
PRIME ROOT RANGE
2147483563 40014 1..2147483562
2147483399 40692 1..2147483398
these two used in conjunction give a total
range of 00000000..FFFFFFFF with a period
of 2147483562*2147483398 (I think)
*************************************************** */
/*
** NOTES ABOUT _RAND taken from earlier implementation
**
** _rand()
**
** returns UNIFORM random numbers in the range 1..2147483646
**
** This is the basic routine.
**
** new = old * root % prime;
**
** The table below shows some convenient values
**
** prime primitive root alternate root range
**
** 17 5 6 1..16
** 257 19 17 1..256
** 32749 182 128 1..32748
** 65537 32749 32719 1..65536
** 2147483647 16807 39373 1..2147483646
** 2147483647 48271 69621 1..2147483646
** 2147483647 397204094 630360016 1..2147483646
**
** This is the best single primitive root for 2^31-1 according to
** Pierre L'Ecuyer; "Random Numbers for Simulation"; CACM 10/90.
**
** 2147483647 41358 1..2147483646
**
** He also lists the prime 4611685301167870637
** and the primitive root 1968402271571654650
** as being pretty good.
**
** As a way of checking the arithmatic in any given implementation,
** the prime 2147483647L and the root 16807L with the initial seed
** of 1L produce the 10000th seed of 1043618065L
**
** There are many other primitive roots available for these
** or other primes. These were chosen as esthetically pleasing.
** Also chosen so that `x*y % p != 1' for roots x and y.
**
** A number (x) is a primitive root of a prime (p)
** iff in the equation:
**
** t = x^i % p
**
** when (i == p-1) and (t == 1)
**
** for all (1 <= i < p-1) then (t >= 2)
**
** The number of primitive roots for a given prime (p) is the number
** of numbers less that (p-1) that are relatively prime to (p-1).
**
** Note: that is how many, not which ones.
**
** if x is a primitive root of p, and y is relatively prime to p-1
** then (x^y % p) is also a primitive root of p.
** also (x^(p-1-y) % p) is a primitive root of p.
**
** The basic algorythm:
** new = old * root % prime;
** has the drawback that intermediate values are larger than the
** 4 byte register size. So an alternate way of implementing
** this expression is used so that no intermediate value exceeds
** the register size.
*/
long
_rand()
{
long test;
ldiv_t temp;
long alpha, beta, delta;
if ((_i_seed1_ == 0L) || (_i_seed2_ == 0L))
setran(0L,0L);
/*
** implement multiplication and modulo in such a way that
** intermediate values all fit in a long int.
*/
/* _i_seed_ = _i_seed_ * PROOT % PRIME */
temp = ldiv(_i_seed1_, _i_quo1_);
alpha = _i_root1_ * temp.rem;
beta = _i_rem1_ * temp.quot;
/*
** normalize the intermediate values
*/
if (alpha > beta) {
_i_seed1_ = alpha - beta;
} else {
_i_seed1_ = alpha - beta + _i_prim1_;
}
temp = ldiv(_i_seed2_, _i_quo2_);
alpha = _i_root2_ * temp.rem;
beta = _i_rem2_ * temp.quot;
if (alpha > beta) {
_i_seed2_ = alpha - beta;
} else {
_i_seed2_ = alpha - beta + _i_prim2_;
}
return( V(_i_seed1_, _i_seed2_) );
}
/* ***************************************************
Name: urand
Purpose: Returns uniform random numbers in a range.
Sample Call: ir = urand(low, high);
Inputs: int low; low end of range
int high; high end of range
Outputs: none
Returns: int ir; uniform random in that range
Algorithm: get a random number and scale it to the range
*************************************************** */
int urand(lo, hi)
int lo, hi;
{
long lr;
if (lo > hi)
return(urand(hi,lo));
if (lo == hi)
return(lo);
lr = _rand() & 0x7FFFFFFF; /* remove sign bit */
return( (int)((lr % ((long)(hi - lo + 1))) + (long)lo) );
}
/* ***************************************************
Name: nrand()
Purpose: returns NORMAL random numbers in a range
Sample Call: ir = nrand(low, high);
Inputs: int low; low end of range
int high; high end of range
Outputs: none
Returns: int ir; normal random integer in range
Algorithm: get some uniforms in the range and average them
*************************************************** */
int nrand(lo, hi)
int lo, hi;
{
long t;
int i;
for(t=i=0; i<NORMAL; i++)
t += (long)urand(lo,hi);
i = (int)((t + (NORMAL/2)) / NORMAL);
return(i);
}
#ifdef TEST
main()
{
int num;
long dd;
int limit;
printf("Testing validity of arithmatic in this implementation\n");
printf("of the EXTENDED STANDARD RANDOM NUMBER GENERATOR\n\n");
setran(1L,1L);
limit = 10000;
for(num=0; num<limit; num++) {
dd = _rand();
if ((num % 100) == 0) {
printf("\rIteration %d /10000", num);
}
}
printf("\rSeed should be %ld,%ld\n and is %ld,%ld.\n",
PRCHK1, PRCHK2, _i_seed1_, _i_seed2_);
if ((_i_seed1_ == PRCHK1) && (_i_seed2_ == PRCHK2)) {
printf ("Arithmatic is performed correctly.\n");
} else {
printf ("Arithmatic is performed INCORRECTLY!\n");
}
exit(0);
}
#endif /* TEST */
#ifdef GTEST
gotoxy(int x, int y)
{
printf("\033[%d;%dH", y, x);
}
clrscr()
{
printf("\033[2J");
}
main()
{
unsigned long count;
unsigned long SDcount;
int num;
long hits[32];
int xx,yy;
long dd;
char buf[20];
clrscr();
gotoxy(1,1); printf("Uniform Random Number Test");
gotoxy(1,23); printf("Total numbers tested");
gotoxy(34,23); printf(", Above midway");
for(num = 0; num < 32; num++) {
hits[num] = 0L;
gotoxy(10+num*2, 20); printf("0");
}
gotoxy(41, 21); printf("^");
setran(0L,0L);
SDcount = 0L;
for(count = 0L; !kbhit(); count++) {
num = urand(0,31);
if (num < 0) num = 0;
if (num > 31) num = 31;
hits[num]++;
xx = 10 + num * 2;
sprintf(buf, "%ld", hits[num]);
yy = 20;
for(dd=hits[num]; dd>0L; dd /= 10L) {
gotoxy(xx,yy--); printf("%ld", dd % 10L);
}
if (num > 15) SDcount++;
gotoxy(22,23); printf("%10lu", count);
gotoxy(48,23); printf("%10lu", SDcount);
if (count>0L) printf(" : %3d%%", (SDcount*100L)/count);
}
if (kbhit()) {
getch();
}
clrscr();
gotoxy(1,1); printf("Normal Random Number Test");
gotoxy(1,2);printf("Uniform numbers used %d at a time",
NORMAL);
gotoxy(1,23); printf("Total numbers tested");
gotoxy(34,23);printf(", Within 1 SD of center");
for(num = 0; num < 32; num++) {
hits[num] = 0L;
gotoxy(10+num*2, 20); printf("0");
}
gotoxy(41, 21); printf("^");
gotoxy(31, 21); printf("|");
gotoxy(51, 21); printf("|");
setran(0L,0L);
SDcount=0L;
for(count = 0L; !kbhit(); count++) {
num = nrand(-1,32);
if (num < 0) num = 0;
if (num > 31) num = 31;
hits[num]++;
xx = 10 + num * 2;
sprintf(buf, "%ld", hits[num]);
yy = 20;
for(dd=hits[num]; dd>0L; dd /= 10L) {
gotoxy(xx,yy--); printf("%ld", dd % 10L);
}
if ((11 <= num) && (num <= 20)) {
SDcount++;
}
gotoxy(22,23); printf("%10ld", count);
gotoxy(57,23); printf("%10ld", SDcount);
if (count>0L) printf(" : %3d%%", (SDcount*100L)/count);
}
if (kbhit()) {
getch();
}
clrscr();
gotoxy(1,23); printf("\n");
}
#endif /* GTEST */
SHAR_EOF
fi
exit 0
# End of shell archive
--
===== jesse@Altos86.Altos.COM tel:1-408-432-6200x4810 fax:-434-0278
"Why are you the way you are? And how long have you been that way?"
-- Tom Cafesjian
Technical Details
Message-ID:
<59@altos86.Altos.COM>Path:
arc1-ca.pugleaf.net!utzoo!news-server.csri.toronto.edu!cs.utexas.edu!usc!apple!altos!altos86!jesseReferences:
<r7bGl?wz@cs.psu.edu>Show Full Headers
Complete Article Headers
Path: arc1-ca.pugleaf.net!utzoo!news-server.csri.toronto.edu!cs.utexas.edu!usc!apple!altos!altos86!jesse From: jesse@altos86.Altos.COM ( Jesse Chisholm) Newsgroups: comp.os.msdos.programmer Subject: Re: Random Numbers Message-ID: <59@altos86.Altos.COM> Date: 7 Mar 91 03:33:39 GMT References: <r7bGl?wz@cs.psu.edu> Reply-To: jesse@altos86.UUCP (Acer - Jesse Chisholm) Distribution: comp Organization: Altos Computer Systems, San Jose, CA Lines: 599