diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/urandom.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/urandom.c | 95 |
1 files changed, 63 insertions, 32 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/urandom.c b/Build/source/libs/mpfr/mpfr-src/src/urandom.c index 1232ff18aa3..a30f9a843dc 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/urandom.c +++ b/Build/source/libs/mpfr/mpfr-src/src/urandom.c @@ -2,7 +2,7 @@ real number between 0 and 1 (exclusive) and round it to the precision of rop according to the given rounding mode. -Copyright 2000-2004, 2006-2016 Free Software Foundation, Inc. +Copyright 2000-2004, 2006-2017 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -37,6 +37,25 @@ random_rounding_bit (gmp_randstate_t rstate) return r & MPFR_LIMB_ONE; } +/* NOTE: The current behavior is to consider "underflow before rounding" + (the significand does not need to be drawn), while the rule in MPFR + is "underflow after rounding". This is unfixable in this 3.1 branch + without changing the behavior of the PRNG (thus breaking the ABI). */ + +/* The mpfr_urandom() function is implemented in the following way for + the generic case. + 1. One determines the exponent exp: 0 with probability 1/2, -1 with + probability 1/4, -2 with probability 1/8, etc. + 2. One draws a 1-ulp interval ]a,b[ containing the exact result (the + interval can be regarded as open since it has the same measure as + the closed interval). + 3. Rounding is done. For the directed rounding modes, the rounded value + is uniquely determined. For rounding to nearest, ]a,m[ and ]m,b[, + where m = (a+b)/2, have the same measure, so that one gets a or b + with equal probabilities. + Note: Only low-level functions are used (except just before a "return"), + so that we do not need MPFR_SAVE_EXPO_*. +*/ int mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) @@ -46,33 +65,39 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) mp_size_t nlimbs; mp_size_t n; mpfr_exp_t exp; - mpfr_exp_t emin; int cnt; int inex; rp = MPFR_MANT (rop); nbits = MPFR_PREC (rop); - nlimbs = MPFR_LIMB_SIZE (rop); MPFR_SET_POS (rop); - exp = 0; - emin = mpfr_get_emin (); - if (MPFR_UNLIKELY (emin > 0)) + + if (MPFR_UNLIKELY (__gmpfr_emin > 0)) { + /* The minimum positive representable number 2^(emin-1) is >= 1, + so that we need to round to +0 or 2^(emin-1). For the directed + rounding modes, the rounded value is uniquely determined. For + rounding to nearest: if emin = 1, one has probability 1/2 for + each; otherwise (i.e. if emin > 1), the rounded value is 0. */ + __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW; if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA - || (emin == 1 && rnd_mode == MPFR_RNDN + || (__gmpfr_emin == 1 && rnd_mode == MPFR_RNDN && random_rounding_bit (rstate))) { - mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); - return +1; + mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode); + MPFR_RET (+1); } else { MPFR_SET_ZERO (rop); - return -1; + MPFR_RET (-1); } } - /* Exponent */ + exp = 0; + MPFR_ASSERTD (exp >= __gmpfr_emin); + + /* Step 1 (exponent). */ #define DRAW_BITS 8 /* we draw DRAW_BITS at a time */ cnt = DRAW_BITS; MPFR_ASSERTN(DRAW_BITS <= GMP_NUMB_BITS); @@ -87,58 +112,64 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) count_leading_zeros (cnt, rp[0]); cnt -= GMP_NUMB_BITS - DRAW_BITS; } + exp -= cnt; /* no integer overflow */ - if (MPFR_UNLIKELY (exp < emin + cnt)) + if (MPFR_UNLIKELY (exp < __gmpfr_emin)) { /* To get here, we have been drawing more than -emin zeros in a row, then return 0 or the smallest representable positive number. - The rounding to nearest mode is subtle: - If exp - cnt == emin - 1, the rounding bit is set, except - if cnt == DRAW_BITS in which case the rounding bit is - outside rp[0] and must be generated. */ + The rounding-to-nearest mode is subtle: We need to round to + the smallest representable positive number iff the exponent + is emin - 1. This condition can be satisfied only if the + current emin is emin - 1. In this case, if cnt != DRAW_BITS, + this in the final emin, so that the condition is satisfied. + But if cnt == DRAW_BITS, we need to draw an additional bit + to determine whether emin == emin - 1 or emin < emin - 1 + (with equal probabilities); the reason is that we return + just below instead of doing more iterations in the "while" + loop to find the final value of emin. */ + __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW; if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA - || (rnd_mode == MPFR_RNDN && cnt == exp - emin - 1 + || (rnd_mode == MPFR_RNDN && exp == __gmpfr_emin - 1 && (cnt != DRAW_BITS || random_rounding_bit (rstate)))) { - mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); - return +1; + mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode); + MPFR_RET (+1); } else { MPFR_SET_ZERO (rop); - return -1; + MPFR_RET (-1); } } - exp -= cnt; + MPFR_ASSERTD (exp >= __gmpfr_emin); } - MPFR_EXP (rop) = exp; /* Warning: may be outside the current - exponent range */ + MPFR_ASSERTD (exp >= __gmpfr_emin); + MPFR_EXP (rop) = exp; /* Warning: may be larger than emax */ - /* Significand: we need generate only nbits-1 bits, since the most - significant is 1 */ + /* Step 2 (significand): we need generate only nbits-1 bits, since the + most significant bit is 1. */ mpfr_rand_raw (rp, rstate, nbits - 1); + nlimbs = MPFR_LIMB_SIZE (rop); n = nlimbs * GMP_NUMB_BITS - nbits; if (MPFR_LIKELY (n != 0)) /* this will put the low bits to zero */ mpn_lshift (rp, rp, nlimbs, n); - - /* Set the msb to 1 since it was fixed by the exponent choice */ rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT; /* Rounding */ if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA || (rnd_mode == MPFR_RNDN && random_rounding_bit (rstate))) { - /* Take care of the exponent range: it may have been reduced */ - if (exp < emin) - mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode); - else if (exp > mpfr_get_emax ()) - mpfr_set_inf (rop, +1); /* overflow, flag set by mpfr_check_range */ + if (MPFR_UNLIKELY (exp > __gmpfr_emax)) + mpfr_set_inf (rop, +1); /* overflow */ else mpfr_nextabove (rop); inex = +1; + /* There is an overflow in the first case and possibly in the second + case. If this occurs, the flag will be set by mpfr_check_range. */ } else inex = -1; |