1 /* mpfr_nrandom (rop, state, rnd_mode) -- Generate a normal deviate with mean 0
2 and variance 1 and round it to the precision of rop according to the given
3 rounding 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 normal distribution with zero mean and unit variance.
27 * This uses Algorithm N given in:
28 * Charles F. F. Karney,
29 * "Sampling exactly from the normal distribution",
30 * ACM Trans. Math. Software 42(1), 3:1-14 (Jan. 2016).
31 * https://dx.doi.org/10.1145/2710016
32 * https://arxiv.org/abs/1303.6257
33 *
34 * Note: the algorithm implemented here has been improved in
35 * Du, Fan and Wei in "An improved exact sampling algorithm for the standard
36 * normal distribution", Computational Statistics, 2021,
37 * https://doi.org/10.1007/s00180-021-01136-w
38 *
39 * The implementation here closely follows the C++ one given in the paper
40 * above. However, here, C is simplified by using gmp_urandomm_ui; the initial
41 * rejection step in H just tests the leading bit of p; and the assignment of
42 * the sign to the deviate using gmp_urandomb_ui. Finally, the C++
43 * implementation benefits from caching temporary random deviates across calls.
44 * This isn't possible in C without additional machinery which would complicate
45 * the interface.
46 *
47 * There are a few "weasel words" regarding the accuracy of this
48 * implementation. The algorithm produces exactly rounded normal deviates
49 * provided that gmp's random number engine delivers truly random bits. If it
50 * did, the algorithm would be perfect; however, this implementation would have
51 * problems, e.g., in that the integer part of the normal deviate is
52 * represented by an unsigned long, whereas in reality the integer part in
53 * unbounded. In this implementation, asserts catch overflow in the integer
54 * part and similar (very, very) unlikely events. In reality, of course, gmp's
55 * random number engine has a finite internal state (19937 bits in the case of
56 * the MT19937 method). This means that these unlikely events in fact won't
57 * occur. If the asserts are triggered, then this is an indication that the
58 * random number engine is defective. (Even if a hardware random number
59 * generator were used, the most likely explanation for the triggering of the
60 * asserts would be that the hardware generator was broken.)
61 */
62
63 #include "random_deviate.h"
64
65 /* Algorithm H: true with probability exp(-1/2). */
66 static int
H(gmp_randstate_t r,mpfr_random_deviate_t p,mpfr_random_deviate_t q)67 H (gmp_randstate_t r, mpfr_random_deviate_t p, mpfr_random_deviate_t q)
68 {
69 /* p and q are temporaries */
70 mpfr_random_deviate_reset (p);
71 if (mpfr_random_deviate_tstbit (p, 1, r))
72 return 1;
73 for (;;)
74 {
75 mpfr_random_deviate_reset (q);
76 if (!mpfr_random_deviate_less (q, p, r))
77 return 0;
78 mpfr_random_deviate_reset (p);
79 if (!mpfr_random_deviate_less (p, q, r))
80 return 1;
81 }
82 }
83
84 /* Step N1: return n >= 0 with prob. exp(-n/2) * (1 - exp(-1/2)). */
85 static unsigned long
G(gmp_randstate_t r,mpfr_random_deviate_t p,mpfr_random_deviate_t q)86 G (gmp_randstate_t r, mpfr_random_deviate_t p, mpfr_random_deviate_t q)
87 {
88 /* p and q are temporaries */
89 unsigned long n = 0;
90
91 while (H (r, p, q))
92 {
93 ++n;
94 /* Catch n wrapping around to 0; for a 32-bit unsigned long, the
95 * probability of this is exp(-2^30)). */
96 MPFR_ASSERTN (n != 0UL);
97 }
98 return n;
99 }
100
101 /* Step N2: true with probability exp(-m*n/2). */
102 static int
P(unsigned long m,unsigned long n,gmp_randstate_t r,mpfr_random_deviate_t p,mpfr_random_deviate_t q)103 P (unsigned long m, unsigned long n, gmp_randstate_t r,
104 mpfr_random_deviate_t p, mpfr_random_deviate_t q)
105 {
106 /* p and q are temporaries. m*n is passed as two separate parameters to deal
107 * with the case where m*n overflows an unsigned long. This may be called
108 * with m = 0 and n = (unsigned long)(-1) and, because m in handled in to the
109 * outer loop, this routine will correctly return 1. */
110 while (m--)
111 {
112 unsigned long k = n;
113 while (k--)
114 {
115 if (!H (r, p, q))
116 return 0;
117 }
118 }
119 return 1;
120 }
121
122 /* Algorithm C: return (-1, 0, 1) with prob (1/m, 1/m, 1-2/m). */
123 static int
C(unsigned long m,gmp_randstate_t r)124 C (unsigned long m, gmp_randstate_t r)
125 {
126 unsigned long n = gmp_urandomm_ui (r, m);
127 return n == 0 ? -1 : (n == 1 ? 0 : 1);
128 }
129
130 /* Algorithm B: true with prob exp(-x * (2*k + x) / (2*k + 2)). */
131 static int
B(unsigned long k,mpfr_random_deviate_t x,gmp_randstate_t r,mpfr_random_deviate_t p,mpfr_random_deviate_t q)132 B (unsigned long k, mpfr_random_deviate_t x, gmp_randstate_t r,
133 mpfr_random_deviate_t p, mpfr_random_deviate_t q)
134 {
135 /* p and q are temporaries */
136
137 unsigned long m = 2 * k + 2;
138 /* n tracks the parity of the loop; s == 1 on first trip through loop. */
139 unsigned n = 0, s = 1;
140 int f;
141
142 /* Check if 2 * k + 2 would overflow; for a 32-bit unsigned long, the
143 * probability of this is exp(-2^61)). */
144 MPFR_ASSERTN (k < ((unsigned long)(-1) >> 1));
145
146 for (;; ++n, s = 0) /* overflow of n is innocuous */
147 {
148 if ( ((f = k ? 0 : C (m, r)) < 0) ||
149 (mpfr_random_deviate_reset (q),
150 !mpfr_random_deviate_less (q, s ? x : p, r)) ||
151 ((f = k ? C (m, r) : f) < 0) ||
152 (f == 0 &&
153 (mpfr_random_deviate_reset (p),
154 !mpfr_random_deviate_less (p, x, r))) )
155 break;
156 mpfr_random_deviate_swap (p, q); /* an efficient way of doing p = q */
157 }
158 return (n & 1U) == 0;
159 }
160
161 /* return a normal random deviate with mean 0 and variance 1 as a MPFR */
162 int
mpfr_nrandom(mpfr_ptr z,gmp_randstate_t r,mpfr_rnd_t rnd)163 mpfr_nrandom (mpfr_ptr z, gmp_randstate_t r, mpfr_rnd_t rnd)
164 {
165 mpfr_random_deviate_t x, p, q;
166 int inex;
167 unsigned long k, j;
168
169 mpfr_random_deviate_init (x);
170 mpfr_random_deviate_init (p);
171 mpfr_random_deviate_init (q);
172 for (;;)
173 {
174 k = G (r, p, q); /* step 1 */
175 if (!P (k, k - 1, r, p, q))
176 continue; /* step 2 */
177 mpfr_random_deviate_reset (x); /* step 3 */
178 for (j = 0; j <= k && B (k, x, r, p, q); ++j); /* step 4 */
179 if (j > k)
180 break;
181 }
182 mpfr_random_deviate_clear (q);
183 mpfr_random_deviate_clear (p);
184 /* steps 5, 6, 7 */
185 inex = mpfr_random_deviate_value (gmp_urandomb_ui (r, 1), k, x, z, r, rnd);
186 mpfr_random_deviate_clear (x);
187 return inex;
188 }
189