Revert "Merge remote branch 'origin/sas' into YaSp"
[scilab.git] / scilab / modules / randlib / src / c / sexpo.c
1 #include "stack-c.h"
2 #include <math.h>
3 #include "grand.h"
4 #include "core_math.h"
5
6 double C2F(sexpo)(void)
7 /*
8 **********************************************************************
9
10
11      (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION
12
13
14 **********************************************************************
15 **********************************************************************
16
17 This source code was taken in the project "freemat"(BSD license)
18 This source code was modified by Gaüzère Sabine according to the 
19 modifications done by JJV
20
21      FOR DETAILS SEE:
22
23                AHRENS, J.H. AND DIETER, U.
24                COMPUTER METHODS FOR SAMPLING FROM THE
25                EXPONENTIAL AND NORMAL DISTRIBUTIONS.
26                COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
27
28      ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
29      'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
30
31      Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
32      SUNIF.  The argument IR thus goes away.
33
34 **********************************************************************
35      Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
36      (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
37 */
38 {
39 /* ** OLD VALUES
40 static double q[8] = {
41     0.69314718055995, 0.93337368751905, 0.98887779618387, 0.99849592529150,
42     0.99982928110614, 0.99998331641007, 0.99999856914388, 0.99999989069256
43 //0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.9999999
44 };
45 */
46 static double q[8] = {
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;
58
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) goto S20;
70     u -= 1.0;
71     if(u > q[0]) goto S60;
72     sexpo = a+u;
73     return sexpo;
74 S60:
75     i = 1;
76     ustar = C2F(ranf)();
77     umin = ustar;
78 S70:
79     ustar = C2F(ranf)();
80     if(ustar < umin) umin = ustar;
81     i += 1;
82     if(u > q[i-1]) goto S70;
83     sexpo = a+umin*q[0];
84     return sexpo;
85 }