d3/d00/slaran_8f_source.html

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