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