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 6Copyright 2011, 2012, 2013 Free Software Foundation, Inc. 7Contributed by the AriC and Caramel projects, INRIA. 8 9This file is part of the GNU MPFR Library. 10 11The GNU MPFR Library is free software; you can redistribute it and/or modify 12it under the terms of the GNU Lesser General Public License as published by 13the Free Software Foundation; either version 3 of the License, or (at your 14option) any later version. 15 16The GNU MPFR Library is distributed in the hope that it will be useful, but 17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 19License for more details. 20 21You should have received a copy of the GNU Lesser General Public License 22along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 23http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2451 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 25 26 27/* #define MPFR_NEED_LONGLONG_H */ 28#include "mpfr-impl.h" 29 30 31int 32mpfr_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