1 /* mpfr_grandom (rop1, rop2, state, rnd_mode) -- Generate up to two
2 pseudorandom real numbers according to a standard normal Gaussian
3 distribution and round it to the precision of rop1, rop2 according
4 to the given rounding mode.
5
6 Copyright 2011-2023 Free Software Foundation, Inc.
7 Contributed by the AriC and Caramba projects, INRIA.
8
9 This file is part of the GNU MPFR Library.
10
11 The GNU MPFR Library is free software; you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 The GNU MPFR Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
23 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
24 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
25
26
27 /* #define MPFR_NEED_LONGLONG_H */
28 #include "mpfr-impl.h"
29
30
31 int
32 mpfr_grandom (mpfr_ptr rop1, mpfr_ptr rop2, gmp_randstate_t rstate,
33 mpfr_rnd_t rnd)
34 {
35 int inex1, inex2, s1, s2;
36 mpz_t x, y, xp, yp, t, a, b, s;
37 mpfr_t sfr, l, r1, r2;
38 mpfr_prec_t tprec, tprec0;
39
40 inex2 = inex1 = 0;
41
42 if (rop2 == NULL) /* only one output requested. */
43 tprec0 = MPFR_PREC (rop1);
44 else
45 tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2));
46
47 tprec0 += 11;
48
49 /* We use "Marsaglia polar method" here (cf.
50 George Marsaglia, Normal (Gaussian) random variables for supercomputers
51 The Journal of Supercomputing, Volume 5, Number 1, 49–55
52 DOI: 10.1007/BF00155857).
53
54 First we draw uniform x and y in [0,1] using mpz_urandomb (in
55 fixed precision), and scale them to [-1, 1].
56 */
57
58 mpz_init (xp);
59 mpz_init (yp);
60 mpz_init (x);
61 mpz_init (y);
62 mpz_init (t);
63 mpz_init (s);
64 mpz_init (a);
65 mpz_init (b);
66 mpfr_init2 (sfr, MPFR_PREC_MIN);
67 mpfr_init2 (l, MPFR_PREC_MIN);
68 mpfr_init2 (r1, MPFR_PREC_MIN);
69 if (rop2 != NULL)
70 mpfr_init2 (r2, MPFR_PREC_MIN);
71
72 mpz_set_ui (xp, 0);
73 mpz_set_ui (yp, 0);
74
75 for (;;)
76 {
77 tprec = tprec0;
78 do
79 {
80 mpz_urandomb (xp, rstate, tprec);
81 mpz_urandomb (yp, rstate, tprec);
82 mpz_mul (a, xp, xp);
83 mpz_mul (b, yp, yp);
84 mpz_add (s, a, b);
85 }
86 while (mpz_sizeinbase (s, 2) > tprec * 2);
87
88 /* now s = x^2 + y^2 < 2^{2tprec} */
89
90 for (;;)
91 {
92 /* we compute (xp+1)^2 + (yp+1)^2 as s + 2xp + 2yp + 2 */
93 mpz_addmul_ui (s, xp, 2);
94 mpz_addmul_ui (s, yp, 2);
95 mpz_add_ui (s, s, 2);
96 /* The case s = 2^(2*tprec) is not possible:
97 (a) if xp and yp have different parities, s is odd
98 (b) if xp and yp are even, (xp+1)^2 and (yp+1)^2 are 1 mod 4,
99 thus s = 2 mod 4 (and tprec >= 1);
100 (c) if xp and yp are odd, if we note x = xp+1, y = yp+1 and
101 p = tprec, we would have x^2 + y^2 = 2^(2p) with x and y even
102 0 < x, y <= 2^p, thus if x' = x/2, y' = y/2 and p'=p-1,
103 we would have x'^2 + y'^2 = 2^(2p') with
104 0 < x', y' <= 2^p', and we conclude by induction. */
105 if (mpz_sizeinbase (s, 2) <= 2 * tprec)
106 goto yeepee;
107 /* Extend by 32 bits: for tprec=12, the probability we get here
108 is 8191/13180825, i.e., about 0.000621 */
109 mpz_mul_2exp (xp, xp, 32);
110 mpz_mul_2exp (yp, yp, 32);
111 mpz_urandomb (x, rstate, 32);
112 mpz_urandomb (y, rstate, 32);
113 mpz_add (xp, xp, x);
114 mpz_add (yp, yp, y);
115 tprec += 32;
116
117 mpz_mul (a, xp, xp);
118 mpz_mul (b, yp, yp);
119 mpz_add (s, a, b);
120 if (mpz_sizeinbase (s, 2) > tprec * 2)
121 break;
122 }
123 }
124 yeepee:
125
126 /* FIXME: compute s with s -= 2x + 2y + 2 */
127 mpz_mul (a, xp, xp);
128 mpz_mul (b, yp, yp);
129 mpz_add (s, a, b);
130 /* Compute the signs of the output */
131 mpz_urandomb (x, rstate, 2);
132 s1 = mpz_tstbit (x, 0);
133 s2 = mpz_tstbit (x, 1);
134 for (;;)
135 {
136 /* s = xp^2 + yp^2 (loop invariant) */
137 mpfr_set_prec (sfr, 2 * tprec);
138 mpfr_set_prec (l, tprec);
139 mpfr_set_z (sfr, s, MPFR_RNDN); /* exact */
140 mpfr_mul_2si (sfr, sfr, -2 * tprec, MPFR_RNDN); /* exact */
141 mpfr_log (l, sfr, MPFR_RNDN);
142 mpfr_neg (l, l, MPFR_RNDN);
143 mpfr_mul_2si (l, l, 1, MPFR_RNDN);
144 mpfr_div (l, l, sfr, MPFR_RNDN);
145 mpfr_sqrt (l, l, MPFR_RNDN);
146
147 mpfr_set_prec (r1, tprec);
148 mpfr_mul_z (r1, l, xp, MPFR_RNDN);
149 mpfr_div_2ui (r1, r1, tprec, MPFR_RNDN); /* exact */
150 if (s1)
151 mpfr_neg (r1, r1, MPFR_RNDN);
152 if (MPFR_CAN_ROUND (r1, tprec - 2, MPFR_PREC (rop1), rnd))
153 {
154 if (rop2 != NULL)
155 {
156 mpfr_set_prec (r2, tprec);
157 mpfr_mul_z (r2, l, yp, MPFR_RNDN);
158 mpfr_div_2ui (r2, r2, tprec, MPFR_RNDN); /* exact */
159 if (s2)
160 mpfr_neg (r2, r2, MPFR_RNDN);
161 if (MPFR_CAN_ROUND (r2, tprec - 2, MPFR_PREC (rop2), rnd))
162 break;
163 }
164 else
165 break;
166 }
167 /* Extend by 32 bits */
168 mpz_mul_2exp (xp, xp, 32);
169 mpz_mul_2exp (yp, yp, 32);
170 mpz_urandomb (x, rstate, 32);
171 mpz_urandomb (y, rstate, 32);
172 mpz_add (xp, xp, x);
173 mpz_add (yp, yp, y);
174 tprec += 32;
175 mpz_mul (a, xp, xp);
176 mpz_mul (b, yp, yp);
177 mpz_add (s, a, b);
178 }
179 inex1 = mpfr_set (rop1, r1, rnd);
180 if (rop2 != NULL)
181 {
182 inex2 = mpfr_set (rop2, r2, rnd);
183 inex2 = mpfr_check_range (rop2, inex2, rnd);
184 }
185 inex1 = mpfr_check_range (rop1, inex1, rnd);
186
187 if (rop2 != NULL)
188 mpfr_clear (r2);
189 mpfr_clear (r1);
190 mpfr_clear (l);
191 mpfr_clear (sfr);
192 mpz_clear (b);
193 mpz_clear (a);
194 mpz_clear (s);
195 mpz_clear (t);
196 mpz_clear (y);
197 mpz_clear (x);
198 mpz_clear (yp);
199 mpz_clear (xp);
200
201 return INEX (inex1, inex2);
202 }