1/* mpfr_urandom (rop, state, rnd_mode) -- Generate a uniform pseudorandom 2 real number between 0 and 1 (exclusive) and round it to the precision of rop 3 according to the given rounding mode. 4 5Copyright 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 6Contributed by the AriC and Caramel projects, INRIA. 7 8This file is part of the GNU MPFR Library. 9 10The GNU MPFR Library is free software; you can redistribute it and/or modify 11it under the terms of the GNU Lesser General Public License as published by 12the Free Software Foundation; either version 3 of the License, or (at your 13option) any later version. 14 15The GNU MPFR Library is distributed in the hope that it will be useful, but 16WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 17or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 18License for more details. 19 20You should have received a copy of the GNU Lesser General Public License 21along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 22http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2351 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 24 25 26#define MPFR_NEED_LONGLONG_H 27#include "mpfr-impl.h" 28 29 30/* generate one random bit */ 31static int 32random_rounding_bit (gmp_randstate_t rstate) 33{ 34 mp_limb_t r; 35 36 mpfr_rand_raw (&r, rstate, 1); 37 return r & MPFR_LIMB_ONE; 38} 39 40 41int 42mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) 43{ 44 mpfr_limb_ptr rp; 45 mpfr_prec_t nbits; 46 mp_size_t nlimbs; 47 mp_size_t n; 48 mpfr_exp_t exp; 49 mpfr_exp_t emin; 50 int cnt; 51 int inex; 52 53 rp = MPFR_MANT (rop); 54 nbits = MPFR_PREC (rop); 55 nlimbs = MPFR_LIMB_SIZE (rop); 56 MPFR_SET_POS (rop); 57 exp = 0; 58 emin = mpfr_get_emin (); 59 if (MPFR_UNLIKELY (emin > 0)) 60 { 61 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA 62 || (emin == 1 && rnd_mode == MPFR_RNDN 63 && random_rounding_bit (rstate))) 64 { 65 mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); 66 return +1; 67 } 68 else 69 { 70 MPFR_SET_ZERO (rop); 71 return -1; 72 } 73 } 74 75 /* Exponent */ 76#define DRAW_BITS 8 /* we draw DRAW_BITS at a time */ 77 cnt = DRAW_BITS; 78 MPFR_ASSERTN(DRAW_BITS <= GMP_NUMB_BITS); 79 while (cnt == DRAW_BITS) 80 { 81 /* generate DRAW_BITS in rp[0] */ 82 mpfr_rand_raw (rp, rstate, DRAW_BITS); 83 if (MPFR_UNLIKELY (rp[0] == 0)) 84 cnt = DRAW_BITS; 85 else 86 { 87 count_leading_zeros (cnt, rp[0]); 88 cnt -= GMP_NUMB_BITS - DRAW_BITS; 89 } 90 91 if (MPFR_UNLIKELY (exp < emin + cnt)) 92 { 93 /* To get here, we have been drawing more than -emin zeros 94 in a row, then return 0 or the smallest representable 95 positive number. 96 97 The rounding to nearest mode is subtle: 98 If exp - cnt == emin - 1, the rounding bit is set, except 99 if cnt == DRAW_BITS in which case the rounding bit is 100 outside rp[0] and must be generated. */ 101 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA 102 || (rnd_mode == MPFR_RNDN && cnt == exp - emin - 1 103 && (cnt != DRAW_BITS || random_rounding_bit (rstate)))) 104 { 105 mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); 106 return +1; 107 } 108 else 109 { 110 MPFR_SET_ZERO (rop); 111 return -1; 112 } 113 } 114 exp -= cnt; 115 } 116 MPFR_EXP (rop) = exp; /* Warning: may be outside the current 117 exponent range */ 118 119 120 /* Significand: we need generate only nbits-1 bits, since the most 121 significant is 1 */ 122 mpfr_rand_raw (rp, rstate, nbits - 1); 123 n = nlimbs * GMP_NUMB_BITS - nbits; 124 if (MPFR_LIKELY (n != 0)) /* this will put the low bits to zero */ 125 mpn_lshift (rp, rp, nlimbs, n); 126 127 /* Set the msb to 1 since it was fixed by the exponent choice */ 128 rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT; 129 130 /* Rounding */ 131 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA 132 || (rnd_mode == MPFR_RNDN && random_rounding_bit (rstate))) 133 { 134 /* Take care of the exponent range: it may have been reduced */ 135 if (exp < emin) 136 mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); 137 else if (exp > mpfr_get_emax ()) 138 mpfr_set_inf (rop, +1); /* overflow, flag set by mpfr_check_range */ 139 else 140 mpfr_nextabove (rop); 141 inex = +1; 142 } 143 else 144 inex = -1; 145 146 return mpfr_check_range (rop, inex, rnd_mode); 147} 148