diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
commit | 6ddc203e5368cf1d4108e64c18b0bf6c7d53d176 (patch) | |
tree | dda1d9b65325699c60418cfcd02fa5d800cde1f8 /Build/source/libs/mpfr/mpfr-src/src/urandom.c | |
parent | bbfe0bab9260ac7a06a69f710b053dcc73cecad2 (diff) |
mpfr-4.0.0
git-svn-id: svn://tug.org/texlive/trunk@46140 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/urandom.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/urandom.c | 156 |
1 files changed, 58 insertions, 98 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/urandom.c b/Build/source/libs/mpfr/mpfr-src/src/urandom.c index a30f9a843dc..816e4b5b9b5 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/urandom.c +++ b/Build/source/libs/mpfr/mpfr-src/src/urandom.c @@ -27,34 +27,27 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" -/* generate one random bit */ -static int -random_rounding_bit (gmp_randstate_t rstate) -{ - mp_limb_t r; - - mpfr_rand_raw (&r, rstate, 1); - 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. +/* The mpfr_urandom() function is implemented in the following way, + so that the exact number (the random value to be rounded) and the + final status of the random generator do not depend on the current + exponent range and on the rounding mode. However, they depend on + the target precision: from the same state of the random generator, + if the precision of the destination is changed, then the value may + be completely different (and the state of the random generator is + different too). 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). + One also draws the rounding bit. This is currently done with a + separate call to mpfr_rand_raw(), but it should be better to draw + the rounding bit as part of the significand; there is space for it + since the MSB is always 1. 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 @@ -65,43 +58,28 @@ 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; + mp_limb_t rbit; int cnt; int inex; + MPFR_SAVE_EXPO_DECL (expo); + + /* We need to extend the exponent range in order to simplify + the case where one rounds upward (we would not be able to + use mpfr_nextabove() in the case emin = max). It could be + partly reimplemented under a simpler form here, but it is + better to make the code shorter and more readable. */ + MPFR_SAVE_EXPO_MARK (expo); rp = MPFR_MANT (rop); nbits = MPFR_PREC (rop); + MPFR_SET_EXP (rop, 0); MPFR_SET_POS (rop); - - 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 - || (__gmpfr_emin == 1 && rnd_mode == MPFR_RNDN - && random_rounding_bit (rstate))) - { - mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode); - MPFR_RET (+1); - } - else - { - MPFR_SET_ZERO (rop); - MPFR_RET (-1); - } - } - 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); - while (cnt == DRAW_BITS) + MPFR_STAT_STATIC_ASSERT (DRAW_BITS <= GMP_NUMB_BITS); + do { /* generate DRAW_BITS in rp[0] */ mpfr_rand_raw (rp, rstate, DRAW_BITS); @@ -112,67 +90,49 @@ 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 < __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: 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 && exp == __gmpfr_emin - 1 - && (cnt != DRAW_BITS || random_rounding_bit (rstate)))) - { - mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode); - MPFR_RET (+1); - } - else - { - MPFR_SET_ZERO (rop); - MPFR_RET (-1); - } - } - MPFR_ASSERTD (exp >= __gmpfr_emin); + /* Any value of exp < MPFR_EMIN_MIN - 1 are equivalent. So, we can + avoid a theoretical integer overflow in the following way. */ + if (MPFR_LIKELY (exp >= MPFR_EMIN_MIN - 1)) + exp -= cnt; /* no integer overflow */ } - - MPFR_ASSERTD (exp >= __gmpfr_emin); - MPFR_EXP (rop) = exp; /* Warning: may be larger than emax */ + while (cnt == DRAW_BITS); + /* We do not want the random generator to depend on the ABI or on the + exponent range. Therefore we do not use MPFR_EMIN_MIN or __gmpfr_emin + in the stop condition. */ /* 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); - rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT; - - /* Rounding */ + if (MPFR_UNLIKELY (nbits == 1)) + { + rp[0] = MPFR_LIMB_HIGHBIT; + } + else + { + 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); + rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT; + } + + /* Rounding bit */ + mpfr_rand_raw (&rbit, rstate, 1); + MPFR_ASSERTD (rbit == 0 || rbit == 1); + + /* Step 3 (rounding). */ if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA - || (rnd_mode == MPFR_RNDN && random_rounding_bit (rstate))) + || (rnd_mode == MPFR_RNDN && rbit != 0)) { - if (MPFR_UNLIKELY (exp > __gmpfr_emax)) - mpfr_set_inf (rop, +1); /* overflow */ - else - mpfr_nextabove (rop); + 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; + { + inex = -1; + } + MPFR_EXP (rop) += exp; /* may be smaller than emin */ + MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (rop, inex, rnd_mode); } |