diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-12-15 09:00:52 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-12-15 09:00:52 +0000 |
commit | b449d9b939964b18e2f45ed8dfa015a41ff0a00b (patch) | |
tree | dbc0082358086974a6b9c6a88bd118effb5e3ec1 | |
parent | c25dc97f895f808eeb2a45d363358c0f6760f339 (diff) | |
download | mpfr-b449d9b939964b18e2f45ed8dfa015a41ff0a00b.tar.gz |
[src/root.c] Fixed the case |x| = 1 and k > 100.
[tests/troot.c] Added non-regression tests.
(merged changesets r11978-11981 from the trunk)
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/3.1@11982 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | src/root.c | 16 | ||||
-rw-r--r-- | tests/troot.c | 21 |
2 files changed, 36 insertions, 1 deletions
diff --git a/src/root.c b/src/root.c index d3b5b171a..ba50b4723 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/tests/troot.c b/tests/troot.c index c6649cf8c..e173db59a 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 (); |