12 #if !defined(RANDOMLIB_MPFRRANDOM_HPP)
13 #define RANDOMLIB_MPFRRANDOM_HPP 1
18 #define HAVE_MPFR (MPFR_VERSION_MAJOR >= 3)
20 #if HAVE_MPFR || defined(DOXYGEN)
25 #if !defined(STATIC_ASSERT)
26 # if defined(__GXX_EXPERIMENTAL_CXX0X__)
27 # define STATIC_ASSERT static_assert
28 # elif defined(_MSC_VER) && _MSC_VER >= 1600
29 # define STATIC_ASSERT static_assert
31 # define STATIC_ASSERT(cond,reason) \
32 { enum{ STATIC_ASSERT_ENUM = 1/int(cond) }; }
54 static const int limb_ = GMP_LIMB_BITS;
55 static const int loglimb_ = (limb_ == 32 ? 5 :
57 (limb_ == 128 ? 7 : -1)));
58 static const int logbits_ = (bits == 1 ? 0 :
65 (bits == 128 ? 7 : -1))))))));
66 static const mp_limb_t mask_ = (bits == limb_ ? ~0UL :
67 ~(~0UL << (bits < limb_ ? bits : 0)));
68 static const int logw_ = loglimb_ - logbits_;
69 static const unsigned wmask_ = ~(~0U << logw_);
76 void AddDigits(gmp_randstate_t r,
long num = 1) {
78 mpz_mul_2exp(_f, _f, num << logbits_);
79 mpz_urandomb(_tt, r, num << logbits_);
84 mp_limb_t Digit(gmp_randstate_t r, mp_size_t k) {
90 (mpz_getlimbn(_f, k >> logw_) >> ((k & wmask_) << logbits_));
94 static int highest_bit_idx(
unsigned long x)
throw() {
98 if (x & 0xffff0000UL) { x >>= 16; r += 16; }
99 if (x & 0x0000ff00UL) { x >>= 8; r += 8; }
100 if (x & 0x000000f0UL) { x >>= 4; r += 4; }
101 if (x & 0x0000000cUL) { x >>= 2; r += 2; }
102 if (x & 0x00000002UL) { r += 1; }
110 STATIC_ASSERT(logbits_ >= 0 && loglimb_ >= 0 && logbits_ <= loglimb_,
111 "MPRFRandom: unsupported value for bits");
112 mpz_init(_f); mpz_init(_tt);
120 { mpz_init(_f); mpz_set(_f, t._f); mpz_init(_tt); }
154 void Init() { mpz_set_ui(_f, 0u); _e = 0; _n = 0; _s = 1; }
158 int Sign()
const throw() {
return _s; }
166 long Floor()
const throw() {
return _s > 0 ? long(_n) : -1 - long(_n); }
170 long Ceiling()
const throw() {
return _s > 0 ? 1 + long(_n) : -long(_n); }
174 unsigned long UInteger()
const throw() {
return _n; }
178 unsigned long Size()
const throw() {
return unsigned(_e); }
186 int ns = k < 0 ? -1 : 1;
189 mpz_mul_2exp(_tt, _tt, _e << logbits_);
190 mpz_sub_ui(_tt, _tt, 1u);
191 mpz_sub(_f, _tt, _f);
194 _n = ns > 0 ? k : -(k + 1);
204 if (
this == &t)
return false;
205 if (_s != t._s)
return _s < t._s;
206 if (_n != t._n)
return (_s < 0) ^ (_n < t._n);
207 for (mp_size_t k = 0; ; ++k) {
208 mp_limb_t x = Digit(r, k);
209 mp_limb_t y = t.Digit(r, k);
210 if (x != y)
return (_s < 0) ^ (x < y);
220 mpz_setbit(_f, (_e << logbits_) - 1);
229 return mpz_tstbit(_f, (_e << logbits_) - 1);
240 if (_n)
return highest_bit_idx(_n);
242 int sgn = mpz_sgn(_f);
244 return mp_size_t(mpz_sizeinbase(_f, 2)) - mp_size_t(_e << logbits_);
256 {
if (_e <= k) AddDigits(r, k - _e + 1); }
284 int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) {
293 round = _s < 0 ? MPFR_RNDU : MPFR_RNDD;
296 round = _s < 0 ? MPFR_RNDD : MPFR_RNDU;
310 expt = -(_e << logbits_);
313 if (round == MPFR_RNDN) {
322 mpfr_prec_t prec = mpfr_get_prec (val);
323 mp_size_t trail = lead - prec;
324 mp_size_t guard = trail + (round == MPFR_RNDN ? 0 : 1);
326 if (guard <= 0)
ExpandTo(r, (-guard) >> logbits_);
341 excess = (_e << logbits_) + expt;
344 mpz_mul_2exp(_tt, _tt, _e << logbits_);
345 mpz_add(_tt, _tt, _f);
347 mpz_tdiv_q_2exp(_tt, _tt, excess);
349 mpz_mul_2exp(_tt, _tt, -excess);
350 if (r || round == MPFR_RNDN)
353 else if (round == MPFR_RNDU)
355 mpz_add_ui(_tt, _tt, 1u);
360 int flag = mpfr_set_z_2exp(val, _tt, expt, round);
362 mpfr_neg (val, val, MPFR_RNDN);
375 mpz_urandomb(_tt, r, 1);
376 return mpz_tstbit(_tt, 0);
383 #endif // RANDOMLIB_MPFRRANDOM_HPP
void ExpandTo(gmp_randstate_t r, mp_size_t k)
unsigned long Size() const
int TestHighBit(gmp_randstate_t r)
MPFRRandom & operator=(const MPFRRandom &t)
unsigned long UInteger() const
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round)
#define STATIC_ASSERT(cond, reason)
int LessThan(gmp_randstate_t r, MPFRRandom &t)
void swap(RandomLib::RandomCanonical< Generator > &r, RandomLib::RandomCanonical< Generator > &s)
MPFRRandom(const MPFRRandom &t)
int operator()(mpfr_t val, mpfr_rnd_t round)
void SetHighBit(gmp_randstate_t r)
int Boolean(gmp_randstate_t r) const
Handling random numbers in MPFR.
mp_size_t LeadingBit(gmp_randstate_t r)