diff options
| -rw-r--r-- | src/compound.c | 101 | ||||
| -rw-r--r-- | tests/tcompound.c | 14 |
2 files changed, 110 insertions, 5 deletions
diff --git a/src/compound.c b/src/compound.c index b53ba2007..97d14f13a 100644 --- a/src/compound.c +++ b/src/compound.c @@ -20,18 +20,51 @@ along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ +#define MPFR_NEED_LONGLONG_H /* needed for MPFR_INT_CEIL_LOG2 */ #include "mpfr-impl.h" -/* put in z the correctly rounded value of (1+x)^n */ +/* assuming |(1+x)^n - 1| < 1/4*ulp(1), return correct rounding, + where s is the sign of n*log2(1+x) */ +static int +mpfr_compound_near_one (mpfr_ptr y, int s, mpfr_rnd_t rnd_mode) +{ + mpfr_set_ui (y, 1, rnd_mode); /* exact */ + if (rnd_mode == MPFR_RNDN || MPFR_RNDF + || (s > 0 && (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD)) + || (s < 0 && (rnd_mode == MPFR_RNDA || rnd_mode == MPFR_RNDU))) + { + /* round toward 1 */ + return -s; + } + else if (s > 0) /* necessarily RNDA or RNDU */ + { + /* round toward +Inf */ + mpfr_nextabove (y); + return +1; + } + else /* necessarily s < 0 and RNDZ or RNDD */ + { + /* round toward 0 */ + mpfr_nextbelow (y); + return -1; + } +} + +/* put in y the correctly rounded value of (1+x)^n */ int mpfr_compound (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode) { - int inex, compared; + int inexact, compared, k, nloop; + mpfr_t t; + mpfr_exp_t e; + mpfr_prec_t prec; + MPFR_ZIV_DECL (loop); + MPFR_SAVE_EXPO_DECL (expo); MPFR_LOG_FUNC (("x[%Pu]=%.*Rg n=%ld rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, n, rnd_mode), - ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex)); + ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inexact)); /* Special cases */ if (MPFR_IS_SINGULAR (x)) @@ -88,6 +121,64 @@ mpfr_compound (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode) } } - inex = 0; - return mpfr_set_si (y, -17, rnd_mode); + MPFR_SAVE_EXPO_MARK (expo); + + prec = MPFR_PREC(y); + prec += MPFR_INT_CEIL_LOG2 (prec) + 6; + + mpfr_init2 (t, prec); + + k = MPFR_INT_CEIL_LOG2((n > 0) ? n : -n); + + MPFR_ZIV_INIT (loop, prec); + for (nloop = 0; ; nloop++) + { + /* we compute (1+x)^n as 2^(n*log2p1(x)) */ + inexact = mpfr_log2p1 (t, x, MPFR_RNDN); + e = MPFR_GET_EXP(t); + /* |t - log2(1+x)| <= 1/2*ulp(t) = 2^(e-prec-1) */ + inexact |= mpfr_mul_ui (t, t, n, MPFR_RNDN); + /* |t - n*log2(1+x)| <= 2^(e2-prec-1) + n*2^(e-prec-1) + <= 2^(e2-prec-1) + 2^(e+k-prec-1) <= 2^(e+k-prec) + where n < 2^k, and e2 is the new exponent of t. */ + MPFR_ASSERTD(MPFR_GET_EXP(t) <= e + k); + e += k; + /* |t - n*log2(1+x)| <= 2^(e-prec) */ + /* detect overflow */ + if (nloop == 0 && mpfr_cmp_ui (t, __gmpfr_emax) >= 0) + { + mpfr_clear (t); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_overflow (y, rnd_mode, 1); + } + /* Detect cases where result is 1 or 1+ulp(1) or 1-1/2*ulp(1): + |2^t - 1| = |exp(t*log(2)) - 1| <= |t|*log(2) < |t| */ + if (nloop == 0 && MPFR_GET_EXP(t) < - (mpfr_exp_t) MPFR_PREC(y)) + { + int signt = MPFR_SIGN(t); + /* since ulp(1) = 2^(1-PREC(y)), we have |t| < 1/4*ulp(1) */ + mpfr_clear (t); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_compound_near_one (y, signt, rnd_mode); + } + inexact |= mpfr_exp2 (t, t, MPFR_RNDA); + /* |t - (1+x)^n| <= ulp(t) + |t|*log(2)*2^(e-prec) + < 2^(EXP(t)-prec) + 2^(EXP(t)+e-prec) */ + e = (e >= 0) ? e + 1 : 1; + /* now |t - (1+x)^n| < 2^(EXP(t)+e-prec) */ + + if (MPFR_LIKELY (inexact == 0 || + MPFR_CAN_ROUND (t, prec - e, MPFR_PREC(y), rnd_mode))) + break; + + MPFR_ZIV_NEXT (loop, prec); + mpfr_set_prec (t, prec); + } + inexact = mpfr_set (y, t, rnd_mode); + + MPFR_ZIV_FREE (loop); + mpfr_clear (t); + + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd_mode); } diff --git a/tests/tcompound.c b/tests/tcompound.c index e66661863..43f481d84 100644 --- a/tests/tcompound.c +++ b/tests/tcompound.c @@ -158,6 +158,20 @@ check_ieee754 (void) exit (1); } + /* hard-coded test: x is the 32-bit nearest approximation of 17/42 */ + mpfr_set_prec (x, 32); + mpfr_set_prec (y, 32); + mpfr_set_ui_2exp (x, 3476878287UL, -33, MPFR_RNDN); + mpfr_compound (y, x, 12, MPFR_RNDN); + mpfr_set_ui_2exp (x, 1981447393UL, -25, MPFR_RNDN); + if (!mpfr_equal_p (y, x)) + { + printf ("Error for compound(3476878287/2^33,12)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); } |