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().
6 Copyright 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
7 Contributed by the Arenaire and Cacao 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"
31 mpfr_rand_raw (mp_ptr mp, gmp_randstate_t rstate, unsigned long int nbits)
35 /* To be sure to avoid the potential allocation of mpz_urandomb */
36 ALLOC(z) = SIZ(z) = (nbits / GMP_NUMB_BITS) + 1;
38 mpz_urandomb(z, rstate, nbits);
42 mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
47 mp_size_t k; /* number of high zero limbs */
52 nbits = MPFR_PREC (rop);
53 nlimbs = MPFR_LIMB_SIZE (rop);
56 /* Uniform non-normalized significand */
57 mpfr_rand_raw (rp, rstate, nlimbs * GMP_NUMB_BITS);
59 /* If nbits isn't a multiple of GMP_NUMB_BITS, mask the low bits */
60 cnt = nlimbs * GMP_NUMB_BITS - nbits;
61 if (MPFR_LIKELY (cnt != 0))
62 rp[0] &= ~MPFR_LIMB_MASK (cnt);
64 /* Count the null significant limbs and remaining limbs */
67 while (nlimbs != 0 && rp[nlimbs - 1] == 0)
74 if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
76 count_leading_zeros (cnt, rp[nlimbs - 1]);
78 if (mpfr_set_exp (rop, exp - cnt))
80 /* If the exponent is not in the current exponent range, we
81 choose to return a NaN as this is probably a user error.
82 Indeed this can happen only if the exponent range has been
83 reduced to a very small interval and/or the precision is
84 huge (very unlikely). */
86 __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
90 mpn_lshift (rp + k, rp, nlimbs, cnt);