summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sqr.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sqr.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sqr.c91
1 files changed, 49 insertions, 42 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sqr.c b/Build/source/libs/mpfr/mpfr-src/src/sqr.c
index 56c21453021..0ab414ba113 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sqr.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sqr.c
@@ -1,6 +1,6 @@
/* mpfr_sqr -- Floating-point square
-Copyright 2004-2019 Free Software Foundation, Inc.
+Copyright 2004-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -26,8 +26,10 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#if !defined(MPFR_GENERIC_ABI) && (GMP_NUMB_BITS == 32 || GMP_NUMB_BITS == 64)
/* Special code for prec(a) < GMP_NUMB_BITS and prec(b) <= GMP_NUMB_BITS.
- Note: this function was copied from mpfr_mul_1 in file mul.c, thus any change
- here should be done also in mpfr_mul_1. */
+ Note: this function was copied from mpfr_mul_1 in file mul.c, thus any
+ change here should be done also in mpfr_mul_1.
+ Although this function works as soon as prec(a) < GMP_NUMB_BITS and
+ prec(b) <= GMP_NUMB_BITS, we use it for prec(a)=prec(b) < GMP_NUMB_BITS. */
static int
mpfr_sqr_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
{
@@ -75,16 +77,17 @@ mpfr_sqr_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
possible for p=53, because the largest significand is 6369051672525772
but its square has only 52 leading ones. For p=24 it is possible,
with b = 11863283, whose square has 24 leading ones. */
- if ((ax == __gmpfr_emin - 1) && (ap[0] == ~mask) &&
+ if (ax == __gmpfr_emin - 1 && ap[0] == ~mask &&
((rnd_mode == MPFR_RNDN && rb) ||
- (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb))))
+ (MPFR_IS_LIKE_RNDA (rnd_mode, 0) && (rb | sb))))
goto rounding; /* no underflow */
/* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2)
we have to change to RNDZ. This corresponds to:
(a) either ax < emin - 1
(b) or ax = emin - 1 and ap[0] = 1000....000 and rb = sb = 0 */
if (rnd_mode == MPFR_RNDN &&
- (ax < __gmpfr_emin - 1 || (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0)))
+ (ax < __gmpfr_emin - 1 ||
+ (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0)))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS);
}
@@ -104,7 +107,7 @@ mpfr_sqr_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
else
goto add_one_ulp;
}
- else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ else if (MPFR_IS_LIKE_RNDZ (rnd_mode, 0))
{
truncate:
MPFR_ASSERTD(ax >= __gmpfr_emin);
@@ -179,7 +182,8 @@ mpfr_sqr_1n (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
(a) either ax < emin - 1
(b) or ax = emin - 1 and ap[0] = 1000....000 and rb = sb = 0 */
if (rnd_mode == MPFR_RNDN &&
- (ax < __gmpfr_emin - 1 || (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0)))
+ (ax < __gmpfr_emin - 1 ||
+ (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0)))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS);
}
@@ -199,7 +203,7 @@ mpfr_sqr_1n (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
else
goto add_one_ulp;
}
- else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ else if (MPFR_IS_LIKE_RNDZ (rnd_mode, 0))
{
truncate:
MPFR_ASSERTD(ax >= __gmpfr_emin);
@@ -292,24 +296,24 @@ mpfr_sqr_2 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
a >= 0.111...111[1]*2^(emin-1), there is no underflow. */
if (MPFR_UNLIKELY(ax < __gmpfr_emin))
{
- /* Note: like for mpfr_sqr_2, the case
+ /* Note: like for mpfr_sqr_1, the case
0.111...111*2^(emin-1) < a < 2^(emin-1) is not possible when emin is
odd, since (modulo a shift) this would imply 1-2^(-p) < a = b^2 < 1,
and this is not possible with 1-2^(-p) <= b < 1.
For emin even, it is possible for some values of p, for example for
p=69 with b=417402170410649030795*2^k. */
- if ((ax == __gmpfr_emin - 1) &&
- (ap[1] == MPFR_LIMB_MAX) &&
- (ap[0] == ~mask) &&
+ if (ax == __gmpfr_emin - 1 &&
+ ap[1] == MPFR_LIMB_MAX &&
+ ap[0] == ~mask &&
((rnd_mode == MPFR_RNDN && rb) ||
- (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb))))
+ (MPFR_IS_LIKE_RNDA (rnd_mode, 0) && (rb | sb))))
goto rounding; /* no underflow */
- /* for RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2)
- we have to change to RNDZ */
+ /* for RNDN, mpfr_underflow always rounds away, thus for
+ |a| <= 2^(emin-2) we have to change to RNDZ */
if (rnd_mode == MPFR_RNDN &&
(ax < __gmpfr_emin - 1 ||
(ap[1] == MPFR_LIMB_HIGHBIT && ap[0] == 0 && (rb | sb) == 0)))
- rnd_mode = MPFR_RNDZ;
+ rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN_POS);
}
@@ -328,7 +332,7 @@ mpfr_sqr_2 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
else
goto add_one_ulp;
}
- else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ else if (MPFR_IS_LIKE_RNDZ (rnd_mode, 0))
{
truncate:
MPFR_ASSERTD(ax >= __gmpfr_emin);
@@ -401,13 +405,13 @@ mpfr_sqr_3 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
sb = sb2 = 1; /* result cannot be exact in that case */
else
{
- mp_limb_t p[6];
- mpn_sqr (p, bp, 3);
- a2 = p[5];
- a1 = p[4];
- a0 = p[3];
- sb = p[2];
- sb2 = p[1] | p[0];
+ mp_limb_t t[6];
+ mpn_sqr (t, bp, 3);
+ a2 = t[5];
+ a1 = t[4];
+ a0 = t[3];
+ sb = t[2];
+ sb2 = t[1] | t[0];
}
if (a2 < MPFR_LIMB_HIGHBIT)
{
@@ -435,12 +439,12 @@ mpfr_sqr_3 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
a >= 0.111...111[1]*2^(emin-1), there is no underflow. */
if (MPFR_UNLIKELY(ax < __gmpfr_emin))
{
- if ((ax == __gmpfr_emin - 1) &&
- (ap[2] == MPFR_LIMB_MAX) &&
- (ap[1] == MPFR_LIMB_MAX) &&
- (ap[0] == ~mask) &&
+ if (ax == __gmpfr_emin - 1 &&
+ ap[2] == MPFR_LIMB_MAX &&
+ ap[1] == MPFR_LIMB_MAX &&
+ ap[0] == ~mask &&
((rnd_mode == MPFR_RNDN && rb) ||
- (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb))))
+ (MPFR_IS_LIKE_RNDA (rnd_mode, 0) && (rb | sb))))
goto rounding; /* no underflow */
/* for RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2)
we have to change to RNDZ */
@@ -467,7 +471,7 @@ mpfr_sqr_3 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, mpfr_prec_t p)
else
goto add_one_ulp;
}
- else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ else if (MPFR_IS_LIKE_RNDZ (rnd_mode, 0))
{
truncate:
MPFR_ASSERTD(ax >= __gmpfr_emin);
@@ -535,19 +539,20 @@ mpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
bq = MPFR_GET_PREC(b);
#if !defined(MPFR_GENERIC_ABI) && (GMP_NUMB_BITS == 32 || GMP_NUMB_BITS == 64)
- if (aq < GMP_NUMB_BITS && bq <= GMP_NUMB_BITS)
- return mpfr_sqr_1 (a, b, rnd_mode, aq);
+ if (aq == bq)
+ {
+ if (aq < GMP_NUMB_BITS)
+ return mpfr_sqr_1 (a, b, rnd_mode, aq);
- if (GMP_NUMB_BITS < aq && aq < 2 * GMP_NUMB_BITS
- && GMP_NUMB_BITS < bq && bq <= 2 * GMP_NUMB_BITS)
- return mpfr_sqr_2 (a, b, rnd_mode, aq);
+ if (GMP_NUMB_BITS < aq && aq < 2 * GMP_NUMB_BITS)
+ return mpfr_sqr_2 (a, b, rnd_mode, aq);
- if (aq == GMP_NUMB_BITS && bq <= GMP_NUMB_BITS)
- return mpfr_sqr_1n (a, b, rnd_mode);
+ if (aq == GMP_NUMB_BITS)
+ return mpfr_sqr_1n (a, b, rnd_mode);
- if (2 * GMP_NUMB_BITS < aq && aq < 3 * GMP_NUMB_BITS
- && 2 * GMP_NUMB_BITS < bq && bq <= 3 * GMP_NUMB_BITS)
- return mpfr_sqr_3 (a, b, rnd_mode, aq);
+ if (2 * GMP_NUMB_BITS < aq && aq < 3 * GMP_NUMB_BITS)
+ return mpfr_sqr_3 (a, b, rnd_mode, aq);
+ }
#endif
ax = 2 * MPFR_GET_EXP (b);
@@ -558,6 +563,8 @@ mpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
2*bn or 2*bn-1 */
if (MPFR_UNLIKELY(bn > MPFR_SQR_THRESHOLD))
+ /* the following line should not be replaced by mpfr_sqr,
+ otherwise we'll get an infinite loop! */
return mpfr_mul (a, b, b, rnd_mode);
MPFR_TMP_MARK(marker);
@@ -585,7 +592,7 @@ mpfr_sqr (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
MPFR_TMP_FREE(marker);
{
- mpfr_exp_t ax2 = ax + (mpfr_exp_t) (b1 - 1 + cc);
+ mpfr_exp_t ax2 = ax + ((int) b1 - 1 + cc);
if (MPFR_UNLIKELY( ax2 > __gmpfr_emax))
return mpfr_overflow (a, rnd_mode, MPFR_SIGN_POS);
if (MPFR_UNLIKELY( ax2 < __gmpfr_emin))