1 /* mpfr_erandom (rop, state, rnd_mode) -- Generate an exponential deviate with
2 mean 1 and round it to the precision of rop according to the given rounding
3 mode.
4
5 Copyright 2013-2023 Free Software Foundation, Inc.
6 Contributed by Charles Karney <charles@karney.com>, SRI International.
7
8 This file is part of the GNU MPFR Library.
9
10 The GNU MPFR Library is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
14
15 The GNU MPFR Library is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 License for more details.
19
20 You should have received a copy of the GNU Lesser General Public License
21 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
22 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
23 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24
25 /*
26 * Sampling from the exponential distribution with unit mean using the method
27 * given in John von Neumann, Various techniques used in connection with random
28 * digits, in A. S. Householder, G. E. Forsythe, and H. H. Germond, editors,
29 * "Monte Carlo Method", number 12 in Applied Mathematics Series, pp. 36-38
30 * (NBS, Washington, DC, 1951), Proceedings of a symposium held June 29-July 1,
31 * 1949, in Los Angeles.
32 *
33 * A modification to this algorithm is given in:
34 * Charles F. F. Karney,
35 * "Sampling exactly from the normal distribution",
36 * ACM Trans. Math. Software 42(1), 3:1-14 (Jan. 2016).
37 * https://dx.doi.org/10.1145/2710016
38 * https://arxiv.org/abs/1303.6257
39 * Although this improves the bit efficiency, in practice, it results in
40 * a slightly slower algorithm for MPFR. So here the original von Neumann
41 * algorithm is used.
42 *
43 * There are a few "weasel words" regarding the accuracy of this
44 * implementation. The algorithm produces exactly rounded exponential deviates
45 * provided that gmp's random number engine delivers truly random bits. If it
46 * did, the algorithm would be perfect; however, this implementation would have
47 * problems, e.g., in that the integer part of the exponential deviate is
48 * represented by an unsigned long, whereas in reality the integer part in
49 * unbounded. In this implementation, asserts catch overflow in the integer
50 * part and similar (very, very) unlikely events. In reality, of course, gmp's
51 * random number engine has a finite internal state (19937 bits in the case of
52 * the MT19937 method). This means that these unlikely events in fact won't
53 * occur. If the asserts are triggered, then this is an indication that the
54 * random number engine is defective. (Even if a hardware random number
55 * generator were used, the most likely explanation for the triggering of the
56 * asserts would be that the hardware generator was broken.)
57 */
58
59 #include "random_deviate.h"
60
61 /* true with prob exp(-x) */
62 static int
E(mpfr_random_deviate_t x,gmp_randstate_t r,mpfr_random_deviate_t p,mpfr_random_deviate_t q)63 E (mpfr_random_deviate_t x, gmp_randstate_t r,
64 mpfr_random_deviate_t p, mpfr_random_deviate_t q)
65 {
66 /* p and q are temporaries */
67 mpfr_random_deviate_reset (p);
68 if (!mpfr_random_deviate_less (p, x, r))
69 return 1;
70 for (;;)
71 {
72 mpfr_random_deviate_reset (q);
73 if (!mpfr_random_deviate_less (q, p, r))
74 return 0;
75 mpfr_random_deviate_reset (p);
76 if (!mpfr_random_deviate_less (p, q, r))
77 return 1;
78 }
79 }
80
81 /* return an exponential random deviate with mean 1 as a MPFR */
82 int
mpfr_erandom(mpfr_ptr z,gmp_randstate_t r,mpfr_rnd_t rnd)83 mpfr_erandom (mpfr_ptr z, gmp_randstate_t r, mpfr_rnd_t rnd)
84 {
85 mpfr_random_deviate_t x, p, q;
86 int inex;
87 unsigned long k = 0;
88
89 mpfr_random_deviate_init (x);
90 mpfr_random_deviate_init (p);
91 mpfr_random_deviate_init (q);
92 while (!E(x, r, p, q))
93 {
94 ++k;
95 /* Catch k wrapping around to 0; for a 32-bit unsigned long, the
96 * probability of this is exp(-2^32)). */
97 MPFR_ASSERTN (k != 0UL);
98 mpfr_random_deviate_reset (x);
99 }
100 mpfr_random_deviate_clear (q);
101 mpfr_random_deviate_clear (p);
102 inex = mpfr_random_deviate_value (0, k, x, z, r, rnd);
103 mpfr_random_deviate_clear (x);
104 return inex;
105 }
106