00001 REAL FUNCTION sexpo() 00002 C**********************************************************************C 00003 C C 00004 C C 00005 C (STANDARD-) E X P O N E N T I A L DISTRIBUTION C 00006 C C 00007 C C 00008 C**********************************************************************C 00009 C**********************************************************************C 00010 C C 00011 C FOR DETAILS SEE: C 00012 C C 00013 C AHRENS, J.H. AND DIETER, U. C 00014 C COMPUTER METHODS FOR SAMPLING FROM THE C 00015 C EXPONENTIAL AND NORMAL DISTRIBUTIONS. C 00016 C COMM. ACM, 15,10 (OCT. 1972), 873 - 882. C 00017 C C 00018 C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM C 00019 C 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C 00020 C C 00021 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C 00022 C SUNIF. The argument IR thus goes away. C 00023 C C 00024 C**********************************************************************C 00025 C 00026 C 00027 C Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N 00028 C (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION 00029 C 00030 C JJV added a Save statement for q (in Data statement) 00031 C .. Local Scalars .. 00032 REAL a,q1,u,umin,ustar 00033 INTEGER i 00034 C .. 00035 C .. Local Arrays .. 00036 REAL q(8) 00037 C .. 00038 C .. External Functions .. 00039 REAL ranf 00040 EXTERNAL ranf 00041 C .. 00042 C .. Equivalences .. 00043 EQUIVALENCE (q(1),q1) 00044 C .. 00045 C .. Save statement .. 00046 SAVE q 00047 C .. 00048 C .. Data statements .. 00049 DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833, 00050 + .9999986,.9999999/ 00051 C .. 00052 C 00053 10 a = 0.0 00054 u = ranf() 00055 GO TO 30 00056 00057 20 a = a + q1 00058 30 u = u + u 00059 C JJV changed the following to reflect the true algorithm and 00060 C JJV prevent unpredictable behavior if U is initially 0.5. 00061 C IF (u.LE.1.0) GO TO 20 00062 IF (u.LT.1.0) GO TO 20 00063 40 u = u - 1.0 00064 IF (u.GT.q1) GO TO 60 00065 50 sexpo = a + u 00066 RETURN 00067 00068 60 i = 1 00069 ustar = ranf() 00070 umin = ustar 00071 70 ustar = ranf() 00072 IF (ustar.LT.umin) umin = ustar 00073 80 i = i + 1 00074 IF (u.GT.q(i)) GO TO 70 00075 90 sexpo = a + umin*q1 00076 RETURN 00077 00078 END
1.7.1
