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