GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sexpo.f
Go to the documentation of this file.
1  REAL FUNCTION sexpo()
2 C**********************************************************************C
3 C C
4 C C
5 C (STANDARD-) E X P O N E N T I A L DISTRIBUTION C
6 C C
7 C C
8 C**********************************************************************C
9 C**********************************************************************C
10 C C
11 C FOR DETAILS SEE: C
12 C C
13 C AHRENS, J.H. AND DIETER, U. C
14 C COMPUTER METHODS FOR SAMPLING FROM THE C
15 C EXPONENTIAL AND NORMAL DISTRIBUTIONS. C
16 C COMM. ACM, 15,10 (OCT. 1972), 873 - 882. C
17 C C
18 C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM C
19 C 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C
20 C C
21 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C
22 C SUNIF. The argument IR thus goes away. C
23 C C
24 C**********************************************************************C
25 C
26 C
27 C Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
28 C (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
29 C
30 C JJV added a Save statement for q (in Data statement)
31 C .. Local Scalars ..
32  REAL a,q1,u,umin,ustar
33  INTEGER i
34 C ..
35 C .. Local Arrays ..
36  REAL q(8)
37 C ..
38 C .. External Functions ..
39  REAL ranf
40  EXTERNAL ranf
41 C ..
42 C .. Equivalences ..
43  equivalence(q(1),q1)
44 C ..
45 C .. Save statement ..
46  SAVE q
47 C ..
48 C .. Data statements ..
49  DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833,
50  + .9999986,.9999999/
51 C ..
52 C
53  10 a = 0.0
54  u = ranf()
55  go to 30
56 
57  20 a = a + q1
58  30 u = u + u
59 C JJV changed the following to reflect the true algorithm and
60 C JJV prevent unpredictable behavior if U is initially 0.5.
61 C IF (u.LE.1.0) GO TO 20
62  IF (u.LT.1.0) go to 20
63  40 u = u - 1.0
64  IF (u.GT.q1) go to 60
65  50 sexpo = a + u
66  RETURN
67 
68  60 i = 1
69  ustar = ranf()
70  umin = ustar
71  70 ustar = ranf()
72  IF (ustar.LT.umin) umin = ustar
73  80 i = i + 1
74  IF (u.GT.q(i)) go to 70
75  90 sexpo = a + umin*q1
76  RETURN
77 
78  END