diff options
| author | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2004-02-16 18:02:42 +0000 |
|---|---|---|
| committer | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2004-02-16 18:02:42 +0000 |
| commit | b6d4dd5a15b2b0a1166c90906259cd82f41cf437 (patch) | |
| tree | 047d78025ea253f20e16e3b21bcb57c199314757 /gmp_op.c | |
| parent | 78ddd7e74f1654c2b6fb7f6185ab9d469671084e (diff) | |
| download | mpfr-b6d4dd5a15b2b0a1166c90906259cd82f41cf437.tar.gz | |
Fix the computing of the error for mpq_add and mpq_sub.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2739 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'gmp_op.c')
| -rw-r--r-- | gmp_op.c | 21 |
1 files changed, 16 insertions, 5 deletions
@@ -105,8 +105,9 @@ mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) int mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) { - mpfr_t t,q; + mpfr_t t,q; mp_prec_t p = MPFR_PREC(y)+10; + mp_exp_t err; int res; if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) @@ -133,9 +134,13 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode) do { mpfr_set_q(q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ mpfr_add(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ - /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) */ - res = mpfr_can_round(t, p-1, GMP_RNDN, GMP_RNDZ, - MPFR_PREC(y) + (rnd_mode == GMP_RNDN) ); + /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) + If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) + <= 2^(EXP(q)-EXP(t)) + If EXP(q)-EXP(t)<0, <= 2^0 */ + err = (mp_exp_t) p - 1 - MAX(MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); + res = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + MPFR_PREC(y) + (rnd_mode == GMP_RNDN) ); if (!res) { p += BITS_PER_MP_LIMB; @@ -154,6 +159,7 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) mpfr_t t,q; mp_prec_t p = MPFR_PREC(y)+10; int res; + mp_exp_t err; if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { @@ -185,7 +191,12 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode) do { mpfr_set_q(q, z, GMP_RNDN); /* Error <= 1/2 ulp(q) */ mpfr_sub(t, x, q, GMP_RNDN); /* Error <= 1/2 ulp(t) */ - res = mpfr_can_round(t, p-1, GMP_RNDN, GMP_RNDZ, + /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) + If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) + <= 2^(EXP(q)-EXP(t)) + If EXP(q)-EXP(t)<0, <= 2^0 */ + err = (mp_exp_t) p - 1 - MAX(MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); + res = mpfr_can_round(t, err, GMP_RNDN, GMP_RNDZ, MPFR_PREC(y) + (rnd_mode == GMP_RNDN) ); if (!res) { |