diff --git a/PATCHES b/PATCHES index aebe097..d5aa55b 100644 --- a/PATCHES +++ b/PATCHES @@ -1 +1,2 @@ +root mpfr_get diff --git a/VERSION b/VERSION index 9fa2663..7a62f99 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.1.6-p1 +3.1.6-p2 diff --git a/src/mpfr.h b/src/mpfr.h index 7c15508..ed718be 100644 --- a/src/mpfr.h +++ b/src/mpfr.h @@ -27,7 +27,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_VERSION_MAJOR 3 #define MPFR_VERSION_MINOR 1 #define MPFR_VERSION_PATCHLEVEL 6 -#define MPFR_VERSION_STRING "3.1.6-p1" +#define MPFR_VERSION_STRING "3.1.6-p2" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff --git a/src/root.c b/src/root.c index d3b5b17..ba50b47 100644 --- a/src/root.c +++ b/src/root.c @@ -107,6 +107,11 @@ mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) MPFR_RET_NAN; } + /* Special case |x| = 1. Note that if x = -1, then k is odd + (NaN results have already been filtered), so that y = -1. */ + if (mpfr_cmpabs (x, __gmpfr_one) == 0) + return mpfr_set (y, x, rnd_mode); + /* General case */ /* For large k, use exp(log(x)/k). The threshold of 100 seems to be quite @@ -188,6 +193,14 @@ mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) Assume all special cases have been eliminated before. In the extended exponent range, overflows/underflows are not possible. Assume x > 0, or x < 0 and k odd. + Also assume |x| <> 1 because log(1) = 0, which does not have an exponent + and would yield a failure in the error bound computation. A priori, this + constraint is quite artificial because if |x| is close enough to 1, then + the exponent of log|x| does not need to be used (in the code, err would + be 1 in such a domain). So this constraint |x| <> 1 could be avoided in + the code. However, this is an exact case easy to detect, so that such a + change would be useless. Values very close to 1 are not an issue, since + an underflow is not possible before the MPFR_GET_EXP. */ static int mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) @@ -219,7 +232,8 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) mpfr_log (t, absx, MPFR_RNDN); /* t = log|x| * (1 + theta) with |theta| <= 2^(-w) */ mpfr_div_ui (t, t, k, MPFR_RNDN); - expt = MPFR_GET_EXP (t); + /* No possible underflow in mpfr_log and mpfr_div_ui. */ + expt = MPFR_GET_EXP (t); /* assumes t <> 0 */ /* t = log|x|/k * (1 + theta) + eps with |theta| <= 2^(-w) and |eps| <= 1/2 ulp(t), thus the total error is bounded by 1.5 * 2^(expt - w) */ diff --git a/src/version.c b/src/version.c index 84fdc4a..3b05c73 100644 --- a/src/version.c +++ b/src/version.c @@ -25,5 +25,5 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., const char * mpfr_get_version (void) { - return "3.1.6-p1"; + return "3.1.6-p2"; } diff --git a/tests/troot.c b/tests/troot.c index c6649cf..e173db5 100644 --- a/tests/troot.c +++ b/tests/troot.c @@ -405,6 +405,26 @@ cmp_pow (void) mpfr_clears (x, y1, y2, (mpfr_ptr) 0); } +static void +bug20171214 (void) +{ + mpfr_t x, y; + int inex; + + mpfr_init2 (x, 805); + mpfr_init2 (y, 837); + mpfr_set_ui (x, 1, MPFR_RNDN); + inex = mpfr_root (y, x, 120, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + MPFR_ASSERTN (mpfr_cmp_ui (y, 1) == 0); + mpfr_set_si (x, -1, MPFR_RNDN); + inex = mpfr_root (y, x, 121, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + MPFR_ASSERTN (mpfr_cmp_si (y, -1) == 0); + mpfr_clear (x); + mpfr_clear (y); +} + int main (void) { @@ -415,6 +435,7 @@ main (void) tests_start_mpfr (); + bug20171214 (); exact_powers (3, 1000); special (); bigint ();