LAPACK 3.3.0
|
00001 REAL FUNCTION SLARAN( ISEED ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Array Arguments .. 00008 INTEGER ISEED( 4 ) 00009 * .. 00010 * 00011 * Purpose 00012 * ======= 00013 * 00014 * SLARAN returns a random real number from a uniform (0,1) 00015 * distribution. 00016 * 00017 * Arguments 00018 * ========= 00019 * 00020 * ISEED (input/output) INTEGER array, dimension (4) 00021 * On entry, the seed of the random number generator; the array 00022 * elements must be between 0 and 4095, and ISEED(4) must be 00023 * odd. 00024 * On exit, the seed is updated. 00025 * 00026 * Further Details 00027 * =============== 00028 * 00029 * This routine uses a multiplicative congruential method with modulus 00030 * 2**48 and multiplier 33952834046453 (see G.S.Fishman, 00031 * 'Multiplicative congruential random number generators with modulus 00032 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for 00033 * b = 48', Math. Comp. 189, pp 331-344, 1990). 00034 * 00035 * 48-bit integers are stored in 4 integer array elements with 12 bits 00036 * per element. Hence the routine is portable across machines with 00037 * integers of 32 bits or more. 00038 * 00039 * ===================================================================== 00040 * 00041 * .. Parameters .. 00042 INTEGER M1, M2, M3, M4 00043 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 ) 00044 REAL ONE 00045 PARAMETER ( ONE = 1.0E+0 ) 00046 INTEGER IPW2 00047 REAL R 00048 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 ) 00049 * .. 00050 * .. Local Scalars .. 00051 INTEGER IT1, IT2, IT3, IT4 00052 REAL RNDOUT 00053 * .. 00054 * .. Intrinsic Functions .. 00055 INTRINSIC MOD, REAL 00056 * .. 00057 * .. Executable Statements .. 00058 10 CONTINUE 00059 * 00060 * multiply the seed by the multiplier modulo 2**48 00061 * 00062 IT4 = ISEED( 4 )*M4 00063 IT3 = IT4 / IPW2 00064 IT4 = IT4 - IPW2*IT3 00065 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3 00066 IT2 = IT3 / IPW2 00067 IT3 = IT3 - IPW2*IT2 00068 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2 00069 IT1 = IT2 / IPW2 00070 IT2 = IT2 - IPW2*IT1 00071 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 + 00072 $ ISEED( 4 )*M1 00073 IT1 = MOD( IT1, IPW2 ) 00074 * 00075 * return updated seed 00076 * 00077 ISEED( 1 ) = IT1 00078 ISEED( 2 ) = IT2 00079 ISEED( 3 ) = IT3 00080 ISEED( 4 ) = IT4 00081 * 00082 * convert 48-bit integer to a real number in the interval (0,1) 00083 * 00084 RNDOUT = R*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R* 00085 $ ( REAL( IT4 ) ) ) ) ) 00086 * 00087 IF (RNDOUT.EQ.1.0) THEN 00088 * If a real number has n bits of precision, and the first 00089 * n bits of the 48-bit integer above happen to be all 1 (which 00090 * will occur about once every 2**n calls), then SLARAN will 00091 * be rounded to exactly 1.0. In IEEE single precision arithmetic, 00092 * this will happen relatively often since n = 24. 00093 * Since SLARAN is not supposed to return exactly 0.0 or 1.0 00094 * (and some callers of SLARAN, such as CLARND, depend on that), 00095 * the statistically correct thing to do in this situation is 00096 * simply to iterate again. 00097 * N.B. the case SLARAN = 0.0 should not be possible. 00098 * 00099 GOTO 10 00100 END IF 00101 * 00102 SLARAN = RNDOUT 00103 RETURN 00104 * 00105 * End of SLARAN 00106 * 00107 END