1 #include <math.h>
2 #include "grand.h"
3 #include "core_math.h"
5 double C2F(sexpo)(void)
6 /*
7 **********************************************************************
10      (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION
13 **********************************************************************
14 **********************************************************************
16 This source code was taken in the project "freemat"(BSD license)
17 This source code was modified by Gaüzère Sabine according to the
18 modifications done by JJV
20      FOR DETAILS SEE:
22                AHRENS, J.H. AND DIETER, U.
23                COMPUTER METHODS FOR SAMPLING FROM THE
24                EXPONENTIAL AND NORMAL DISTRIBUTIONS.
25                COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
27      ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
28      'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
30      Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
31      SUNIF.  The argument IR thus goes away.
33 **********************************************************************
34      Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
35      (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
36 */
37 {
38     /* ** OLD VALUES
39     static double q[8] = {
40         0.69314718055995, 0.93337368751905, 0.98887779618387, 0.99849592529150,
41         0.99982928110614, 0.99998331641007, 0.99999856914388, 0.99999989069256
42     //0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.9999999
43     };
44     */
45     static double q[8] =
46     {
47         0.69314718246459960938,
48         0.93337368965148925781,
49         0.98887777328491210938,
50         0.99849587678909301758,
51         0.99982929229736328125,
52         0.99998331069946289062,
53         0.99999862909317016602,
54         0.99999988079071044922
55     };
56     static int i;
57     static double sexpo, a, u, ustar, umin;
59     a = 0.0;
60     u = C2F(ranf)();
61     goto S30;
62 S20:
63     a += q[0];
64 S30:
65     u += u;
66     //  JJV changed the following to reflect the true algorithm and
67     //  JJV prevent unpredictable behavior if U is initially 0.5.
68     //  IF (u.LE.1.0) GO TO 20
69     if (u < 1.0)
70     {
71         goto S20;
72     }
73     u -= 1.0;
74     if (u > q[0])
75     {
76         goto S60;
77     }
78     sexpo = a + u;
79     return sexpo;
80 S60:
81     i = 1;
82     ustar = C2F(ranf)();
83     umin = ustar;
84 S70:
85     ustar = C2F(ranf)();
86     if (ustar < umin)
87     {
88         umin = ustar;
89     }
90     i += 1;
91     if (u > q[i - 1])
92     {
93         goto S70;
94     }
95     sexpo = a + umin * q[0];
96     return sexpo;
97 }