diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2016-03-07 13:18:57 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2016-03-07 13:18:57 +0000 |
commit | f7fd80054d9fb1620b6c65db276325f35d86f7ee (patch) | |
tree | cc7052c8a052eb2e0981b9c47f9b3902afd51561 /Build/source/libs/mpfr/mpfr-src/src/rem1.c | |
parent | ab4c67522df470475142ea5802c2951d0319ccac (diff) |
mpfr 3.1.4
git-svn-id: svn://tug.org/texlive/trunk@39960 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/rem1.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/rem1.c | 46 |
1 files changed, 37 insertions, 9 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/rem1.c b/Build/source/libs/mpfr/mpfr-src/src/rem1.c index a977a0924c3..1299396e32b 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/rem1.c +++ b/Build/source/libs/mpfr/mpfr-src/src/rem1.c @@ -2,8 +2,8 @@ mpfr_fmod -- compute the floating-point remainder of x/y mpfr_remquo and mpfr_remainder -- argument reduction functions -Copyright 2007-2015 Free Software Foundation, Inc. -Contributed by the AriC and Caramel projects, INRIA. +Copyright 2007-2016 Free Software Foundation, Inc. +Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -59,6 +59,7 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q, mpfr_exp_t ex, ey; int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x); mpz_t mx, my, r; + int tiny = 0; MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ); @@ -109,13 +110,27 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q, if (ex <= ey) { /* q = x/y = mx/(my*2^(ey-ex)) */ - mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */ - if (rnd_q == MPFR_RNDZ) - /* 0 <= |r| <= |my|, r has the same sign as mx */ - mpz_tdiv_qr (mx, r, mx, my); + + /* First detect cases where q=0, to avoid creating a huge number + my*2^(ey-ex): if sx = mpz_sizeinbase (mx, 2) and sy = + mpz_sizeinbase (my, 2), we have x < 2^(ex + sx) and + y >= 2^(ey + sy - 1), thus if ex + sx <= ey + sy - 1 + the quotient is 0 */ + if (ex + (mpfr_exp_t) mpz_sizeinbase (mx, 2) < + ey + (mpfr_exp_t) mpz_sizeinbase (my, 2)) + { + tiny = 1; + mpz_set (r, mx); + mpz_set_ui (mx, 0); + } else - /* 0 <= |r| <= |my|, r has the same sign as my */ - mpz_fdiv_qr (mx, r, mx, my); + { + mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */ + + /* since mx > 0 and my > 0, we can use mpz_tdiv_qr in all cases */ + mpz_tdiv_qr (mx, r, mx, my); + /* 0 <= |r| <= |my|, r has the same sign as mx */ + } if (rnd_q == MPFR_RNDN) q_is_odd = mpz_tstbit (mx, 0); @@ -181,7 +196,20 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q, /* FIXME: the comparison 2*r < my could be done more efficiently at the mpn level */ mpz_mul_2exp (r, r, 1); - compare = mpz_cmpabs (r, my); + /* if tiny=1, we should compare r with my*2^(ey-ex) */ + if (tiny) + { + if (ex + (mpfr_exp_t) mpz_sizeinbase (r, 2) < + ey + (mpfr_exp_t) mpz_sizeinbase (my, 2)) + compare = 0; /* r*2^ex < my*2^ey */ + else + { + mpz_mul_2exp (my, my, ey - ex); + compare = mpz_cmpabs (r, my); + } + } + else + compare = mpz_cmpabs (r, my); mpz_fdiv_q_2exp (r, r, 1); compare = ((compare > 0) || ((rnd_q == MPFR_RNDN) && (compare == 0) && q_is_odd)); |