1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
REAL FUNCTION URAND(IY)
INTEGER IY
C
C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND
C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2. THE INTEGER IY
C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL
C TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY
C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED
C IN THE INTERVAL (0,1).
C
INTEGER IA,IC,ITWO,M2,M,MIC
DOUBLE PRECISION HALFM
REAL S
DOUBLE PRECISION DATAN,DSQRT
DATA M2/0/,ITWO/2/
IF (M2 .NE. 0) GO TO 20
C
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH
C
M = 1
10 M2 = M
M = ITWO*M2
IF (M .GT. M2) GO TO 10
HALFM = M2
C
C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
C
IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5
IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1
MIC = (M2 - IC) + M2
C
C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
C
S = 0.5/HALFM
C
C COMPUTE NEXT RANDOM NUMBER
C
20 IY = IY*IA
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW
C INTEGER OVERFLOW ON ADDITION
C
IF (IY .GT. MIC) IY = (IY - M2) - M2
C
IY = IY + IC
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE
C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION
C
IF (IY/2 .GT. M2) IY = (IY - M2) - M2
C
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER
C OVERFLOW AFFECTS THE SIGN BIT
C
IF (IY .LT. 0) IY = (IY + M2) + M2
URAND = FLOAT(IY)*S
RETURN
END
|
