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