Hitmap 1.3
Main Page
Related Pages
Modules
Namespaces
Data Structures
Files
File List
Globals
All
Data Structures
Namespaces
Files
Functions
Variables
Typedefs
Friends
Macros
Groups
Pages
examples
NAS_IS
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
A
#define A
Definition:
mg.c:64
randlc
double randlc(double *X, double *A)
Definition:
c_randlc.c:40
Generated on Thu Oct 11 2018 12:23:26 for Hitmap 1.3 by
1.8.5