diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-02-05 15:19:37 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-02-05 15:19:37 +0000 |
commit | 0dc50bdf2e2d85356b094cdc940b520535a15cea (patch) | |
tree | 4c318d3e9a77b8c608d6ae6b87a3786fe4a41269 /src/root.c | |
parent | 4e0435970ebe847e7ad433ee92020ba296f33146 (diff) | |
download | mpfr-0dc50bdf2e2d85356b094cdc940b520535a15cea.tar.gz |
added a test for exact powers for mpfr_root
and fixed mpfr_root for negative x (and odd k)
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@9970 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src/root.c')
-rw-r--r-- | src/root.c | 39 |
1 files changed, 23 insertions, 16 deletions
diff --git a/src/root.c b/src/root.c index 7946c0dce..3d78fe758 100644 --- a/src/root.c +++ b/src/root.c @@ -186,22 +186,27 @@ mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) /* Compute y <- x^(1/k) using exp(log(x)/k). Assume all special cases have been eliminated before. - In the extended exponent range, overflows/underflows are not possible. */ + In the extended exponent range, overflows/underflows are not possible. + Assume x > 0, or x < 0 and k odd. +*/ static int mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) { - int inexact, exact_root = 0; + int neg, inexact, exact_root = 0; mpfr_prec_t w; /* working precision */ - mpfr_t t; + mpfr_t absx, t; MPFR_GROUP_DECL(group); MPFR_TMP_DECL(marker); MPFR_ZIV_DECL(loop); MPFR_SAVE_EXPO_DECL (expo); + + neg = MPFR_IS_NEG (x); + MPFR_TMP_INIT_ABS (absx, x); MPFR_TMP_MARK(marker); w = MPFR_PREC(y) + 10; /* Take some guard bits to prepare for the 'expt' lost bits below. - If x < 2^k, then log(x) < k, thus taking log2(k) bits should be fine. */ + If |x| < 2^k, then log|x| < k, thus taking log2(k) bits should be fine. */ if (MPFR_GET_EXP(x) > 0) w += MPFR_INT_CEIL_LOG2 (MPFR_GET_EXP(x)); MPFR_GROUP_INIT_1(group, w, t); @@ -212,30 +217,30 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) mpfr_exp_t expt; unsigned int err; - mpfr_log (t, x, MPFR_RNDN); - /* t = log(x) * (1 + theta) with |theta| <= 2^(-w) */ + 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); - /* t = log(x)/k * (1 + theta) + eps with |theta| <= 2^(-w) + /* 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) */ mpfr_exp (t, t, MPFR_RNDN); - /* t = x^(1/k) * exp(tau) * (1 + theta1) with + /* t = |x|^(1/k) * exp(tau) * (1 + theta1) with |tau| <= 1.5 * 2^(expt - w) and |theta1| <= 2^(-w). For |tau| <= 0.5 we have |exp(tau)-1| < 4/3*tau, thus for w >= expt + 2 we have: - t = x^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with + t = |x|^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with |theta1|, |theta2| <= 2^(-w). If expt+2 > 0, as long as w >= 1, we have: - t = x^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w). + t = |x|^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w). For expt+2 = 0, we have: - t = x^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w). + t = |x|^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w). Finally for expt+2 < 0 we have: - t = x^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w). + t = |x|^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w). */ err = (expt + 2 > 0) ? expt + 3 : (expt + 2 == 0) ? 2 : 1; - /* now t = x^(1/k) * (1 + 2^(err-w)) thus the error is at most + /* now t = |x|^(1/k) * (1 + 2^(err-w)) thus the error is at most 2^(EXP(t) - w + err) */ if (MPFR_LIKELY (MPFR_CAN_ROUND(t, w - err, MPFR_PREC(y), rnd_mode))) break; @@ -250,9 +255,10 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) mpfr_init2 (zk, MPFR_PREC(x)); mpfr_set (z, t, MPFR_RNDN); inexact = mpfr_pow_ui (zk, z, k, MPFR_RNDN); - exact_root = !inexact && mpfr_equal_p (zk, x); + exact_root = !inexact && mpfr_equal_p (zk, absx); if (exact_root) /* z is the exact root, thus round z directly */ - inexact = mpfr_set (y, z, rnd_mode); + inexact = (neg == 0) ? mpfr_set (y, z, rnd_mode) + : mpfr_neg (y, z, rnd_mode); mpfr_clear (zk); mpfr_clear (z); if (exact_root) @@ -265,7 +271,8 @@ mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode) MPFR_ZIV_FREE (loop); if (!exact_root) - inexact = mpfr_set (y, t, rnd_mode); + inexact = (neg == 0) ? mpfr_set (y, t, rnd_mode) + : mpfr_neg (y, t, rnd_mode); MPFR_GROUP_CLEAR(group); MPFR_TMP_FREE(marker); |