1 SUBROUTINE DLARUV( ISEED, N, X )
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 DOUBLE PRECISION X( N )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLARUV returns a vector of n random real numbers from a uniform (0,1)
20 * distribution (n <= 128).
21 *
22 * This is an auxiliary routine called by DLARNV and ZLARNV.
23 *
24 * Arguments
25 * =========
26 *
27 * ISEED (input/output) INTEGER array, dimension (4)
28 * On entry, the seed of the random number generator; the array
29 * elements must be between 0 and 4095, and ISEED(4) must be
30 * odd.
31 * On exit, the seed is updated.
32 *
33 * N (input) INTEGER
34 * The number of random numbers to be generated. N <= 128.
35 *
36 * X (output) DOUBLE PRECISION array, dimension (N)
37 * The generated random numbers.
38 *
39 * Further Details
40 * ===============
41 *
42 * This routine uses a multiplicative congruential method with modulus
43 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
44 * 'Multiplicative congruential random number generators with modulus
45 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
46 * b = 48', Math. Comp. 189, pp 331-344, 1990).
47 *
48 * 48-bit integers are stored in 4 integer array elements with 12 bits
49 * per element. Hence the routine is portable across machines with
50 * integers of 32 bits or more.
51 *
52 * =====================================================================
53 *
54 * .. Parameters ..
55 DOUBLE PRECISION ONE
56 PARAMETER ( ONE = 1.0D0 )
57 INTEGER LV, IPW2
58 DOUBLE PRECISION R
59 PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J
63 * ..
64 * .. Local Arrays ..
65 INTEGER MM( LV, 4 )
66 * ..
67 * .. Intrinsic Functions ..
68 INTRINSIC DBLE, MIN, MOD
69 * ..
70 * .. Data statements ..
71 DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508,
72 $ 2549 /
73 DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754,
74 $ 1145 /
75 DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766,
76 $ 2253 /
77 DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572,
78 $ 305 /
79 DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893,
80 $ 3301 /
81 DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307,
82 $ 1065 /
83 DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297,
84 $ 3133 /
85 DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966,
86 $ 2913 /
87 DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758,
88 $ 3285 /
89 DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598,
90 $ 1241 /
91 DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406,
92 $ 1197 /
93 DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922,
94 $ 3729 /
95 DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038,
96 $ 2501 /
97 DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934,
98 $ 1673 /
99 DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091,
100 $ 541 /
101 DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451,
102 $ 2753 /
103 DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580,
104 $ 949 /
105 DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958,
106 $ 2361 /
107 DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055,
108 $ 1165 /
109 DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507,
110 $ 4081 /
111 DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078,
112 $ 2725 /
113 DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273,
114 $ 3305 /
115 DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17,
116 $ 3069 /
117 DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854,
118 $ 3617 /
119 DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916,
120 $ 3733 /
121 DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971,
122 $ 409 /
123 DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889,
124 $ 2157 /
125 DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831,
126 $ 1361 /
127 DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621,
128 $ 3973 /
129 DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541,
130 $ 1865 /
131 DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893,
132 $ 2525 /
133 DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736,
134 $ 1409 /
135 DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992,
136 $ 3445 /
137 DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787,
138 $ 3577 /
139 DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125,
140 $ 77 /
141 DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364,
142 $ 3761 /
143 DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460,
144 $ 2149 /
145 DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257,
146 $ 1449 /
147 DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574,
148 $ 3005 /
149 DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912,
150 $ 225 /
151 DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216,
152 $ 85 /
153 DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248,
154 $ 3673 /
155 DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401,
156 $ 3117 /
157 DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124,
158 $ 3089 /
159 DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762,
160 $ 1349 /
161 DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149,
162 $ 2057 /
163 DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245,
164 $ 413 /
165 DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166,
166 $ 65 /
167 DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466,
168 $ 1845 /
169 DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018,
170 $ 697 /
171 DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399,
172 $ 3085 /
173 DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190,
174 $ 3441 /
175 DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879,
176 $ 1573 /
177 DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153,
178 $ 3689 /
179 DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320,
180 $ 2941 /
181 DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18,
182 $ 929 /
183 DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712,
184 $ 533 /
185 DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159,
186 $ 2841 /
187 DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318,
188 $ 4077 /
189 DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091,
190 $ 721 /
191 DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443,
192 $ 2821 /
193 DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510,
194 $ 2249 /
195 DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449,
196 $ 2397 /
197 DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956,
198 $ 2817 /
199 DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201,
200 $ 245 /
201 DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137,
202 $ 1913 /
203 DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399,
204 $ 1997 /
205 DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321,
206 $ 3121 /
207 DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271,
208 $ 997 /
209 DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667,
210 $ 1833 /
211 DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703,
212 $ 2877 /
213 DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629,
214 $ 1633 /
215 DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365,
216 $ 981 /
217 DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431,
218 $ 2009 /
219 DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113,
220 $ 941 /
221 DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922,
222 $ 2449 /
223 DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554,
224 $ 197 /
225 DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184,
226 $ 2441 /
227 DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099,
228 $ 285 /
229 DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228,
230 $ 1473 /
231 DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012,
232 $ 2741 /
233 DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921,
234 $ 3129 /
235 DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452,
236 $ 909 /
237 DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901,
238 $ 2801 /
239 DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572,
240 $ 421 /
241 DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309,
242 $ 4073 /
243 DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171,
244 $ 2813 /
245 DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817,
246 $ 2337 /
247 DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039,
248 $ 1429 /
249 DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696,
250 $ 1177 /
251 DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256,
252 $ 1901 /
253 DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715,
254 $ 81 /
255 DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077,
256 $ 1669 /
257 DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019,
258 $ 2633 /
259 DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497,
260 $ 2269 /
261 DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101,
262 $ 129 /
263 DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717,
264 $ 1141 /
265 DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51,
266 $ 249 /
267 DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981,
268 $ 3917 /
269 DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978,
270 $ 2481 /
271 DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813,
272 $ 3941 /
273 DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881,
274 $ 2217 /
275 DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76,
276 $ 2749 /
277 DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846,
278 $ 3041 /
279 DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694,
280 $ 1877 /
281 DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682,
282 $ 345 /
283 DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124,
284 $ 2861 /
285 DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660,
286 $ 1809 /
287 DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997,
288 $ 3141 /
289 DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479,
290 $ 2825 /
291 DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141,
292 $ 157 /
293 DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886,
294 $ 2881 /
295 DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514,
296 $ 3637 /
297 DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301,
298 $ 1465 /
299 DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604,
300 $ 2829 /
301 DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888,
302 $ 2161 /
303 DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836,
304 $ 3365 /
305 DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990,
306 $ 361 /
307 DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058,
308 $ 2685 /
309 DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692,
310 $ 3745 /
311 DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194,
312 $ 2325 /
313 DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20,
314 $ 3609 /
315 DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285,
316 $ 3821 /
317 DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046,
318 $ 3537 /
319 DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107,
320 $ 517 /
321 DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508,
322 $ 3017 /
323 DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525,
324 $ 2141 /
325 DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801,
326 $ 1537 /
327 * ..
328 * .. Executable Statements ..
329 *
330 I1 = ISEED( 1 )
331 I2 = ISEED( 2 )
332 I3 = ISEED( 3 )
333 I4 = ISEED( 4 )
334 *
335 DO 10 I = 1, MIN( N, LV )
336 *
337 20 CONTINUE
338 *
339 * Multiply the seed by i-th power of the multiplier modulo 2**48
340 *
341 IT4 = I4*MM( I, 4 )
342 IT3 = IT4 / IPW2
343 IT4 = IT4 - IPW2*IT3
344 IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 )
345 IT2 = IT3 / IPW2
346 IT3 = IT3 - IPW2*IT2
347 IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 )
348 IT1 = IT2 / IPW2
349 IT2 = IT2 - IPW2*IT1
350 IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) +
351 $ I4*MM( I, 1 )
352 IT1 = MOD( IT1, IPW2 )
353 *
354 * Convert 48-bit integer to a real number in the interval (0,1)
355 *
356 X( I ) = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
357 $ DBLE( IT4 ) ) ) )
358 *
359 IF (X( I ).EQ.1.0D0) THEN
360 * If a real number has n bits of precision, and the first
361 * n bits of the 48-bit integer above happen to be all 1 (which
362 * will occur about once every 2**n calls), then X( I ) will
363 * be rounded to exactly 1.0.
364 * Since X( I ) is not supposed to return exactly 0.0 or 1.0,
365 * the statistically correct thing to do in this situation is
366 * simply to iterate again.
367 * N.B. the case X( I ) = 0.0 should not be possible.
368 I1 = I1 + 2
369 I2 = I2 + 2
370 I3 = I3 + 2
371 I4 = I4 + 2
372 GOTO 20
373 END IF
374 *
375 10 CONTINUE
376 *
377 * Return final value of seed
378 *
379 ISEED( 1 ) = IT1
380 ISEED( 2 ) = IT2
381 ISEED( 3 ) = IT3
382 ISEED( 4 ) = IT4
383 RETURN
384 *
385 * End of DLARUV
386 *
387 END
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 DOUBLE PRECISION X( N )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLARUV returns a vector of n random real numbers from a uniform (0,1)
20 * distribution (n <= 128).
21 *
22 * This is an auxiliary routine called by DLARNV and ZLARNV.
23 *
24 * Arguments
25 * =========
26 *
27 * ISEED (input/output) INTEGER array, dimension (4)
28 * On entry, the seed of the random number generator; the array
29 * elements must be between 0 and 4095, and ISEED(4) must be
30 * odd.
31 * On exit, the seed is updated.
32 *
33 * N (input) INTEGER
34 * The number of random numbers to be generated. N <= 128.
35 *
36 * X (output) DOUBLE PRECISION array, dimension (N)
37 * The generated random numbers.
38 *
39 * Further Details
40 * ===============
41 *
42 * This routine uses a multiplicative congruential method with modulus
43 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
44 * 'Multiplicative congruential random number generators with modulus
45 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
46 * b = 48', Math. Comp. 189, pp 331-344, 1990).
47 *
48 * 48-bit integers are stored in 4 integer array elements with 12 bits
49 * per element. Hence the routine is portable across machines with
50 * integers of 32 bits or more.
51 *
52 * =====================================================================
53 *
54 * .. Parameters ..
55 DOUBLE PRECISION ONE
56 PARAMETER ( ONE = 1.0D0 )
57 INTEGER LV, IPW2
58 DOUBLE PRECISION R
59 PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J
63 * ..
64 * .. Local Arrays ..
65 INTEGER MM( LV, 4 )
66 * ..
67 * .. Intrinsic Functions ..
68 INTRINSIC DBLE, MIN, MOD
69 * ..
70 * .. Data statements ..
71 DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508,
72 $ 2549 /
73 DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754,
74 $ 1145 /
75 DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766,
76 $ 2253 /
77 DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572,
78 $ 305 /
79 DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893,
80 $ 3301 /
81 DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307,
82 $ 1065 /
83 DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297,
84 $ 3133 /
85 DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966,
86 $ 2913 /
87 DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758,
88 $ 3285 /
89 DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598,
90 $ 1241 /
91 DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406,
92 $ 1197 /
93 DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922,
94 $ 3729 /
95 DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038,
96 $ 2501 /
97 DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934,
98 $ 1673 /
99 DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091,
100 $ 541 /
101 DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451,
102 $ 2753 /
103 DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580,
104 $ 949 /
105 DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958,
106 $ 2361 /
107 DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055,
108 $ 1165 /
109 DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507,
110 $ 4081 /
111 DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078,
112 $ 2725 /
113 DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273,
114 $ 3305 /
115 DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17,
116 $ 3069 /
117 DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854,
118 $ 3617 /
119 DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916,
120 $ 3733 /
121 DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971,
122 $ 409 /
123 DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889,
124 $ 2157 /
125 DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831,
126 $ 1361 /
127 DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621,
128 $ 3973 /
129 DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541,
130 $ 1865 /
131 DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893,
132 $ 2525 /
133 DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736,
134 $ 1409 /
135 DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992,
136 $ 3445 /
137 DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787,
138 $ 3577 /
139 DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125,
140 $ 77 /
141 DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364,
142 $ 3761 /
143 DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460,
144 $ 2149 /
145 DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257,
146 $ 1449 /
147 DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574,
148 $ 3005 /
149 DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912,
150 $ 225 /
151 DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216,
152 $ 85 /
153 DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248,
154 $ 3673 /
155 DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401,
156 $ 3117 /
157 DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124,
158 $ 3089 /
159 DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762,
160 $ 1349 /
161 DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149,
162 $ 2057 /
163 DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245,
164 $ 413 /
165 DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166,
166 $ 65 /
167 DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466,
168 $ 1845 /
169 DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018,
170 $ 697 /
171 DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399,
172 $ 3085 /
173 DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190,
174 $ 3441 /
175 DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879,
176 $ 1573 /
177 DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153,
178 $ 3689 /
179 DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320,
180 $ 2941 /
181 DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18,
182 $ 929 /
183 DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712,
184 $ 533 /
185 DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159,
186 $ 2841 /
187 DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318,
188 $ 4077 /
189 DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091,
190 $ 721 /
191 DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443,
192 $ 2821 /
193 DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510,
194 $ 2249 /
195 DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449,
196 $ 2397 /
197 DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956,
198 $ 2817 /
199 DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201,
200 $ 245 /
201 DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137,
202 $ 1913 /
203 DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399,
204 $ 1997 /
205 DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321,
206 $ 3121 /
207 DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271,
208 $ 997 /
209 DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667,
210 $ 1833 /
211 DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703,
212 $ 2877 /
213 DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629,
214 $ 1633 /
215 DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365,
216 $ 981 /
217 DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431,
218 $ 2009 /
219 DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113,
220 $ 941 /
221 DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922,
222 $ 2449 /
223 DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554,
224 $ 197 /
225 DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184,
226 $ 2441 /
227 DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099,
228 $ 285 /
229 DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228,
230 $ 1473 /
231 DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012,
232 $ 2741 /
233 DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921,
234 $ 3129 /
235 DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452,
236 $ 909 /
237 DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901,
238 $ 2801 /
239 DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572,
240 $ 421 /
241 DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309,
242 $ 4073 /
243 DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171,
244 $ 2813 /
245 DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817,
246 $ 2337 /
247 DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039,
248 $ 1429 /
249 DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696,
250 $ 1177 /
251 DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256,
252 $ 1901 /
253 DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715,
254 $ 81 /
255 DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077,
256 $ 1669 /
257 DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019,
258 $ 2633 /
259 DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497,
260 $ 2269 /
261 DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101,
262 $ 129 /
263 DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717,
264 $ 1141 /
265 DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51,
266 $ 249 /
267 DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981,
268 $ 3917 /
269 DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978,
270 $ 2481 /
271 DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813,
272 $ 3941 /
273 DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881,
274 $ 2217 /
275 DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76,
276 $ 2749 /
277 DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846,
278 $ 3041 /
279 DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694,
280 $ 1877 /
281 DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682,
282 $ 345 /
283 DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124,
284 $ 2861 /
285 DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660,
286 $ 1809 /
287 DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997,
288 $ 3141 /
289 DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479,
290 $ 2825 /
291 DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141,
292 $ 157 /
293 DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886,
294 $ 2881 /
295 DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514,
296 $ 3637 /
297 DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301,
298 $ 1465 /
299 DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604,
300 $ 2829 /
301 DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888,
302 $ 2161 /
303 DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836,
304 $ 3365 /
305 DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990,
306 $ 361 /
307 DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058,
308 $ 2685 /
309 DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692,
310 $ 3745 /
311 DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194,
312 $ 2325 /
313 DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20,
314 $ 3609 /
315 DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285,
316 $ 3821 /
317 DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046,
318 $ 3537 /
319 DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107,
320 $ 517 /
321 DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508,
322 $ 3017 /
323 DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525,
324 $ 2141 /
325 DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801,
326 $ 1537 /
327 * ..
328 * .. Executable Statements ..
329 *
330 I1 = ISEED( 1 )
331 I2 = ISEED( 2 )
332 I3 = ISEED( 3 )
333 I4 = ISEED( 4 )
334 *
335 DO 10 I = 1, MIN( N, LV )
336 *
337 20 CONTINUE
338 *
339 * Multiply the seed by i-th power of the multiplier modulo 2**48
340 *
341 IT4 = I4*MM( I, 4 )
342 IT3 = IT4 / IPW2
343 IT4 = IT4 - IPW2*IT3
344 IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 )
345 IT2 = IT3 / IPW2
346 IT3 = IT3 - IPW2*IT2
347 IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 )
348 IT1 = IT2 / IPW2
349 IT2 = IT2 - IPW2*IT1
350 IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) +
351 $ I4*MM( I, 1 )
352 IT1 = MOD( IT1, IPW2 )
353 *
354 * Convert 48-bit integer to a real number in the interval (0,1)
355 *
356 X( I ) = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
357 $ DBLE( IT4 ) ) ) )
358 *
359 IF (X( I ).EQ.1.0D0) THEN
360 * If a real number has n bits of precision, and the first
361 * n bits of the 48-bit integer above happen to be all 1 (which
362 * will occur about once every 2**n calls), then X( I ) will
363 * be rounded to exactly 1.0.
364 * Since X( I ) is not supposed to return exactly 0.0 or 1.0,
365 * the statistically correct thing to do in this situation is
366 * simply to iterate again.
367 * N.B. the case X( I ) = 0.0 should not be possible.
368 I1 = I1 + 2
369 I2 = I2 + 2
370 I3 = I3 + 2
371 I4 = I4 + 2
372 GOTO 20
373 END IF
374 *
375 10 CONTINUE
376 *
377 * Return final value of seed
378 *
379 ISEED( 1 ) = IT1
380 ISEED( 2 ) = IT2
381 ISEED( 3 ) = IT3
382 ISEED( 4 ) = IT4
383 RETURN
384 *
385 * End of DLARUV
386 *
387 END