 Algorithmique Numérique Distribuée Public GIT Repository
1 /*
2  *    FUNCTION RANDLC (X, A)
3  *
4  *  This routine returns a uniform pseudorandom double precision number in the
5  *  range (0, 1) by using the linear congruential generator
6  *
7  *  x_{k+1} = a x_k  (mod 2^46)
8  *
9  *  where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
10  *  before repeating.  The argument A is the same as 'a' in the above formula,
11  *  and X is the same as x_0.  A and X must be odd double precision integers
12  *  in the range (1, 2^46).  The returned value RANDLC is normalized to be
13  *  between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
14  *  the new seed x_1, so that subsequent calls to RANDLC using the same
15  *  arguments will generate a continuous sequence.
16  *
17  *  This routine should produce the same results on any computer with at least
18  *  48 mantissa bits in double precision floating point data.  On Cray systems,
19  *  double precision should be disabled.
20  *
21  *  David H. Bailey     October 26, 1990
22  *
23  *     IMPLICIT DOUBLE PRECISION (A-H, O-Z)
24  *     SAVE KS, R23, R46, T23, T46
25  *     DATA KS/0/
26  *
27  *  If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
28  *  T23 = 2 ^ 23, and T46 = 2 ^ 46.  These are computed in loops, rather than
29  *  by merely using the ** operator, in order to insure that the results are
30  *  exact on all systems.  This code assumes that 0.5D0 is represented exactly.
31  */
32 double  randlc(double *X, double*A)
33 {
34   static int        KS=0;
35   static double  R23, R46, T23, T46;
36   double    T1, T2, T3, T4;
37   double    A1, A2;
38   double    X1, X2;
39   double    Z;
40   int       i, j;
42   if (KS == 0) {
43     R23 = 1.0;
44     R46 = 1.0;
45     T23 = 1.0;
46     T46 = 1.0;
48     for (i=1; i<=23; i++) {
49       R23 = 0.50 * R23;
50       T23 = 2.0 * T23;
51     }
52     for (i=1; i<=46; i++) {
53       R46 = 0.50 * R46;
54       T46 = 2.0 * T46;
55     }
56     KS = 1;
57   }
59 /*  Break A into two parts such that A = 2^23 * A1 + A2 and set X = N.  */
60   T1 = R23 * *A;
61   j  = T1;
62   A1 = j;
63   A2 = *A - T23 * A1;
65 /*  Break X into two parts such that X = 2^23 * X1 + X2, compute
66     Z = A1 * X2 + A2 * X1  (mod 2^23), and then X = 2^23 * Z + A2 * X2  (mod 2^46). */
67   T1 = R23 * *X;
68   j  = T1;
69   X1 = j;
70   X2 = *X - T23 * X1;
71   T1 = A1 * X2 + A2 * X1;
73   j  = R23 * T1;
74   T2 = j;
75   Z = T1 - T23 * T2;
76   T3 = T23 * Z + A2 * X2;
77   j  = R46 * T3;
78   T4 = j;
79   *X = T3 - T46 * T4;
80   return(R46 * *X);
81 }