/*							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 union { double d; unsigned short s[4];} unkans;/* This function implements the three * congruential generators. */ int 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. */return 0;}/*	drand.c * * Random double precision floating point number between 1 and 2. * * C callable: *	drand( &x ); */int drand( a )double *a;{unsigned short r;#ifdef DECunsigned short s, t;#endif/* This algorithm of Wichmann and Hill computes a floating point * result: */ranwh();unkans.d = sx/30269.0  +  sy/30307.0  +  sz/30323.0;r = unkans.d;unkans.d -= r;unkans.d += 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 */#elseranwh();r = sx * sy + sz;#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 = unkans.s[2];t = s & 07;	/* save these bits to put in at the bottom */s &= 0177770;s |= (r >> 13) & 07;unkans.s[2] = s;t |= r << 3;unkans.s[3] = t;#endif#ifdef IBMPCunkans.s[0] = r;#endif#ifdef MIEEEunkans.s[3] = r;#endif*a = unkans.d;return 0;}