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