/*-------------------------------------------------------------------------*/
/*  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 <math.h>

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);
}