GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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*4 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
ranf
real function ranf()
Definition: ranf.f:2
sexpo
real function sexpo()
Definition: sexpo.f:2