1/* mpfr_urandomb (rop, state, nbits) -- Generate a uniform pseudorandom 2 real number between 0 (inclusive) and 1 (exclusive) of size NBITS, 3 using STATE as the random state previously initialized by a call to 4 gmp_randinit_lc_2exp_size(). 5 6Copyright 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 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/* generate nbits random bits into mp[], assuming mp was allocated to contain 31 a sufficient number of limbs */ 32void 33mpfr_rand_raw (mpfr_limb_ptr mp, gmp_randstate_t rstate, 34 mpfr_prec_t nbits) 35{ 36 mpz_t z; 37 38 MPFR_ASSERTN (nbits >= 1); 39 /* To be sure to avoid the potential allocation of mpz_urandomb */ 40 ALLOC(z) = SIZ(z) = MPFR_PREC2LIMBS (nbits); 41 PTR(z) = mp; 42#if __MPFR_GMP(5,0,0) 43 /* Check for integer overflow (unless mp_bitcnt_t is signed, 44 but according to the GMP manual, this shouldn't happen). 45 Note: mp_bitcnt_t has been introduced in GMP 5.0.0. */ 46 MPFR_ASSERTN ((mp_bitcnt_t) -1 < 0 || nbits <= (mp_bitcnt_t) -1); 47#endif 48 mpz_urandomb (z, rstate, nbits); 49} 50 51int 52mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate) 53{ 54 mpfr_limb_ptr rp; 55 mpfr_prec_t nbits; 56 mp_size_t nlimbs; 57 mp_size_t k; /* number of high zero limbs */ 58 mpfr_exp_t exp; 59 int cnt; 60 61 rp = MPFR_MANT (rop); 62 nbits = MPFR_PREC (rop); 63 nlimbs = MPFR_LIMB_SIZE (rop); 64 MPFR_SET_POS (rop); 65 cnt = nlimbs * GMP_NUMB_BITS - nbits; 66 67 /* Uniform non-normalized significand */ 68 /* generate exactly nbits so that the random generator stays in the same 69 state, independent of the machine word size GMP_NUMB_BITS */ 70 mpfr_rand_raw (rp, rstate, nbits); 71 if (MPFR_LIKELY (cnt != 0)) /* this will put the low bits to zero */ 72 mpn_lshift (rp, rp, nlimbs, cnt); 73 74 /* Count the null significant limbs and remaining limbs */ 75 exp = 0; 76 k = 0; 77 while (nlimbs != 0 && rp[nlimbs - 1] == 0) 78 { 79 k ++; 80 nlimbs --; 81 exp -= GMP_NUMB_BITS; 82 } 83 84 if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */ 85 { 86 count_leading_zeros (cnt, rp[nlimbs - 1]); 87 /* Normalization */ 88 if (mpfr_set_exp (rop, exp - cnt)) 89 { 90 /* If the exponent is not in the current exponent range, we 91 choose to return a NaN as this is probably a user error. 92 Indeed this can happen only if the exponent range has been 93 reduced to a very small interval and/or the precision is 94 huge (very unlikely). */ 95 MPFR_SET_NAN (rop); 96 __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */ 97 return 1; 98 } 99 if (cnt != 0) 100 mpn_lshift (rp + k, rp, nlimbs, cnt); 101 if (k != 0) 102 MPN_ZERO (rp, k); 103 } 104 else 105 MPFR_SET_ZERO (rop); 106 107 return 0; 108} 109