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