00001 /* dlaran.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 doublereal dlaran_(integer *iseed) 00017 { 00018 /* System generated locals */ 00019 doublereal ret_val; 00020 00021 /* Local variables */ 00022 integer it1, it2, it3, it4; 00023 doublereal rndout; 00024 00025 00026 /* -- LAPACK auxiliary routine (version 3.1) -- */ 00027 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00028 /* November 2006 */ 00029 00030 /* .. Array Arguments .. */ 00031 /* .. */ 00032 00033 /* Purpose */ 00034 /* ======= */ 00035 00036 /* DLARAN returns a random real number from a uniform (0,1) */ 00037 /* distribution. */ 00038 00039 /* Arguments */ 00040 /* ========= */ 00041 00042 /* ISEED (input/output) INTEGER array, dimension (4) */ 00043 /* On entry, the seed of the random number generator; the array */ 00044 /* elements must be between 0 and 4095, and ISEED(4) must be */ 00045 /* odd. */ 00046 /* On exit, the seed is updated. */ 00047 00048 /* Further Details */ 00049 /* =============== */ 00050 00051 /* This routine uses a multiplicative congruential method with modulus */ 00052 /* 2**48 and multiplier 33952834046453 (see G.S.Fishman, */ 00053 /* 'Multiplicative congruential random number generators with modulus */ 00054 /* 2**b: an exhaustive analysis for b = 32 and a partial analysis for */ 00055 /* b = 48', Math. Comp. 189, pp 331-344, 1990). */ 00056 00057 /* 48-bit integers are stored in 4 integer array elements with 12 bits */ 00058 /* per element. Hence the routine is portable across machines with */ 00059 /* integers of 32 bits or more. */ 00060 00061 /* ===================================================================== */ 00062 00063 /* .. Parameters .. */ 00064 /* .. */ 00065 /* .. Local Scalars .. */ 00066 /* .. */ 00067 /* .. Intrinsic Functions .. */ 00068 /* .. */ 00069 /* .. Executable Statements .. */ 00070 /* Parameter adjustments */ 00071 --iseed; 00072 00073 /* Function Body */ 00074 L10: 00075 00076 /* multiply the seed by the multiplier modulo 2**48 */ 00077 00078 it4 = iseed[4] * 2549; 00079 it3 = it4 / 4096; 00080 it4 -= it3 << 12; 00081 it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508; 00082 it2 = it3 / 4096; 00083 it3 -= it2 << 12; 00084 it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322; 00085 it1 = it2 / 4096; 00086 it2 -= it1 << 12; 00087 it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4] 00088 * 494; 00089 it1 %= 4096; 00090 00091 /* return updated seed */ 00092 00093 iseed[1] = it1; 00094 iseed[2] = it2; 00095 iseed[3] = it3; 00096 iseed[4] = it4; 00097 00098 /* convert 48-bit integer to a real number in the interval (0,1) */ 00099 00100 rndout = ((doublereal) it1 + ((doublereal) it2 + ((doublereal) it3 + ( 00101 doublereal) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4) 00102 * 2.44140625e-4; 00103 00104 if (rndout == 1.) { 00105 /* If a real number has n bits of precision, and the first */ 00106 /* n bits of the 48-bit integer above happen to be all 1 (which */ 00107 /* will occur about once every 2**n calls), then DLARAN will */ 00108 /* be rounded to exactly 1.0. */ 00109 /* Since DLARAN is not supposed to return exactly 0.0 or 1.0 */ 00110 /* (and some callers of DLARAN, such as CLARND, depend on that), */ 00111 /* the statistically correct thing to do in this situation is */ 00112 /* simply to iterate again. */ 00113 /* N.B. the case DLARAN = 0.0 should not be possible. */ 00114 00115 goto L10; 00116 } 00117 00118 ret_val = rndout; 00119 return ret_val; 00120 00121 /* End of DLARAN */ 00122 00123 } /* dlaran_ */