1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-2002 the R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/
19 *
20 * SYNOPSIS
21 *
22 * #include <Rmath.h>
23 * double exp_rand(void);
24 *
25 * DESCRIPTION
26 *
27 * Random variates from the standard exponential distribution.
28 *
29 * REFERENCE
30 *
31 * Ahrens, J.H. and Dieter, U. (1972).
32 * Computer methods for sampling from the exponential and
33 * normal distributions.
34 * Comm. ACM, 15, 873-882.
35 */
36
37 #include "nmath.h"
38
exp_rand(void)39 double exp_rand(void)
40 {
41 /* q[k-1] = sum(log(2)^k / k!) k=1,..,n, */
42 /* The highest n (here 16) is determined by q[n-1] = 1.0 */
43 /* within standard precision */
44 const static double q[] =
45 {
46 0.6931471805599453,
47 0.9333736875190459,
48 0.9888777961838675,
49 0.9984959252914960,
50 0.9998292811061389,
51 0.9999833164100727,
52 0.9999985691438767,
53 0.9999998906925558,
54 0.9999999924734159,
55 0.9999999995283275,
56 0.9999999999728814,
57 0.9999999999985598,
58 0.9999999999999289,
59 0.9999999999999968,
60 0.9999999999999999,
61 1.0000000000000000
62 };
63
64 double a = 0.;
65 double u = unif_rand(); /* precaution if u = 0 is ever returned */
66 while(u <= 0. || u >= 1.) u = unif_rand();
67 for (;;) {
68 u += u;
69 if (u > 1.)
70 break;
71 a += q[0];
72 }
73 u -= 1.;
74
75 if (u <= q[0])
76 return a + u;
77
78 int i = 0;
79 double ustar = unif_rand(), umin = ustar;
80 do {
81 ustar = unif_rand();
82 if (umin > ustar)
83 umin = ustar;
84 i++;
85 } while (u > q[i]);
86 return a + umin * q[0];
87 }
88