diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/rem1.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/rem1.c | 24 |
1 files changed, 18 insertions, 6 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/rem1.c b/Build/source/libs/mpfr/mpfr-src/src/rem1.c index 7c3c5591f18..9e75e428798 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/rem1.c +++ b/Build/source/libs/mpfr/mpfr-src/src/rem1.c @@ -2,7 +2,7 @@ mpfr_fmod -- compute the floating-point remainder of x/y mpfr_remquo and mpfr_remainder -- argument reduction functions -Copyright 2007-2019 Free Software Foundation, Inc. +Copyright 2007-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -22,7 +22,7 @@ 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. */ -# include "mpfr-impl.h" +#include "mpfr-impl.h" /* we return as many bits as we can, keeping just one bit for the sign */ # define WANTED_BITS (sizeof(long) * CHAR_BIT - 1) @@ -100,9 +100,11 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q, mpz_abs (my, my); q_is_odd = 0; - /* divide my by 2^k if possible to make operations mod my easier */ + /* Divide my by 2^k if possible to make operations mod my easier. + Since my comes from a regular MPFR number, due to the constraints on the + exponent and the precision, there can be no integer overflow below. */ { - unsigned long k = mpz_scan1 (my, 0); + mpfr_exp_t k = mpz_scan1 (my, 0); ey += k; mpz_fdiv_q_2exp (my, my, k); } @@ -157,8 +159,18 @@ mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q, which is obtained by dividing by 2Y. */ mpz_mul_2exp (my, my, 1); /* 2Y */ - mpz_set_ui (r, 2); - mpz_powm_ui (r, r, ex - ey, my); /* 2^(ex-ey) mod my */ + /* Warning: up to GMP 6.2.0, mpz_powm_ui is not optimized when BASE^EXP + has about the same size as MOD, in which case it should first compute + BASE^EXP exactly, then reduce it modulo MOD: + https://gmplib.org/list-archives/gmp-bugs/2020-February/004736.html + Thus when 2^(ex-ey) is less than my^3, we use this algorithm. */ + if (ex - ey > 3 * mpz_sizeinbase (my, 2)) + { + mpz_set_ui (r, 2); + mpz_powm_ui (r, r, ex - ey, my); /* 2^(ex-ey) mod my */ + } + else + mpz_ui_pow_ui (r, 2, ex - ey); mpz_mul (r, r, mx); mpz_mod (r, r, my); |