1 DOUBLE PRECISION FUNCTION DLARAN( ISEED )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Array Arguments ..
8 INTEGER ISEED( 4 )
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DLARAN returns a random real number from a uniform (0,1)
15 * distribution.
16 *
17 * Arguments
18 * =========
19 *
20 * ISEED (input/output) INTEGER array, dimension (4)
21 * On entry, the seed of the random number generator; the array
22 * elements must be between 0 and 4095, and ISEED(4) must be
23 * odd.
24 * On exit, the seed is updated.
25 *
26 * Further Details
27 * ===============
28 *
29 * This routine uses a multiplicative congruential method with modulus
30 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
31 * 'Multiplicative congruential random number generators with modulus
32 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
33 * b = 48', Math. Comp. 189, pp 331-344, 1990).
34 *
35 * 48-bit integers are stored in 4 integer array elements with 12 bits
36 * per element. Hence the routine is portable across machines with
37 * integers of 32 bits or more.
38 *
39 * =====================================================================
40 *
41 * .. Parameters ..
42 INTEGER M1, M2, M3, M4
43 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 )
44 DOUBLE PRECISION ONE
45 PARAMETER ( ONE = 1.0D+0 )
46 INTEGER IPW2
47 DOUBLE PRECISION R
48 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 )
49 * ..
50 * .. Local Scalars ..
51 INTEGER IT1, IT2, IT3, IT4
52 DOUBLE PRECISION RNDOUT
53 * ..
54 * .. Intrinsic Functions ..
55 INTRINSIC DBLE, MOD
56 * ..
57 * .. Executable Statements ..
58 10 CONTINUE
59 *
60 * multiply the seed by the multiplier modulo 2**48
61 *
62 IT4 = ISEED( 4 )*M4
63 IT3 = IT4 / IPW2
64 IT4 = IT4 - IPW2*IT3
65 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3
66 IT2 = IT3 / IPW2
67 IT3 = IT3 - IPW2*IT2
68 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2
69 IT1 = IT2 / IPW2
70 IT2 = IT2 - IPW2*IT1
71 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 +
72 $ ISEED( 4 )*M1
73 IT1 = MOD( IT1, IPW2 )
74 *
75 * return updated seed
76 *
77 ISEED( 1 ) = IT1
78 ISEED( 2 ) = IT2
79 ISEED( 3 ) = IT3
80 ISEED( 4 ) = IT4
81 *
82 * convert 48-bit integer to a real number in the interval (0,1)
83 *
84 RNDOUT = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
85 $ ( DBLE( IT4 ) ) ) ) )
86 *
87 IF (RNDOUT.EQ.1.0D+0) THEN
88 * If a real number has n bits of precision, and the first
89 * n bits of the 48-bit integer above happen to be all 1 (which
90 * will occur about once every 2**n calls), then DLARAN will
91 * be rounded to exactly 1.0.
92 * Since DLARAN is not supposed to return exactly 0.0 or 1.0
93 * (and some callers of DLARAN, such as CLARND, depend on that),
94 * the statistically correct thing to do in this situation is
95 * simply to iterate again.
96 * N.B. the case DLARAN = 0.0 should not be possible.
97 *
98 GOTO 10
99 END IF
100 *
101 DLARAN = RNDOUT
102 RETURN
103 *
104 * End of DLARAN
105 *
106 END
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Array Arguments ..
8 INTEGER ISEED( 4 )
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DLARAN returns a random real number from a uniform (0,1)
15 * distribution.
16 *
17 * Arguments
18 * =========
19 *
20 * ISEED (input/output) INTEGER array, dimension (4)
21 * On entry, the seed of the random number generator; the array
22 * elements must be between 0 and 4095, and ISEED(4) must be
23 * odd.
24 * On exit, the seed is updated.
25 *
26 * Further Details
27 * ===============
28 *
29 * This routine uses a multiplicative congruential method with modulus
30 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
31 * 'Multiplicative congruential random number generators with modulus
32 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
33 * b = 48', Math. Comp. 189, pp 331-344, 1990).
34 *
35 * 48-bit integers are stored in 4 integer array elements with 12 bits
36 * per element. Hence the routine is portable across machines with
37 * integers of 32 bits or more.
38 *
39 * =====================================================================
40 *
41 * .. Parameters ..
42 INTEGER M1, M2, M3, M4
43 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 )
44 DOUBLE PRECISION ONE
45 PARAMETER ( ONE = 1.0D+0 )
46 INTEGER IPW2
47 DOUBLE PRECISION R
48 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 )
49 * ..
50 * .. Local Scalars ..
51 INTEGER IT1, IT2, IT3, IT4
52 DOUBLE PRECISION RNDOUT
53 * ..
54 * .. Intrinsic Functions ..
55 INTRINSIC DBLE, MOD
56 * ..
57 * .. Executable Statements ..
58 10 CONTINUE
59 *
60 * multiply the seed by the multiplier modulo 2**48
61 *
62 IT4 = ISEED( 4 )*M4
63 IT3 = IT4 / IPW2
64 IT4 = IT4 - IPW2*IT3
65 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3
66 IT2 = IT3 / IPW2
67 IT3 = IT3 - IPW2*IT2
68 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2
69 IT1 = IT2 / IPW2
70 IT2 = IT2 - IPW2*IT1
71 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 +
72 $ ISEED( 4 )*M1
73 IT1 = MOD( IT1, IPW2 )
74 *
75 * return updated seed
76 *
77 ISEED( 1 ) = IT1
78 ISEED( 2 ) = IT2
79 ISEED( 3 ) = IT3
80 ISEED( 4 ) = IT4
81 *
82 * convert 48-bit integer to a real number in the interval (0,1)
83 *
84 RNDOUT = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
85 $ ( DBLE( IT4 ) ) ) ) )
86 *
87 IF (RNDOUT.EQ.1.0D+0) THEN
88 * If a real number has n bits of precision, and the first
89 * n bits of the 48-bit integer above happen to be all 1 (which
90 * will occur about once every 2**n calls), then DLARAN will
91 * be rounded to exactly 1.0.
92 * Since DLARAN is not supposed to return exactly 0.0 or 1.0
93 * (and some callers of DLARAN, such as CLARND, depend on that),
94 * the statistically correct thing to do in this situation is
95 * simply to iterate again.
96 * N.B. the case DLARAN = 0.0 should not be possible.
97 *
98 GOTO 10
99 END IF
100 *
101 DLARAN = RNDOUT
102 RETURN
103 *
104 * End of DLARAN
105 *
106 END