1/* mpfr_random2 -- Generate a positive random mpfr_t of specified size, with
2   long runs of consecutive ones and zeros in the binary representation.
3
4Copyright 1999, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
5Contributed by the Arenaire and Cacao projects, INRIA.
6
7This file is part of the GNU MPFR Library.
8
9The GNU MPFR Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MPFR Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24#include "mpfr-test.h"
25
26#define LOGBITS_PER_BLOCK 4
27#if GMP_NUMB_BITS < 32
28#define BITS_PER_RANDCALL GMP_NUMB_BITS
29#else
30#define BITS_PER_RANDCALL 32
31#endif
32
33void
34mpfr_random2 (mpfr_ptr x, mp_size_t size, mpfr_exp_t exp,
35              gmp_randstate_t rstate)
36{
37  mp_size_t xn, k, ri;
38  unsigned long sh;
39  mp_limb_t *xp;
40  mp_limb_t elimb, ran, acc;
41  int ran_nbits, bit_pos, nb;
42
43  if (MPFR_UNLIKELY(size == 0))
44    {
45      MPFR_SET_ZERO (x);
46      MPFR_SET_POS (x);
47      return ;
48    }
49  else if (size > 0)
50    {
51      MPFR_SET_POS (x);
52    }
53  else
54    {
55      MPFR_SET_NEG (x);
56      size = -size;
57    }
58
59  xn = MPFR_LIMB_SIZE (x);
60  xp = MPFR_MANT (x);
61  if (size > xn)
62    size = xn;
63  k = xn - size;
64
65  /* Code extracted from GMP, function mpn_random2, to avoid the use
66     of GMP's internal random state in MPFR */
67
68  mpfr_rand_raw (&elimb, rstate, BITS_PER_RANDCALL);
69  ran = elimb;
70
71  /* Start off at a random bit position in the most significant limb.  */
72  bit_pos = GMP_NUMB_BITS - 1;
73  ran >>= 6;                            /* Ideally   log2(GMP_NUMB_BITS) */
74  ran_nbits = BITS_PER_RANDCALL - 6;    /* Ideally - log2(GMP_NUMB_BITS) */
75
76  /* Bit 0 of ran chooses string of ones/string of zeroes.
77     Make most significant limb be non-zero by setting bit 0 of RAN.  */
78  ran |= 1;
79  ri = xn - 1;
80
81  acc = 0;
82  while (ri >= k)
83    {
84      if (ran_nbits < LOGBITS_PER_BLOCK + 1)
85        {
86          mpfr_rand_raw (&elimb, rstate, BITS_PER_RANDCALL);
87          ran = elimb;
88          ran_nbits = BITS_PER_RANDCALL;
89        }
90
91      nb = (ran >> 1) % (1 << LOGBITS_PER_BLOCK) + 1;
92      if ((ran & 1) != 0)
93        {
94          /* Generate a string of nb ones.  */
95          if (nb > bit_pos)
96            {
97              xp[ri--] = acc | (((mp_limb_t) 2 << bit_pos) - 1);
98              bit_pos += GMP_NUMB_BITS;
99              bit_pos -= nb;
100              acc = ((~(mp_limb_t) 1) << bit_pos) & GMP_NUMB_MASK;
101            }
102          else
103            {
104              bit_pos -= nb;
105              acc |= (((mp_limb_t) 2 << nb) - 2) << bit_pos;
106            }
107        }
108      else
109        {
110          /* Generate a string of nb zeroes.  */
111          if (nb > bit_pos)
112            {
113              xp[ri--] = acc;
114              acc = 0;
115              bit_pos += GMP_NUMB_BITS;
116            }
117          bit_pos -= nb;
118        }
119      ran_nbits -= LOGBITS_PER_BLOCK + 1;
120      ran >>= LOGBITS_PER_BLOCK + 1;
121    }
122
123  /* Set mandatory most significant bit.  */
124  /* xp[xn - 1] |= MPFR_LIMB_HIGHBIT; */
125
126  if (k != 0)
127    {
128      /* Clear last limbs */
129      MPN_ZERO (xp, k);
130    }
131  else
132    {
133      /* Mask off non significant bits in the low limb.  */
134      MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (x));
135      xp[0] &= ~MPFR_LIMB_MASK (sh);
136    }
137
138  /* Generate random exponent.  */
139  mpfr_rand_raw (&elimb, RANDS, GMP_NUMB_BITS);
140  exp = ABS (exp);
141  MPFR_SET_EXP (x, elimb % (2 * exp + 1) - exp);
142
143  return ;
144}
145