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, 2012, 2013 Free Software Foundation, Inc. 7 Contributed by the AriC and Caramel 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 http://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 { 44 tprec0 = MPFR_PREC (rop1); 45 } 46 else 47 { 48 tprec0 = MAX (MPFR_PREC (rop1), MPFR_PREC (rop2)); 49 } 50 51 tprec0 += 11; 52 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). 57 58 First we draw uniform x and y in [0,1] using mpz_urandomb (in 59 fixed precision), and scale them to [-1, 1]. 60 */ 61 62 mpz_init (xp); 63 mpz_init (yp); 64 mpz_init (x); 65 mpz_init (y); 66 mpz_init (t); 67 mpz_init (s); 68 mpz_init (a); 69 mpz_init (b); 70 mpfr_init2 (sfr, MPFR_PREC_MIN); 71 mpfr_init2 (l, MPFR_PREC_MIN); 72 mpfr_init2 (r1, MPFR_PREC_MIN); 73 if (rop2 != NULL) 74 mpfr_init2 (r2, MPFR_PREC_MIN); 75 76 mpz_set_ui (xp, 0); 77 mpz_set_ui (yp, 0); 78 79 for (;;) 80 { 81 tprec = tprec0; 82 do 83 { 84 mpz_urandomb (xp, rstate, tprec); 85 mpz_urandomb (yp, rstate, tprec); 86 mpz_mul (a, xp, xp); 87 mpz_mul (b, yp, yp); 88 mpz_add (s, a, b); 89 } 90 while (mpz_sizeinbase (s, 2) > tprec * 2); /* x^2 + y^2 <= 2^{2tprec} */ 91 92 for (;;) 93 { 94 /* FIXME: compute s as s += 2x + 2y + 2 */ 95 mpz_add_ui (a, xp, 1); 96 mpz_add_ui (b, yp, 1); 97 mpz_mul (a, a, a); 98 mpz_mul (b, b, b); 99 mpz_add (s, a, b); 100 if ((mpz_sizeinbase (s, 2) <= 2 * tprec) || 101 ((mpz_sizeinbase (s, 2) == 2 * tprec + 1) && 102 (mpz_scan1 (s, 0) == 2 * tprec))) 103 goto yeepee; 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); 109 mpz_add (xp, xp, x); 110 mpz_add (yp, yp, y); 111 tprec += 32; 112 113 mpz_mul (a, xp, xp); 114 mpz_mul (b, yp, yp); 115 mpz_add (s, a, b); 116 if (mpz_sizeinbase (s, 2) > tprec * 2) 117 break; 118 } 119 } 120 yeepee: 121 122 /* FIXME: compute s with s -= 2x + 2y + 2 */ 123 mpz_mul (a, xp, xp); 124 mpz_mul (b, yp, yp); 125 mpz_add (s, a, b); 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); 130 for (;;) 131 { 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); 142 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 */ 146 if (s1) 147 mpfr_neg (r1, r1, MPFR_RNDN); 148 if (MPFR_CAN_ROUND (r1, tprec - 2, MPFR_PREC (rop1), rnd)) 149 { 150 if (rop2 != NULL) 151 { 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 */ 155 if (s2) 156 mpfr_neg (r2, r2, MPFR_RNDN); 157 if (MPFR_CAN_ROUND (r2, tprec - 2, MPFR_PREC (rop2), rnd)) 158 break; 159 } 160 else 161 break; 162 } 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); 168 mpz_add (xp, xp, x); 169 mpz_add (yp, yp, y); 170 tprec += 32; 171 mpz_mul (a, xp, xp); 172 mpz_mul (b, yp, yp); 173 mpz_add (s, a, b); 174 } 175 inex1 = mpfr_set (rop1, r1, rnd); 176 if (rop2 != NULL) 177 { 178 inex2 = mpfr_set (rop2, r2, rnd); 179 inex2 = mpfr_check_range (rop2, inex2, rnd); 180 } 181 inex1 = mpfr_check_range (rop1, inex1, rnd); 182 183 if (rop2 != NULL) 184 mpfr_clear (r2); 185 mpfr_clear (r1); 186 mpfr_clear (l); 187 mpfr_clear (sfr); 188 mpz_clear (b); 189 mpz_clear (a); 190 mpz_clear (s); 191 mpz_clear (t); 192 mpz_clear (y); 193 mpz_clear (x); 194 mpz_clear (yp); 195 mpz_clear (xp); 196 197 return INEX (inex1, inex2); 198 } 199