diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2010-08-18 11:55:19 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2010-08-18 11:55:19 +0000 |
commit | 8614f8fbe5b58ce91e8b6112736ac56a0e81c0c1 (patch) | |
tree | 78b463345494e39a37224eb2c4cd119ac8f9a120 /src/gmp_op.c | |
parent | 4d94693ff5087575daef1befa851bb1fc0c21199 (diff) | |
download | mpfr-8614f8fbe5b58ce91e8b6112736ac56a0e81c0c1.tar.gz |
[src/gmp_op.c] Added function mpfr_muldiv_z (currently static -- should
it be in the API?) that computes y = RND(x*n/d), where n and d are mpz
integers. Changed mpfr_mul_q and mpfr_div_q to use this function.
Note: the code of the general case is currently the same as the old
mpfr_mul_q/mpfr_div_q code, thus needs to be fixed.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@7095 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src/gmp_op.c')
-rw-r--r-- | src/gmp_op.c | 89 |
1 files changed, 47 insertions, 42 deletions
diff --git a/src/gmp_op.c b/src/gmp_op.c index a738c2427..65f49158b 100644 --- a/src/gmp_op.c +++ b/src/gmp_op.c @@ -98,35 +98,56 @@ mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z) /* FIXME [VL] (for mpfr_mul_q and mpfr_div_q): an intermediate overflow doesn't necessarily imply an overflow on the final result. Moreover - the exponent range should be extended in the usual way. - To fix this, I think that these functions should call a common - function mpfr_muldiv_z: - res = mpfr_muldiv_z (y, x, mpq_numref(z), mpq_denref(z), rnd_mode); - res = mpfr_muldiv_z (y, x, mpq_denref(z), mpq_numref(z), rnd_mode); - respectively, so that all the work isn't done twice. */ - -int -mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) + the exponent range should be extended in the usual way (unless the + computations are done in mpn, mpz or mpq). */ + +/* Compute y = RND(x*n/d), where n and d are mpz integers. + An integer 0 is assumed to have a positive sign. + This function is used by mpfr_mul_q and mpfr_div_q. + Note: the status of the rational 0/(-1) is not clear (if there is + a signed infinity, there should be a signed zero). But infinities + are not currently supported/documented in GMP, and if the rational + is canonicalized as it should be, the case 0/(-1) cannot occur. */ +static int +mpfr_muldiv_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr n, mpz_srcptr d, + mpfr_rnd_t rnd_mode) { - mpfr_t tmp; - int res; - mpfr_prec_t p; - - if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) - return mpfr_mul_ui (y, x, 0, rnd_mode); + if (MPFR_UNLIKELY (mpz_sgn (n) == 0)) + { + if (MPFR_UNLIKELY (mpz_sgn (d) == 0)) + MPFR_SET_NAN (y); + else + { + mpfr_mul_ui (y, x, 0, MPFR_RNDN); /* exact: +0, -0 or NaN */ + if (MPFR_UNLIKELY (mpz_sgn (d) < 0)) + MPFR_CHANGE_SIGN (y); + } + return 0; + } + else if (MPFR_UNLIKELY (mpz_sgn (d) == 0)) + { + mpfr_div_ui (y, x, 0, MPFR_RNDN); /* exact: +Inf, -Inf or NaN */ + if (MPFR_UNLIKELY (mpz_sgn (n) < 0)) + MPFR_CHANGE_SIGN (y); + return 0; + } else { - MPFR_MPZ_SIZEINBASE2 (p, mpq_numref (z)); + mpfr_prec_t p; + mpfr_t tmp; + int res; + + MPFR_MPZ_SIZEINBASE2 (p, n); mpfr_init2 (tmp, MPFR_PREC (x) + p); - res = mpfr_mul_z (tmp, x, mpq_numref(z), MPFR_RNDN ); + res = mpfr_mul_z (tmp, x, n, MPFR_RNDN); if (MPFR_UNLIKELY (res != 0)) { - /* overflow case */ + /* overflow case - FIXME */ MPFR_ASSERTD (mpfr_inf_p (tmp)); mpfr_set (y, tmp, MPFR_RNDN); /* exact */ } else - res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode); + res = mpfr_div_z (y, tmp, d, rnd_mode); mpfr_clear (tmp); return res; @@ -134,31 +155,15 @@ mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) } int -mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) +mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) { - mpfr_t tmp; - int res; - mpfr_prec_t p; - - if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) - return mpfr_div_ui (y, x, 0, rnd_mode); - else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0)) - p = 0; - else - MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z)); - mpfr_init2 (tmp, MPFR_PREC(x) + p); - res = mpfr_mul_z (tmp, x, mpq_denref(z), MPFR_RNDN ); - if (MPFR_UNLIKELY (res != 0)) - { - /* overflow case */ - MPFR_ASSERTD (mpfr_inf_p (tmp)); - mpfr_set (y, tmp, MPFR_RNDN); /* exact */ - } - else - res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode); + return mpfr_muldiv_z (y, x, mpq_numref (z), mpq_denref (z), rnd_mode); +} - mpfr_clear (tmp); - return res; +int +mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) +{ + return mpfr_muldiv_z (y, x, mpq_denref (z), mpq_numref (z), rnd_mode); } int |