1/* 2*/ 3#if defined(USE_POW) 4#define r23 pow(0.5, 23.0) 5#define r46 (r23*r23) 6#define t23 pow(2.0, 23.0) 7#define t46 (t23*t23) 8#else 9#define r23 (0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5) 10#define r46 (r23*r23) 11#define t23 (2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0) 12#define t46 (t23*t23) 13#endif 14 15/*c--------------------------------------------------------------------- 16c---------------------------------------------------------------------*/ 17 18double randlc (double *x, double a) { 19 20/*c--------------------------------------------------------------------- 21c---------------------------------------------------------------------*/ 22 23/*c--------------------------------------------------------------------- 24c 25c This routine returns a uniform pseudorandom double precision number in the 26c range (0, 1) by using the linear congruential generator 27c 28c x_{k+1} = a x_k (mod 2^46) 29c 30c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers 31c before repeating. The argument A is the same as 'a' in the above formula, 32c and X is the same as x_0. A and X must be odd double precision integers 33c in the range (1, 2^46). The returned value RANDLC is normalized to be 34c between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain 35c the new seed x_1, so that subsequent calls to RANDLC using the same 36c arguments will generate a continuous sequence. 37c 38c This routine should produce the same results on any computer with at least 39c 48 mantissa bits in double precision floating point data. On 64 bit 40c systems, double precision should be disabled. 41c 42c David H. Bailey October 26, 1990 43c 44c---------------------------------------------------------------------*/ 45 46 double t1,t2,t3,t4,a1,a2,x1,x2,z; 47 48/*c--------------------------------------------------------------------- 49c Break A into two parts such that A = 2^23 * A1 + A2. 50c---------------------------------------------------------------------*/ 51 t1 = r23 * a; 52 a1 = (int)t1; 53 a2 = a - t23 * a1; 54 55/*c--------------------------------------------------------------------- 56c Break X into two parts such that X = 2^23 * X1 + X2, compute 57c Z = A1 * X2 + A2 * X1 (mod 2^23), and then 58c X = 2^23 * Z + A2 * X2 (mod 2^46). 59c---------------------------------------------------------------------*/ 60 t1 = r23 * (*x); 61 x1 = (int)t1; 62 x2 = (*x) - t23 * x1; 63 t1 = a1 * x2 + a2 * x1; 64 t2 = (int)(r23 * t1); 65 z = t1 - t23 * t2; 66 t3 = t23 * z + a2 * x2; 67 t4 = (int)(r46 * t3); 68 (*x) = t3 - t46 * t4; 69 70 return (r46 * (*x)); 71} 72 73/*c--------------------------------------------------------------------- 74c---------------------------------------------------------------------*/ 75 76void vranlc (int n, double *x_seed, double a, double y[]) { 77 78/*c--------------------------------------------------------------------- 79c---------------------------------------------------------------------*/ 80 81/*c--------------------------------------------------------------------- 82c 83c This routine generates N uniform pseudorandom double precision numbers in 84c the range (0, 1) by using the linear congruential generator 85c 86c x_{k+1} = a x_k (mod 2^46) 87c 88c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers 89c before repeating. The argument A is the same as 'a' in the above formula, 90c and X is the same as x_0. A and X must be odd double precision integers 91c in the range (1, 2^46). The N results are placed in Y and are normalized 92c to be between 0 and 1. X is updated to contain the new seed, so that 93c subsequent calls to VRANLC using the same arguments will generate a 94c continuous sequence. If N is zero, only initialization is performed, and 95c the variables X, A and Y are ignored. 96c 97c This routine is the standard version designed for scalar or RISC systems. 98c However, it should produce the same results on any single processor 99c computer with at least 48 mantissa bits in double precision floating point 100c data. On 64 bit systems, double precision should be disabled. 101c 102c---------------------------------------------------------------------*/ 103 104 int i; 105 double x,t1,t2,t3,t4,a1,a2,x1,x2,z; 106 107/*c--------------------------------------------------------------------- 108c Break A into two parts such that A = 2^23 * A1 + A2. 109c---------------------------------------------------------------------*/ 110 t1 = r23 * a; 111 a1 = (int)t1; 112 a2 = a - t23 * a1; 113 x = *x_seed; 114 115/*c--------------------------------------------------------------------- 116c Generate N results. This loop is not vectorizable. 117c---------------------------------------------------------------------*/ 118 for (i = 1; i <= n; i++) { 119 120/*c--------------------------------------------------------------------- 121c Break X into two parts such that X = 2^23 * X1 + X2, compute 122c Z = A1 * X2 + A2 * X1 (mod 2^23), and then 123c X = 2^23 * Z + A2 * X2 (mod 2^46). 124c---------------------------------------------------------------------*/ 125 t1 = r23 * x; 126 x1 = (int)t1; 127 x2 = x - t23 * x1; 128 t1 = a1 * x2 + a2 * x1; 129 t2 = (int)(r23 * t1); 130 z = t1 - t23 * t2; 131 t3 = t23 * z + a2 * x2; 132 t4 = (int)(r46 * t3); 133 x = t3 - t46 * t4; 134 y[i] = r46 * x; 135 } 136 *x_seed = x; 137} 138