1 REAL FUNCTION SLARAN( 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 * SLARAN 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 REAL ONE
45 PARAMETER ( ONE = 1.0E+0 )
46 INTEGER IPW2
47 REAL R
48 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 )
49 * ..
50 * .. Local Scalars ..
51 INTEGER IT1, IT2, IT3, IT4
52 REAL RNDOUT
53 * ..
54 * .. Intrinsic Functions ..
55 INTRINSIC MOD, REAL
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*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R*
85 $ ( REAL( IT4 ) ) ) ) )
86 *
87 IF (RNDOUT.EQ.1.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 SLARAN will
91 * be rounded to exactly 1.0. In IEEE single precision arithmetic,
92 * this will happen relatively often since n = 24.
93 * Since SLARAN is not supposed to return exactly 0.0 or 1.0
94 * (and some callers of SLARAN, such as CLARND, depend on that),
95 * the statistically correct thing to do in this situation is
96 * simply to iterate again.
97 * N.B. the case SLARAN = 0.0 should not be possible.
98 *
99 GOTO 10
100 END IF
101 *
102 SLARAN = RNDOUT
103 RETURN
104 *
105 * End of SLARAN
106 *
107 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 * SLARAN 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 REAL ONE
45 PARAMETER ( ONE = 1.0E+0 )
46 INTEGER IPW2
47 REAL R
48 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 )
49 * ..
50 * .. Local Scalars ..
51 INTEGER IT1, IT2, IT3, IT4
52 REAL RNDOUT
53 * ..
54 * .. Intrinsic Functions ..
55 INTRINSIC MOD, REAL
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*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R*
85 $ ( REAL( IT4 ) ) ) ) )
86 *
87 IF (RNDOUT.EQ.1.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 SLARAN will
91 * be rounded to exactly 1.0. In IEEE single precision arithmetic,
92 * this will happen relatively often since n = 24.
93 * Since SLARAN is not supposed to return exactly 0.0 or 1.0
94 * (and some callers of SLARAN, such as CLARND, depend on that),
95 * the statistically correct thing to do in this situation is
96 * simply to iterate again.
97 * N.B. the case SLARAN = 0.0 should not be possible.
98 *
99 GOTO 10
100 END IF
101 *
102 SLARAN = RNDOUT
103 RETURN
104 *
105 * End of SLARAN
106 *
107 END