diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2008-08-11 08:09:14 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2008-08-11 08:09:14 +0000 |
commit | 985b630fe5ed82efc6376652cecd82552031d8f3 (patch) | |
tree | 10d2442c33ad9867c038dbc25015894c8e059489 /pow_ui.c | |
parent | 6851316430ab3b8c175175ad8d4c179d051d7c8f (diff) | |
download | mpfr-985b630fe5ed82efc6376652cecd82552031d8f3.tar.gz |
Merged vlefevre branch:
svn merge -c-5445 .
svn merge -r5436:HEAD .../mpfr/branches/vlefevre
* pow.c:
- Moved the general case from mpfr_pow() to a new internal function
mpfr_pow_general().
- In this function (from old code), avoid unnecessary overflow test
if the intermediate result is not an infinity (which was the case
of underflow with non-zero result, thus not an overflow).
- Fixed a double-rounding problem that occurred in this function in
some underflow cases when rescaling the result.
- Added log messages.
* mpfr-impl.h: added mpfr_pow_general prototype.
* pow_z.c:
- The underflow case of mpfr_pow_pos_z() in rounding to nearest,
which was incorrect, is now handled by calling mpfr_pow_general(),
which can scale the result thus decide whether the rounded result
should be 0 or nextabove(0). To avoid the exact cases of x^y with
y integer (not supported by mpfr_pow_general()), rounding is done
in precision 2 (this is also faster!).
- Fixed underflow-related bug (case exact result = 2^(emin-2), in
rounding to nearest).
- Added log messages.
* pow_ui.c:
- Swapped parameters x and y for consistency (-> y = x^n).
- Fixed the internal overflows and underflows (which could yield
spurious overflows/underflows and incorrect results) by using
mpfr_pow_z.
* tests/tpow_all.c:
- Test flags in test_others and cmpres; cmpres argument z1 can now
be a null pointer (if unknown pure FP value, thus not tested).
- Added a test of 2^(emin - i/4) with 0 <= i <= 12, that triggered
the bugs mentioned above (and now fixed).
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@5505 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'pow_ui.c')
-rw-r--r-- | pow_ui.c | 94 |
1 files changed, 50 insertions, 44 deletions
@@ -24,9 +24,9 @@ MA 02110-1301, USA. */ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -/* sets x to y^n, and return 0 if exact, non-zero otherwise */ +/* sets y to x^n, and return 0 if exact, non-zero otherwise */ int -mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) +mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mp_rnd_t rnd) { unsigned long m; mpfr_t res; @@ -37,58 +37,61 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) MPFR_ZIV_DECL (loop); MPFR_BLOCK_DECL (flags); - /* y^0 = 1 for any y, even a NaN */ + MPFR_LOG_FUNC (("x[%#R]=%R n=%lu rnd=%d", x, x, n, rnd), + ("y[%#R]=%R inexact=%d", y, y, inexact)); + + /* x^0 = 1 for any x, even a NaN */ if (MPFR_UNLIKELY (n == 0)) - return mpfr_set_ui (x, 1, rnd); + return mpfr_set_ui (y, 1, rnd); - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (y))) + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) { - if (MPFR_IS_NAN (y)) + if (MPFR_IS_NAN (x)) { - MPFR_SET_NAN (x); + MPFR_SET_NAN (y); MPFR_RET_NAN; } - else if (MPFR_IS_INF (y)) + else if (MPFR_IS_INF (x)) { /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ - if ((MPFR_IS_NEG (y)) && ((n & 1) == 1)) - MPFR_SET_NEG (x); + if (MPFR_IS_NEG (x) && (n & 1) == 1) + MPFR_SET_NEG (y); else - MPFR_SET_POS (x); - MPFR_SET_INF (x); + MPFR_SET_POS (y); + MPFR_SET_INF (y); MPFR_RET (0); } - else /* y is zero */ + else /* x is zero */ { - MPFR_ASSERTD (MPFR_IS_ZERO (y)); + MPFR_ASSERTD (MPFR_IS_ZERO (x)); /* 0^n = 0 for any n */ - MPFR_SET_ZERO (x); - if (MPFR_IS_POS (y) || ((n & 1) == 0)) - MPFR_SET_POS (x); + MPFR_SET_ZERO (y); + if (MPFR_IS_POS (x) || (n & 1) == 0) + MPFR_SET_POS (y); else - MPFR_SET_NEG (x); + MPFR_SET_NEG (y); MPFR_RET (0); } } else if (MPFR_UNLIKELY (n <= 2)) { if (n < 2) - /* y^1 = y */ - return mpfr_set (x, y, rnd); + /* x^1 = x */ + return mpfr_set (y, x, rnd); else - /* y^2 = sqr(y) */ - return mpfr_sqr (x, y, rnd); + /* x^2 = sqr(x) */ + return mpfr_sqr (y, x, rnd); } /* Augment exponent range */ MPFR_SAVE_EXPO_MARK (expo); /* setup initial precision */ - prec = MPFR_PREC (x) + 3 + BITS_PER_MP_LIMB - + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); + prec = MPFR_PREC (y) + 3 + BITS_PER_MP_LIMB + + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)); mpfr_init2 (res, prec); - rnd1 = MPFR_IS_POS (y) ? GMP_RNDU : GMP_RNDD; /* away */ + rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */ MPFR_ZIV_INIT (loop, prec); for (;;) @@ -100,17 +103,17 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) /* now 2^(i-1) <= n < 2^i */ MPFR_ASSERTD (prec > (mpfr_prec_t) i); err = prec - 1 - (mpfr_prec_t) i; - /* First step: compute square from y */ + /* First step: compute square from x */ MPFR_BLOCK (flags, - inexact = mpfr_mul (res, y, y, GMP_RNDU); + inexact = mpfr_mul (res, x, x, GMP_RNDU); MPFR_ASSERTD (i >= 2); if (n & (1UL << (i-2))) - inexact |= mpfr_mul (res, res, y, rnd1); + inexact |= mpfr_mul (res, res, x, rnd1); for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--) { inexact |= mpfr_mul (res, res, res, GMP_RNDU); if (n & (1UL << i)) - inexact |= mpfr_mul (res, res, y, rnd1); + inexact |= mpfr_mul (res, res, x, rnd1); }); /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2, and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1. @@ -122,7 +125,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) */ if (MPFR_LIKELY (inexact == 0 || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags) - || MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd))) + || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd))) break; /* Actualisation of the precision */ MPFR_ZIV_NEXT (loop, prec); @@ -130,26 +133,29 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) } MPFR_ZIV_FREE (loop); - /* Check Overflow */ - if (MPFR_OVERFLOW (flags)) - { - mpfr_clear (res); - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_overflow (x, rnd, - (n % 2) ? MPFR_SIGN (y) : MPFR_SIGN_POS); - } - /* Check Underflow */ - else if (MPFR_UNDERFLOW (flags)) + if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))) { + mpz_t z; + + /* Internal overflow or underflow. However the approximation error has + * not been taken into account. So, let's solve this problem by using + * mpfr_pow_z, which can handle it. This case could be improved in the + * future, without having to use mpfr_pow_z. + */ + MPFR_LOG_MSG (("Internal overflow or underflow," + " let's use mpfr_pow_z.\n", 0)); mpfr_clear (res); MPFR_SAVE_EXPO_FREE (expo); - return mpfr_underflow (x, rnd == GMP_RNDN ? GMP_RNDZ : rnd, - (n % 2) ? MPFR_SIGN(y) : MPFR_SIGN_POS); + mpz_init (z); + mpz_set_ui (z, n); + inexact = mpfr_pow_z (y, x, z, rnd); + mpz_clear (z); + return inexact; } - inexact = mpfr_set (x, res, rnd); + inexact = mpfr_set (y, res, rnd); mpfr_clear (res); MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (x, inexact, rnd); + return mpfr_check_range (y, inexact, rnd); } |