diff options
| -rw-r--r-- | doc/mpfr.texi | 9 | ||||
| -rw-r--r-- | src/urandom.c | 77 | ||||
| -rw-r--r-- | tests/turandom.c | 162 |
3 files changed, 220 insertions, 28 deletions
diff --git a/doc/mpfr.texi b/doc/mpfr.texi index 8de7cad68..167fb6f13 100644 --- a/doc/mpfr.texi +++ b/doc/mpfr.texi @@ -2740,6 +2740,15 @@ created using the GMP @code{gmp_randinit} function (see the GMP manual). Note: the note for @code{mpfr_urandomb} holds too. In addition, the exponent range and the rounding mode might have a side effect on the next random state. + +The rule for the underflow flag is here ``underflow before rounding'' +instead of the usual ``underflow after rounding''. The reason is that +the exponent is drawn first, and if it is smaller than the minimum +exponent, the significand is not drawn. To fix the behavior on the +underflow flag, one would have to draw the significand in some cases, +meaning that the behavior of the random generator would change, thus +it would break the ABI for the MPFR 3.1 branch. However, the observed +behavior always corresponds to an existing number. @end deftypefun @deftypefun int mpfr_grandom (mpfr_t @var{rop1}, mpfr_t @var{rop2}, gmp_randstate_t @var{state}, mpfr_rnd_t @var{rnd}) diff --git a/src/urandom.c b/src/urandom.c index 0145be52c..0283bdccb 100644 --- a/src/urandom.c +++ b/src/urandom.c @@ -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); @@ -89,7 +114,7 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) } exp -= cnt; /* no integer overflow */ - if (MPFR_UNLIKELY (exp < emin)) + 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 @@ -99,46 +124,46 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode) If exp == 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. */ + __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW; if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA - || (rnd_mode == MPFR_RNDN && 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); } } + 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; diff --git a/tests/turandom.c b/tests/turandom.c index ec5ff8ff9..bedd98ca4 100644 --- a/tests/turandom.c +++ b/tests/turandom.c @@ -38,6 +38,7 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, long count = 0; int i; int inex = 1; + unsigned int ex_flags, flags; size_tab = (nbtests >= 1000 ? nbtests / 50 : 20); tab = (int *) calloc (size_tab, sizeof(int)); @@ -60,7 +61,10 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, for (k = 0; k < nbtests; k++) { + mpfr_clear_flags (); + ex_flags = MPFR_FLAGS_INEXACT; i = mpfr_urandom (x, RANDS, rnd); + flags = __gmpfr_flags; inex = (i != 0) && inex; /* check that lower bits are zero */ if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh) && !MPFR_IS_ZERO (x)) @@ -76,6 +80,17 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, mpfr_print_binary (x); puts (""); exit (1); } + /* check the flags (an underflow is theoretically possible, but + impossible in practice due to the huge exponent range) */ + if (flags != ex_flags) + { + printf ("Error: mpfr_urandom() returns incorrect flags.\n"); + printf ("Expected "); + flags_out (ex_flags); + printf ("Got "); + flags_out (flags); + exit (1); + } d = mpfr_get_d1 (x); av += d; var += d*d; i = (int)(size_tab * d); @@ -98,7 +113,20 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, for (k = 0; k < 5; k++) { set_emin (k+1); + mpfr_clear_flags (); + ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; inex = mpfr_urandom (x, RANDS, rnd); + flags = __gmpfr_flags; + if (flags != ex_flags) + { + printf ("Error: mpfr_urandom() returns incorrect flags" + " for emin = %d.\n", k+1); + printf ("Expected "); + flags_out (ex_flags); + printf ("Got "); + flags_out (flags); + exit (1); + } if (( (rnd == MPFR_RNDZ || rnd == MPFR_RNDD) && (!MPFR_IS_ZERO (x) || inex != -1)) || ((rnd == MPFR_RNDU || rnd == MPFR_RNDA) @@ -130,7 +158,7 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, printf ("Average = %.5f\nVariance = %.5f\n", av, var); printf ("Repartition for urandom with rounding mode %s. " "Each integer should be close to %d.\n", - mpfr_print_rnd_mode (rnd), (int)th); + mpfr_print_rnd_mode (rnd), (int) th); for (k = 0; k < size_tab; k++) { @@ -209,6 +237,134 @@ bug20170123 (void) mpfr_set_emin (emin); } +static void +underflow_tests (void) +{ + mpfr_t x; + mpfr_exp_t emin; + int i, k; + int inex; + int rnd; + unsigned int ex_flags, flags; + + emin = mpfr_get_emin (); + mpfr_init2 (x, 4); + ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; /* if underflow */ + for (i = 2; i >= -4; i--) + { + mpfr_set_emin (i); + RND_LOOP (rnd) + for (k = 0; k < 100; k++) + { + mpfr_clear_flags (); + inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd); + flags = __gmpfr_flags; + MPFR_ASSERTN (mpfr_inexflag_p ()); + if (MPFR_IS_NEG (x)) + { + printf ("Error in underflow_tests: got a negative sign" + " for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + if (MPFR_IS_ZERO (x)) + { + if (rnd == MPFR_RNDU || rnd == MPFR_RNDA) + { + printf ("Error in underflow_tests: the value cannot" + " be 0 for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + if (flags != ex_flags) + { + printf ("Error in underflow_tests: incorrect flags" + " for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + printf ("Expected "); + flags_out (ex_flags); + printf ("Got "); + flags_out (flags); + exit (1); + } + } + if (inex == 0 || (MPFR_IS_ZERO (x) && inex > 0)) + { + printf ("Error in underflow_tests: incorrect inex (%d)" + " for i=%d rnd=%s k=%d.\n", inex, + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + } + } + mpfr_clear (x); + mpfr_set_emin (emin); +} + +static void +overflow_tests (void) +{ + mpfr_t x; + mpfr_exp_t emax; + int i, k; + int inex; + int rnd; + unsigned int ex_flags, flags; + + emax = mpfr_get_emax (); + mpfr_init2 (x, 4); + ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; /* if overflow */ + for (i = -4; i <= 0; i++) + { + mpfr_set_emax (i); + RND_LOOP (rnd) + for (k = 0; k < 100; k++) + { + mpfr_clear_flags (); + inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd); + flags = __gmpfr_flags; + MPFR_ASSERTN (mpfr_inexflag_p ()); + if (MPFR_IS_NEG (x)) + { + printf ("Error in overflow_tests: got a negative sign" + " for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + if (MPFR_IS_INF (x)) + { + if (rnd == MPFR_RNDD || rnd == MPFR_RNDZ) + { + printf ("Error in overflow_tests: the value cannot" + " be +inf for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + if (flags != ex_flags) + { + printf ("Error in overflow_tests: incorrect flags" + " for i=%d rnd=%s k=%d.\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + printf ("Expected "); + flags_out (ex_flags); + printf ("Got "); + flags_out (flags); + exit (1); + } + } + if (inex == 0 || (MPFR_IS_INF (x) && inex < 0)) + { + printf ("Error in overflow_tests: incorrect inex (%d)" + " for i=%d rnd=%s k=%d.\n", inex, + i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k); + exit (1); + } + } + } + mpfr_clear (x); + mpfr_set_emax (emax); +} + int main (int argc, char *argv[]) { @@ -260,8 +416,10 @@ main (int argc, char *argv[]) } } - bug20100914 (); + underflow_tests (); + overflow_tests (); + bug20100914 (); bug20170123 (); tests_end_mpfr (); |