Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
c_randlc.c
Go to the documentation of this file.
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 
33 
34 /*****************************************************************/
35 /************* R A N D L C ************/
36 /************* ************/
37 /************* portable random number generator ************/
38 /*****************************************************************/
39 
40 double randlc( double *X, double *A )
41 {
42  static int KS=0;
43  static double R23, R46, T23, T46;
44  double T1, T2, T3, T4;
45  double A1;
46  double A2;
47  double X1;
48  double X2;
49  double Z;
50  int i, j;
51 
52  if (KS == 0)
53  {
54  R23 = 1.0;
55  R46 = 1.0;
56  T23 = 1.0;
57  T46 = 1.0;
58 
59  for (i=1; i<=23; i++)
60  {
61  R23 = 0.50 * R23;
62  T23 = 2.0 * T23;
63  }
64  for (i=1; i<=46; i++)
65  {
66  R46 = 0.50 * R46;
67  T46 = 2.0 * T46;
68  }
69  KS = 1;
70  }
71 
72 /* Break A into two parts such that A = 2^23 * A1 + A2 and set X = N. */
73 
74  T1 = R23 * *A;
75  j = (int) T1;
76  A1 = j;
77  A2 = *A - T23 * A1;
78 
79 /* Break X into two parts such that X = 2^23 * X1 + X2, compute
80  Z = A1 * X2 + A2 * X1 (mod 2^23), and then
81  X = 2^23 * Z + A2 * X2 (mod 2^46). */
82 
83  T1 = R23 * *X;
84  j = (int) T1;
85  X1 = j;
86  X2 = *X - T23 * X1;
87  T1 = A1 * X2 + A2 * X1;
88 
89  j = (int) (R23 * T1);
90  T2 = j;
91  Z = T1 - T23 * T2;
92  T3 = T23 * Z + A2 * X2;
93  j = (int) (R46 * T3);
94  T4 = j;
95  *X = T3 - T46 * T4;
96  return(R46 * *X);
97 }
98 
#define A
Definition: mg.c:64
double randlc(double *X, double *A)
Definition: c_randlc.c:40