Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : DRAND.C
*
* Pseudorandom number generator
*
*
*
* SYNOPSIS:
*
* double y, drand();
*
* drand( &y );
*
*
*
* DESCRIPTION:
*
* Yields a random number 1.0 <= y < 2.0.
*
* The three-generator congruential algorithm by Brian
* Wichmann and David Hill (BYTE magazine, March, 1987,
* pp 127-8) is used. The period, given by them, is
* 6953607871644.
*
* Versions invoked by the different arithmetic compile
* time options DEC, IBMPC, and MIEEE, produce
* approximately the same sequences, differing only in the
* least significant bits of the numbers. The UNK option
* implements the algorithm as recommended in the BYTE
* article. It may be used on all computers. However,
* the low order bits of a double precision number may
* not be adequately random, and may vary due to arithmetic
* implementation details on different computers.
*
* The other compile options generate an additional random
* integer that overwrites the low order bits of the double
* precision number. This reduces the period by a factor of
* two but tends to overcome the problems mentioned.
*
*/
#include "mconf.h"
/* Three-generator random number algorithm
* of Brian Wichmann and David Hill
* BYTE magazine, March, 1987 pp 127-8
*
* The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
* This program was written by Steve Moshier.
*/
static int sx = 1;
static int sy = 10000;
static int sz = 3000;
static double unkans = 0.0;
/* This function implements the three
* congruential generators.
*/
ranwh()
{
int r, s;
/* sx = sx * 171 mod 30269 */
r = sx/177;
s = sx - 177 * r;
sx = 171 * s - 2 * r;
if( sx < 0 )
sx += 30269;
/* sy = sy * 172 mod 30307 */
r = sy/176;
s = sy - 176 * r;
sy = 172 * s - 35 * r;
if( sy < 0 )
sy += 30307;
/* sz = 170 * sz mod 30323 */
r = sz/178;
s = sz - 178 * r;
sz = 170 * s - 63 * r;
if( sz < 0 )
sz += 30323;
/* The results are in static sx, sy, sz. */
}
/* drand.c
*
* Random double precision floating point number between 1 and 2.
*
* C callable:
* drand( &x );
*/
drand( a )
double *a;
{
unsigned short r;
unsigned short *p;
#ifdef DEC
unsigned short s, t;
#endif
/* This algorithm of Wichmann and Hill computes a floating point
* result:
*/
ranwh();
unkans = sx/30269.0 + sy/30307.0 + sz/30323.0;
r = unkans;
unkans -= r;
unkans += 1.0;
/* if UNK option, do nothing further.
* Otherwise, make a random 16 bit integer
* to overwrite the least significant word
* of unkans.
*/
#ifdef UNK
/* do nothing */
#else
ranwh();
r = sx * sy + sz;
p = (unsigned short *)&unkans;
#endif
#ifdef DEC
/* To make the numbers as similar as possible
* in all arithmetics, the random integer has
* to be inserted 3 bits higher up in a DEC number.
* An alternative would be put it 3 bits lower down
* in all the other number types.
*/
s = *(p + 2);
t = s & 07; /* save these bits to put in at the bottom */
s &= 0177770;
s |= (r >> 13) & 07;
*(p + 2) = s;
t |= r << 3;
*(p + 3) = t;
#endif
#ifdef IBMPC
*p = r;
#endif
#ifdef MIEEE
*(p + 3) = r;
#endif
*a = unkans;
}
/* drand.c
*
* Pseudorandom number generator
*
*
*
* SYNOPSIS:
*
* double y, drand();
*
* drand( &y );
*
*
*
* DESCRIPTION:
*
* Yields a random number 1.0 <= y < 2.0.
*
* The three-generator congruential algorithm by Brian
* Wichmann and David Hill (BYTE magazine, March, 1987,
* pp 127-8) is used. The period, given by them, is
* 6953607871644.
*
* Versions invoked by the different arithmetic compile
* time options DEC, IBMPC, and MIEEE, produce
* approximately the same sequences, differing only in the
* least significant bits of the numbers. The UNK option
* implements the algorithm as recommended in the BYTE
* article. It may be used on all computers. However,
* the low order bits of a double precision number may
* not be adequately random, and may vary due to arithmetic
* implementation details on different computers.
*
* The other compile options generate an additional random
* integer that overwrites the low order bits of the double
* precision number. This reduces the period by a factor of
* two but tends to overcome the problems mentioned.
*
*/
#include "mconf.h"
/* Three-generator random number algorithm
* of Brian Wichmann and David Hill
* BYTE magazine, March, 1987 pp 127-8
*
* The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
*/
static int sx = 1;
static int sy = 10000;
static int sz = 3000;
static double unkans = 0.0;
/* This function implements the three
* congruential generators.
*/
ranwh()
{
int r, s;
/* sx = sx * 171 mod 30269 */
r = sx/177;
s = sx - 177 * r;
sx = 171 * s - 2 * r;
if( sx < 0 )
sx += 30269;
/* sy = sy * 172 mod 30307 */
r = sy/176;
s = sy - 176 * r;
sy = 172 * s - 35 * r;
if( sy < 0 )
sy += 30307;
/* sz = 170 * sz mod 30323 */
r = sz/178;
s = sz - 178 * r;
sz = 170 * s - 63 * r;
if( sz < 0 )
sz += 30323;
/* The results are in static sx, sy, sz. */
}
/* drand.c
*
* Random double precision floating point number between 1 and 2.
*
* C callable:
* drand( &x );
*/
drand( a )
double *a;
{
unsigned short r;
unsigned short *p;
#ifdef DEC
unsigned short s, t;
#endif
/* This algorithm of Wichmann and Hill computes a floating point
* result:
*/
ranwh();
unkans = sx/30269.0 + sy/30307.0 + sz/30323.0;
r = unkans;
unkans -= r;
unkans += 1.0;
/* if UNK option, do nothing further.
* Otherwise, make a random 16 bit integer
* to overwrite the least significant word
* of unkans.
*/
#ifdef UNK
/* do nothing */
#else
ranwh();
r = sx * sy + sz;
p = (unsigned short *)&unkans;
#endif
#ifdef DEC
/* To make the numbers as similar as possible
* in all arithmetics, the random integer has
* to be inserted 3 bits higher up in a DEC number.
* An alternative would be put it 3 bits lower down
* in all the other number types.
*/
s = *(p + 2);
t = s & 07; /* save these bits to put in at the bottom */
s &= 0177770;
s |= (r >> 13) & 07;
*(p + 2) = s;
t |= r << 3;
*(p + 3) = t;
#endif
#ifdef IBMPC
*p = r;
#endif
#ifdef MIEEE
*(p + 3) = r;
#endif
*a = unkans;
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/