/*-------------------------------------------------------------------------*/
/* Uniform random numbers: Multiplicative congruential generator for */
/* uniformly distributed random numbers on (0, 1) */
/* */
/* Modul m = 2**54, Multiplier a = 2.783.377.640.906.189 */
/* Period = 2**52, = a1 * 2**27 + a2 */
/* where a1 = 20.737.779, a2 = 59.760.077 */
/*-------------------------------------------------------------------------*/
/* */
/* The multiplier a of the generator was the winner of a search among one */
/* million possible multipliers near to the well-known golden section */
/* factor f = .5*(sqrt(5)-1)*2^52. The multiplier a was choosen because of */
/* its good lattice properties up to dimension 8 measured by both Beyer */
/* quotients (Afflerbach/Grothe,1985, Computing 35, 269-276) and */
/* standardized values of the spectral test according to Fishman (1990), */
/* Math.Comput. 54,331-344. The search was carried out by J.H. Ahrens(1992).*/
/* The generator (denoted as krand) was implemented and tested as described*/
/* in Stadlober/Kremer (1992), Lecture Notes in Econ.Math.Systems 374, */
/* 154-162. */
/* The values of the measures from dimension 2 to dimension 8 are: */
/* Beyer-quotients: .87, .71, .72, .73, .66, .82, .80; Minimum = .6626 */
/* Spectral-test: .87, .76, .79, .78, .75, .73, .69; Minimum = .6870 */
/* */
/*-------------------------------------------------------------------------*/
/* */
/* The seed-values z_hi and z_lo has to be defined as follows: */ */
/* printf("Input seed-high (0 ... 2^27 - 1) --> "); */
/* scanf("%ld", &z_hi); */
/* printf("Input seed-low (0 ... 2^25 - 1) --> "); */
/* scanf("%ld", &z_lo); z_lo = z_lo*4L + 1L; */
/*-------------------------------------------------------------------------*/
#include
double drand(unsigned long int *z_hi, unsigned long int *z_lo)
{
union ux {
double d;
int i[4];
} x;
unsigned long int i_x;
double d_x;
x.d = (double)(*z_lo) * 59760077.0; x.i[3] -= 0X01b0;
i_x = x.d;
*z_hi = (*z_hi * 59760077L + *z_lo * 20737779L + i_x) & 0X07FFFFFFL;
d_x = x.d -= (double)i_x; x.i[3] += 0X01b0;
*z_lo = x.d;
x.d = (double)*z_hi + d_x; x.i[3] -= 0X01b0;
return (x.d);
}