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-2004, 2006-2023 Free Software Foundation, Inc.
5Contributed by the AriC and Caramba 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
21https://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#if GMP_NUMB_BITS < 16
27#define LOGBITS_PER_BLOCK 3
28#else
29#define LOGBITS_PER_BLOCK 4
30#endif
31
32#if GMP_NUMB_BITS < 32
33#define BITS_PER_RANDCALL GMP_NUMB_BITS
34#else
35#define BITS_PER_RANDCALL 32
36#endif
37
38/* exp is the maximum exponent in absolute value */
39void
40mpfr_random2 (mpfr_ptr x, mp_size_t size, mpfr_exp_t exp,
41              gmp_randstate_t rstate)
42{
43  mp_size_t xn, k, ri;
44  unsigned long sh;
45  mp_limb_t *xp;
46  mp_limb_t elimb, ran, acc;
47  int ran_nbits, bit_pos, nb;
48
49  if (MPFR_UNLIKELY(size == 0))
50    {
51      MPFR_SET_ZERO (x);
52      MPFR_SET_POS (x);
53      return ;
54    }
55  else if (size > 0)
56    {
57      MPFR_SET_POS (x);
58    }
59  else
60    {
61      MPFR_SET_NEG (x);
62      size = -size;
63    }
64
65  xn = MPFR_LIMB_SIZE (x);
66  xp = MPFR_MANT (x);
67  if (size > xn)
68    size = xn;
69  k = xn - size;
70
71  /* Code extracted from GMP, function mpn_random2, to avoid the use
72     of GMP's internal random state in MPFR */
73
74  mpfr_rand_raw (&elimb, rstate, BITS_PER_RANDCALL);
75  ran = elimb;
76
77  /* Start off at a random bit position in the most significant limb.  */
78  bit_pos = GMP_NUMB_BITS - 1;
79  ran >>= MPFR_LOG2_GMP_NUMB_BITS;      /* Ideally   log2(GMP_NUMB_BITS) */
80  ran_nbits = BITS_PER_RANDCALL - MPFR_LOG2_GMP_NUMB_BITS; /* Ideally - log2(GMP_NUMB_BITS) */
81
82  /* Bit 0 of ran chooses string of ones/string of zeroes.
83     Make most significant limb be non-zero by setting bit 0 of RAN.  */
84  ran |= 1;
85  ri = xn - 1;
86
87  acc = 0;
88  while (ri >= k)
89    {
90      if (ran_nbits < LOGBITS_PER_BLOCK + 1)
91        {
92          mpfr_rand_raw (&elimb, rstate, BITS_PER_RANDCALL);
93          ran = elimb;
94          ran_nbits = BITS_PER_RANDCALL;
95        }
96
97      nb = (ran >> 1) % (1 << LOGBITS_PER_BLOCK) + 1;
98      if ((ran & 1) != 0)
99        {
100          MPFR_ASSERTN (bit_pos < GMP_NUMB_BITS);
101          /* Generate a string of nb ones.  */
102          if (nb > bit_pos)
103            {
104              xp[ri--] = acc | MPFR_LIMB_MASK (bit_pos + 1);
105              bit_pos += GMP_NUMB_BITS;
106              bit_pos -= nb;
107              acc = MPFR_LIMB_LSHIFT (MPFR_LIMB_MAX << 1, bit_pos);
108            }
109          else
110            {
111              bit_pos -= nb;
112              acc |= MPFR_LIMB_LSHIFT (MPFR_LIMB_MASK (nb) << 1, bit_pos);
113            }
114        }
115      else
116        {
117          /* Generate a string of nb zeroes.  */
118          if (nb > bit_pos)
119            {
120              xp[ri--] = acc;
121              acc = 0;
122              bit_pos += GMP_NUMB_BITS;
123            }
124          bit_pos -= nb;
125        }
126      ran_nbits -= LOGBITS_PER_BLOCK + 1;
127      ran >>= LOGBITS_PER_BLOCK + 1;
128    }
129
130  /* Set mandatory most significant bit.  */
131  /* xp[xn - 1] |= MPFR_LIMB_HIGHBIT; */
132
133  if (k != 0)
134    {
135      /* Clear last limbs */
136      MPN_ZERO (xp, k);
137    }
138  else
139    {
140      /* Mask off non significant bits in the low limb.  */
141      MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (x));
142      xp[0] &= ~MPFR_LIMB_MASK (sh);
143    }
144
145  /* Generate random exponent.  */
146  mpfr_rand_raw (&elimb, RANDS, GMP_NUMB_BITS);
147  MPFR_ASSERTN (exp >= 0 && exp <= MPFR_EMAX_MAX);
148  exp = (mpfr_exp_t) (elimb % (2 * exp + 1)) - exp;
149  MPFR_SET_EXP (x, exp);
150
151  return ;
152}
153