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.
6 Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
7 Contributed by the AriC and Caramel projects, INRIA.
9 This file is part of the GNU MPFR Library.
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.
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.
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 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
24 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
27 /* #define MPFR_NEED_LONGLONG_H */
28 #include "mpfr-impl.h"
32 mpfr_grandom (mpfr_ptr rop1, mpfr_ptr rop2, gmp_randstate_t rstate,
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;
42 if (rop2 == NULL) /* only one output requested. */
44 tprec0 = MPFR_PREC (rop1);
48 tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2));
53 /* We use "Marsaglia polar method" here (cf.
54 George Marsaglia, Normal (Gaussian) random variables for supercomputers
55 The Journal of Supercomputing, Volume 5, Number 1, 49–55
56 DOI: 10.1007/BF00155857).
58 First we draw uniform x and y in [0,1] using mpz_urandomb (in
59 fixed precision), and scale them to [-1, 1].
70 mpfr_init2 (sfr, MPFR_PREC_MIN);
71 mpfr_init2 (l, MPFR_PREC_MIN);
72 mpfr_init2 (r1, MPFR_PREC_MIN);
74 mpfr_init2 (r2, MPFR_PREC_MIN);
84 mpz_urandomb (xp, rstate, tprec);
85 mpz_urandomb (yp, rstate, tprec);
90 while (mpz_sizeinbase (s, 2) > tprec * 2); /* x^2 + y^2 <= 2^{2tprec} */
94 /* FIXME: compute s as s += 2x + 2y + 2 */
95 mpz_add_ui (a, xp, 1);
96 mpz_add_ui (b, yp, 1);
100 if ((mpz_sizeinbase (s, 2) <= 2 * tprec) ||
101 ((mpz_sizeinbase (s, 2) == 2 * tprec + 1) &&
102 (mpz_scan1 (s, 0) == 2 * tprec)))
104 /* Extend by 32 bits */
105 mpz_mul_2exp (xp, xp, 32);
106 mpz_mul_2exp (yp, yp, 32);
107 mpz_urandomb (x, rstate, 32);
108 mpz_urandomb (y, rstate, 32);
116 if (mpz_sizeinbase (s, 2) > tprec * 2)
122 /* FIXME: compute s with s -= 2x + 2y + 2 */
126 /* Compute the signs of the output */
127 mpz_urandomb (x, rstate, 2);
128 s1 = mpz_tstbit (x, 0);
129 s2 = mpz_tstbit (x, 1);
132 /* s = xp^2 + yp^2 (loop invariant) */
133 mpfr_set_prec (sfr, 2 * tprec);
134 mpfr_set_prec (l, tprec);
135 mpfr_set_z (sfr, s, MPFR_RNDN); /* exact */
136 mpfr_mul_2si (sfr, sfr, -2 * tprec, MPFR_RNDN); /* exact */
137 mpfr_log (l, sfr, MPFR_RNDN);
138 mpfr_neg (l, l, MPFR_RNDN);
139 mpfr_mul_2si (l, l, 1, MPFR_RNDN);
140 mpfr_div (l, l, sfr, MPFR_RNDN);
141 mpfr_sqrt (l, l, MPFR_RNDN);
143 mpfr_set_prec (r1, tprec);
144 mpfr_mul_z (r1, l, xp, MPFR_RNDN);
145 mpfr_div_2ui (r1, r1, tprec, MPFR_RNDN); /* exact */
147 mpfr_neg (r1, r1, MPFR_RNDN);
148 if (MPFR_CAN_ROUND (r1, tprec - 2, MPFR_PREC (rop1), rnd))
152 mpfr_set_prec (r2, tprec);
153 mpfr_mul_z (r2, l, yp, MPFR_RNDN);
154 mpfr_div_2ui (r2, r2, tprec, MPFR_RNDN); /* exact */
156 mpfr_neg (r2, r2, MPFR_RNDN);
157 if (MPFR_CAN_ROUND (r2, tprec - 2, MPFR_PREC (rop2), rnd))
163 /* Extend by 32 bits */
164 mpz_mul_2exp (xp, xp, 32);
165 mpz_mul_2exp (yp, yp, 32);
166 mpz_urandomb (x, rstate, 32);
167 mpz_urandomb (y, rstate, 32);
175 inex1 = mpfr_set (rop1, r1, rnd);
178 inex2 = mpfr_set (rop2, r2, rnd);
179 inex2 = mpfr_check_range (rop2, inex2, rnd);
181 inex1 = mpfr_check_range (rop1, inex1, rnd);
197 return INEX (inex1, inex2);