summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c36
1 files changed, 28 insertions, 8 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c b/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c
index 3e5bb571d88..a95d3caa548 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c
@@ -1,6 +1,6 @@
/* mpfr_rec_sqrt -- inverse square root
-Copyright 2008-2019 Free Software Foundation, Inc.
+Copyright 2008-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -165,7 +165,7 @@ mpfr_mpn_rec_sqrt (mpfr_limb_ptr x, mpfr_prec_t p,
MPFR_ASSERTD((a[an - 1] & MPFR_LIMB_HIGHBIT) != 0);
/* We should have enough bits in one limb and GMP_NUMB_BITS should be even.
Since that does not depend on MPFR, we always check this. */
- MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS >= 12 && (GMP_NUMB_BITS & 1) == 0);
+ MPFR_STAT_STATIC_ASSERT ((GMP_NUMB_BITS & 1) == 0);
/* {a, an} and {x, n} should not overlap */
MPFR_ASSERTD((a + an <= x) || (x + n <= a));
MPFR_ASSERTD(p >= 11);
@@ -179,15 +179,34 @@ mpfr_mpn_rec_sqrt (mpfr_limb_ptr x, mpfr_prec_t p,
if (p == 11) /* should happen only from recursive calls */
{
unsigned long i, ab, ac;
- mp_limb_t t;
+ unsigned int t;
/* take the 12+as most significant bits of A */
+#if GMP_NUMB_BITS >= 16
i = a[an - 1] >> (GMP_NUMB_BITS - (12 + as));
+#else
+ MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 8);
+ {
+ unsigned int a12 = a[an - 1] << 8;
+ if (an >= 2)
+ a12 |= a[an - 2];
+ MPFR_ASSERTN(GMP_NUMB_BITS >= 4 + as);
+ i = a12 >> (GMP_NUMB_BITS - (4 + as));
+ }
+#endif
/* if one wants faithful rounding for p=11, replace #if 0 by #if 1 */
ab = i >> 4;
ac = (ab & 0x3F0) | (i & 0x0F);
- t = (mp_limb_t) T1[ab - 0x80] + (mp_limb_t) T2[ac - 0x80];
- x[0] = t << (GMP_NUMB_BITS - p);
+ t = T1[ab - 0x80] + T2[ac - 0x80]; /* fits on 16 bits */
+#if GMP_NUMB_BITS >= 16
+ /* x has only one limb */
+ x[0] = (mp_limb_t) t << (GMP_NUMB_BITS - p);
+#else
+ MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 8);
+ MPFR_ASSERTD (1024 <= t && t <= 2047);
+ x[1] = t >> 3; /* 128 <= x[1] <= 255 */
+ x[0] = MPFR_LIMB_LSHIFT(t, 5);
+#endif
}
else /* p >= 12 */
{
@@ -225,7 +244,7 @@ mpfr_mpn_rec_sqrt (mpfr_limb_ptr x, mpfr_prec_t p,
/* we need h+1+as bits of a */
ahn = LIMB_SIZE(h + 1 + as); /* number of high limbs of A
- needed for the recursive call*/
+ needed for the recursive call */
if (MPFR_UNLIKELY(ahn > an))
ahn = an;
mpfr_mpn_rec_sqrt (x + ln, h, a + an - ahn, ahn * GMP_NUMB_BITS, as);
@@ -373,9 +392,10 @@ mpfr_mpn_rec_sqrt (mpfr_limb_ptr x, mpfr_prec_t p,
thus we round u to nearest at bit pl-1 of u[0] */
if (pl > 0)
{
- cu = mpn_add_1 (u, u, un, u[0] & (MPFR_LIMB_ONE << (pl - 1)));
+ cu = mpn_add_1 (u, u, un,
+ u[0] & MPFR_LIMB_LSHIFT(MPFR_LIMB_ONE, pl - 1));
/* mask bits 0..pl-1 of u[0] */
- u[0] &= ~MPFR_LIMB_MASK(pl);
+ u[0] &= MPFR_LIMB(~MPFR_LIMB_MASK(pl));
}
else /* round bit is in u[-1] */
cu = mpn_add_1 (u, u, un, u[-1] >> (GMP_NUMB_BITS - 1));