Commit e3e61396 authored by Vincent Lefèvre's avatar Vincent Lefèvre
Browse files

[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)

[From SVN r11982 (branches/3.1)]
parent 30630a4d
......@@ -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) */
......
......@@ -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 ();
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment