Revert "Remove gmp from AOS"
This reverts commit f37c97684f0910a3f241394549392f00145ab0f7.
We need gmp for SymEngine for symbolicmanipultion in C++
Change-Id: Ia13216d1715cf96944f7b4f422b7a799f921d4a4
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/third_party/gmp/mpn/generic/rootrem.c b/third_party/gmp/mpn/generic/rootrem.c
new file mode 100644
index 0000000..a79099e
--- /dev/null
+++ b/third_party/gmp/mpn/generic/rootrem.c
@@ -0,0 +1,515 @@
+/* mpn_rootrem(rootp,remp,ap,an,nth) -- Compute the nth root of {ap,an}, and
+ store the truncated integer part at rootp and the remainder at remp.
+
+ Contributed by Paul Zimmermann (algorithm) and
+ Paul Zimmermann and Torbjorn Granlund (implementation).
+ Marco Bodrato wrote logbased_root to seed the loop.
+
+ THE FUNCTIONS IN THIS FILE ARE INTERNAL, AND HAVE MUTABLE INTERFACES. IT'S
+ ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT'S ALMOST
+ GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 2002, 2005, 2009-2012, 2015 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of either:
+
+ * the GNU Lesser General Public License as published by the Free
+ Software Foundation; either version 3 of the License, or (at your
+ option) any later version.
+
+or
+
+ * the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any
+ later version.
+
+or both in parallel, as here.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received copies of the GNU General Public License and the
+GNU Lesser General Public License along with the GNU MP Library. If not,
+see https://www.gnu.org/licenses/. */
+
+/* FIXME:
+ This implementation is not optimal when remp == NULL, since the complexity
+ is M(n), whereas it should be M(n/k) on average.
+*/
+
+#include <stdio.h> /* for NULL */
+
+#include "gmp-impl.h"
+#include "longlong.h"
+
+static mp_size_t mpn_rootrem_internal (mp_ptr, mp_ptr, mp_srcptr, mp_size_t,
+ mp_limb_t, int);
+
+#define MPN_RSHIFT(rp,up,un,cnt) \
+ do { \
+ if ((cnt) != 0) \
+ mpn_rshift (rp, up, un, cnt); \
+ else \
+ { \
+ MPN_COPY_INCR (rp, up, un); \
+ } \
+ } while (0)
+
+#define MPN_LSHIFT(cy,rp,up,un,cnt) \
+ do { \
+ if ((cnt) != 0) \
+ cy = mpn_lshift (rp, up, un, cnt); \
+ else \
+ { \
+ MPN_COPY_DECR (rp, up, un); \
+ cy = 0; \
+ } \
+ } while (0)
+
+
+/* Put in {rootp, ceil(un/k)} the kth root of {up, un}, rounded toward zero.
+ If remp <> NULL, put in {remp, un} the remainder.
+ Return the size (in limbs) of the remainder if remp <> NULL,
+ or a non-zero value iff the remainder is non-zero when remp = NULL.
+ Assumes:
+ (a) up[un-1] is not zero
+ (b) rootp has at least space for ceil(un/k) limbs
+ (c) remp has at least space for un limbs (in case remp <> NULL)
+ (d) the operands do not overlap.
+
+ The auxiliary memory usage is 3*un+2 if remp = NULL,
+ and 2*un+2 if remp <> NULL. FIXME: This is an incorrect comment.
+*/
+mp_size_t
+mpn_rootrem (mp_ptr rootp, mp_ptr remp,
+ mp_srcptr up, mp_size_t un, mp_limb_t k)
+{
+ ASSERT (un > 0);
+ ASSERT (up[un - 1] != 0);
+ ASSERT (k > 1);
+
+ if (UNLIKELY (k == 2))
+ return mpn_sqrtrem (rootp, remp, up, un);
+ /* (un-1)/k > 2 <=> un > 3k <=> (un + 2)/3 > k */
+ if (remp == NULL && (un + 2) / 3 > k)
+ /* Pad {up,un} with k zero limbs. This will produce an approximate root
+ with one more limb, allowing us to compute the exact integral result. */
+ {
+ mp_ptr sp, wp;
+ mp_size_t rn, sn, wn;
+ TMP_DECL;
+ TMP_MARK;
+ wn = un + k;
+ sn = (un - 1) / k + 2; /* ceil(un/k) + 1 */
+ TMP_ALLOC_LIMBS_2 (wp, wn, /* will contain the padded input */
+ sp, sn); /* approximate root of padded input */
+ MPN_COPY (wp + k, up, un);
+ MPN_FILL (wp, k, 0);
+ rn = mpn_rootrem_internal (sp, NULL, wp, wn, k, 1);
+ /* The approximate root S = {sp,sn} is either the correct root of
+ {sp,sn}, or 1 too large. Thus unless the least significant limb of
+ S is 0 or 1, we can deduce the root of {up,un} is S truncated by one
+ limb. (In case sp[0]=1, we can deduce the root, but not decide
+ whether it is exact or not.) */
+ MPN_COPY (rootp, sp + 1, sn - 1);
+ TMP_FREE;
+ return rn;
+ }
+ else
+ {
+ return mpn_rootrem_internal (rootp, remp, up, un, k, 0);
+ }
+}
+
+#define LOGROOT_USED_BITS 8
+#define LOGROOT_NEEDS_TWO_CORRECTIONS 1
+#define LOGROOT_RETURNED_BITS (LOGROOT_USED_BITS + LOGROOT_NEEDS_TWO_CORRECTIONS)
+/* Puts in *rootp some bits of the k^nt root of the number
+ 2^bitn * 1.op ; where op represents the "fractional" bits.
+
+ The returned value is the number of bits of the root minus one;
+ i.e. an approximation of the root will be
+ (*rootp) * 2^(retval-LOGROOT_RETURNED_BITS+1).
+
+ Currently, only LOGROOT_USED_BITS bits of op are used (the implicit
+ one is not counted).
+ */
+static unsigned
+logbased_root (mp_ptr rootp, mp_limb_t op, mp_bitcnt_t bitn, mp_limb_t k)
+{
+ /* vlog=vector(256,i,floor((log(256+i)/log(2)-8)*256)-(i>255)) */
+ static const
+ unsigned char vlog[] = {1, 2, 4, 5, 7, 8, 9, 11, 12, 14, 15, 16, 18, 19, 21, 22,
+ 23, 25, 26, 27, 29, 30, 31, 33, 34, 35, 37, 38, 39, 40, 42, 43,
+ 44, 46, 47, 48, 49, 51, 52, 53, 54, 56, 57, 58, 59, 61, 62, 63,
+ 64, 65, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 80, 81, 82,
+ 83, 84, 85, 87, 88, 89, 90, 91, 92, 93, 94, 96, 97, 98, 99, 100,
+ 101, 102, 103, 104, 105, 106, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
+ 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 131, 132, 133, 134,
+ 135, 136, 137, 138, 139, 140, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
+ 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 162, 163, 164,
+ 165, 166, 167, 168, 169, 170, 171, 172, 173, 173, 174, 175, 176, 177, 178, 179,
+ 180, 181, 181, 182, 183, 184, 185, 186, 187, 188, 188, 189, 190, 191, 192, 193,
+ 194, 194, 195, 196, 197, 198, 199, 200, 200, 201, 202, 203, 204, 205, 205, 206,
+ 207, 208, 209, 209, 210, 211, 212, 213, 214, 214, 215, 216, 217, 218, 218, 219,
+ 220, 221, 222, 222, 223, 224, 225, 225, 226, 227, 228, 229, 229, 230, 231, 232,
+ 232, 233, 234, 235, 235, 236, 237, 238, 239, 239, 240, 241, 242, 242, 243, 244,
+ 245, 245, 246, 247, 247, 248, 249, 250, 250, 251, 252, 253, 253, 254, 255, 255};
+
+ /* vexp=vector(256,i,floor(2^(8+i/256)-256)-(i>255)) */
+ static const
+ unsigned char vexp[] = {0, 1, 2, 2, 3, 4, 4, 5, 6, 7, 7, 8, 9, 9, 10, 11,
+ 12, 12, 13, 14, 14, 15, 16, 17, 17, 18, 19, 20, 20, 21, 22, 23,
+ 23, 24, 25, 26, 26, 27, 28, 29, 30, 30, 31, 32, 33, 33, 34, 35,
+ 36, 37, 37, 38, 39, 40, 41, 41, 42, 43, 44, 45, 45, 46, 47, 48,
+ 49, 50, 50, 51, 52, 53, 54, 55, 55, 56, 57, 58, 59, 60, 61, 61,
+ 62, 63, 64, 65, 66, 67, 67, 68, 69, 70, 71, 72, 73, 74, 75, 75,
+ 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 86, 87, 88, 89, 90,
+ 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106,
+ 107, 108, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 119, 120, 121, 122,
+ 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138,
+ 139, 140, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 154, 155, 156,
+ 157, 158, 159, 160, 161, 163, 164, 165, 166, 167, 168, 169, 171, 172, 173, 174,
+ 175, 176, 178, 179, 180, 181, 182, 183, 185, 186, 187, 188, 189, 191, 192, 193,
+ 194, 196, 197, 198, 199, 200, 202, 203, 204, 205, 207, 208, 209, 210, 212, 213,
+ 214, 216, 217, 218, 219, 221, 222, 223, 225, 226, 227, 229, 230, 231, 232, 234,
+ 235, 236, 238, 239, 240, 242, 243, 245, 246, 247, 249, 250, 251, 253, 254, 255};
+ mp_bitcnt_t retval;
+
+ if (UNLIKELY (bitn > (~ (mp_bitcnt_t) 0) >> LOGROOT_USED_BITS))
+ {
+ /* In the unlikely case, we use two divisions and a modulo. */
+ retval = bitn / k;
+ bitn %= k;
+ bitn = (bitn << LOGROOT_USED_BITS |
+ vlog[op >> (GMP_NUMB_BITS - LOGROOT_USED_BITS)]) / k;
+ }
+ else
+ {
+ bitn = (bitn << LOGROOT_USED_BITS |
+ vlog[op >> (GMP_NUMB_BITS - LOGROOT_USED_BITS)]) / k;
+ retval = bitn >> LOGROOT_USED_BITS;
+ bitn &= (CNST_LIMB (1) << LOGROOT_USED_BITS) - 1;
+ }
+ ASSERT(bitn < CNST_LIMB (1) << LOGROOT_USED_BITS);
+ *rootp = CNST_LIMB(1) << (LOGROOT_USED_BITS - ! LOGROOT_NEEDS_TWO_CORRECTIONS)
+ | vexp[bitn] >> ! LOGROOT_NEEDS_TWO_CORRECTIONS;
+ return retval;
+}
+
+/* if approx is non-zero, does not compute the final remainder */
+static mp_size_t
+mpn_rootrem_internal (mp_ptr rootp, mp_ptr remp, mp_srcptr up, mp_size_t un,
+ mp_limb_t k, int approx)
+{
+ mp_ptr qp, rp, sp, wp, scratch;
+ mp_size_t qn, rn, sn, wn, nl, bn;
+ mp_limb_t save, save2, cy, uh;
+ mp_bitcnt_t unb; /* number of significant bits of {up,un} */
+ mp_bitcnt_t xnb; /* number of significant bits of the result */
+ mp_bitcnt_t b, kk;
+ mp_bitcnt_t sizes[GMP_NUMB_BITS + 1];
+ int ni;
+ int perf_pow;
+ unsigned ulz, snb, c, logk;
+ TMP_DECL;
+
+ /* MPN_SIZEINBASE_2EXP(unb, up, un, 1); --unb; */
+ uh = up[un - 1];
+ count_leading_zeros (ulz, uh);
+ ulz = ulz - GMP_NAIL_BITS + 1; /* Ignore the first 1. */
+ unb = (mp_bitcnt_t) un * GMP_NUMB_BITS - ulz;
+ /* unb is the (truncated) logarithm of the input U in base 2*/
+
+ if (unb < k) /* root is 1 */
+ {
+ rootp[0] = 1;
+ if (remp == NULL)
+ un -= (*up == CNST_LIMB (1)); /* Non-zero iif {up,un} > 1 */
+ else
+ {
+ mpn_sub_1 (remp, up, un, CNST_LIMB (1));
+ un -= (remp [un - 1] == 0); /* There should be at most one zero limb,
+ if we demand u to be normalized */
+ }
+ return un;
+ }
+ /* if (unb - k < k/2 + k/16) // root is 2 */
+
+ if (ulz == GMP_NUMB_BITS)
+ uh = up[un - 2];
+ else
+ uh = (uh << ulz & GMP_NUMB_MASK) | up[un - 1 - (un != 1)] >> (GMP_NUMB_BITS - ulz);
+ ASSERT (un != 1 || up[un - 1 - (un != 1)] >> (GMP_NUMB_BITS - ulz) == 1);
+
+ xnb = logbased_root (rootp, uh, unb, k);
+ snb = LOGROOT_RETURNED_BITS - 1;
+ /* xnb+1 is the number of bits of the root R */
+ /* snb+1 is the number of bits of the current approximation S */
+
+ kk = k * xnb; /* number of truncated bits in the input */
+
+ /* FIXME: Should we skip the next two loops when xnb <= snb ? */
+ for (uh = (k - 1) / 2, logk = 3; (uh >>= 1) != 0; ++logk )
+ ;
+ /* logk = ceil(log(k)/log(2)) + 1 */
+
+ /* xnb is the number of remaining bits to determine in the kth root */
+ for (ni = 0; (sizes[ni] = xnb) > snb; ++ni)
+ {
+ /* invariant: here we want xnb+1 total bits for the kth root */
+
+ /* if c is the new value of xnb, this means that we'll go from a
+ root of c+1 bits (say s') to a root of xnb+1 bits.
+ It is proved in the book "Modern Computer Arithmetic" by Brent
+ and Zimmermann, Chapter 1, that
+ if s' >= k*beta, then at most one correction is necessary.
+ Here beta = 2^(xnb-c), and s' >= 2^c, thus it suffices that
+ c >= ceil((xnb + log2(k))/2). */
+ if (xnb > logk)
+ xnb = (xnb + logk) / 2;
+ else
+ --xnb; /* add just one bit at a time */
+ }
+
+ *rootp >>= snb - xnb;
+ kk -= xnb;
+
+ ASSERT_ALWAYS (ni < GMP_NUMB_BITS + 1);
+ /* We have sizes[0] = b > sizes[1] > ... > sizes[ni] = 0 with
+ sizes[i] <= 2 * sizes[i+1].
+ Newton iteration will first compute sizes[ni-1] extra bits,
+ then sizes[ni-2], ..., then sizes[0] = b. */
+
+ TMP_MARK;
+ /* qp and wp need enough space to store S'^k where S' is an approximate
+ root. Since S' can be as large as S+2, the worst case is when S=2 and
+ S'=4. But then since we know the number of bits of S in advance, S'
+ can only be 3 at most. Similarly for S=4, then S' can be 6 at most.
+ So the worst case is S'/S=3/2, thus S'^k <= (3/2)^k * S^k. Since S^k
+ fits in un limbs, the number of extra limbs needed is bounded by
+ ceil(k*log2(3/2)/GMP_NUMB_BITS). */
+ /* THINK: with the use of logbased_root, maybe the constant is
+ 258/256 instead of 3/2 ? log2(258/256) < 1/89 < 1/64 */
+#define EXTRA 2 + (mp_size_t) (0.585 * (double) k / (double) GMP_NUMB_BITS)
+ TMP_ALLOC_LIMBS_3 (scratch, un + 1, /* used by mpn_div_q */
+ qp, un + EXTRA, /* will contain quotient and remainder
+ of R/(k*S^(k-1)), and S^k */
+ wp, un + EXTRA); /* will contain S^(k-1), k*S^(k-1),
+ and temporary for mpn_pow_1 */
+
+ if (remp == NULL)
+ rp = scratch; /* will contain the remainder */
+ else
+ rp = remp;
+ sp = rootp;
+
+ sn = 1; /* Initial approximation has one limb */
+
+ for (b = xnb; ni != 0; --ni)
+ {
+ /* 1: loop invariant:
+ {sp, sn} is the current approximation of the root, which has
+ exactly 1 + sizes[ni] bits.
+ {rp, rn} is the current remainder
+ {wp, wn} = {sp, sn}^(k-1)
+ kk = number of truncated bits of the input
+ */
+
+ /* Since each iteration treats b bits from the root and thus k*b bits
+ from the input, and we already considered b bits from the input,
+ we now have to take another (k-1)*b bits from the input. */
+ kk -= (k - 1) * b; /* remaining input bits */
+ /* {rp, rn} = floor({up, un} / 2^kk) */
+ rn = un - kk / GMP_NUMB_BITS;
+ MPN_RSHIFT (rp, up + kk / GMP_NUMB_BITS, rn, kk % GMP_NUMB_BITS);
+ rn -= rp[rn - 1] == 0;
+
+ /* 9: current buffers: {sp,sn}, {rp,rn} */
+
+ for (c = 0;; c++)
+ {
+ /* Compute S^k in {qp,qn}. */
+ /* W <- S^(k-1) for the next iteration,
+ and S^k = W * S. */
+ wn = mpn_pow_1 (wp, sp, sn, k - 1, qp);
+ mpn_mul (qp, wp, wn, sp, sn);
+ qn = wn + sn;
+ qn -= qp[qn - 1] == 0;
+
+ perf_pow = 1;
+ /* if S^k > floor(U/2^kk), the root approximation was too large */
+ if (qn > rn || (qn == rn && (perf_pow=mpn_cmp (qp, rp, rn)) > 0))
+ MPN_DECR_U (sp, sn, 1);
+ else
+ break;
+ }
+
+ /* 10: current buffers: {sp,sn}, {rp,rn}, {qp,qn}, {wp,wn} */
+
+ /* sometimes two corrections are needed with logbased_root*/
+ ASSERT (c <= 1 + LOGROOT_NEEDS_TWO_CORRECTIONS);
+ ASSERT_ALWAYS (rn >= qn);
+
+ b = sizes[ni - 1] - sizes[ni]; /* number of bits to compute in the
+ next iteration */
+ bn = b / GMP_NUMB_BITS; /* lowest limb from high part of rp[], after shift */
+
+ kk = kk - b;
+ /* nl is the number of limbs in U which contain bits [kk,kk+b-1] */
+ nl = 1 + (kk + b - 1) / GMP_NUMB_BITS - (kk / GMP_NUMB_BITS);
+ /* nl = 1 + floor((kk + b - 1) / GMP_NUMB_BITS)
+ - floor(kk / GMP_NUMB_BITS)
+ <= 1 + (kk + b - 1) / GMP_NUMB_BITS
+ - (kk - GMP_NUMB_BITS + 1) / GMP_NUMB_BITS
+ = 2 + (b - 2) / GMP_NUMB_BITS
+ thus since nl is an integer:
+ nl <= 2 + floor(b/GMP_NUMB_BITS) <= 2 + bn. */
+
+ /* 11: current buffers: {sp,sn}, {rp,rn}, {wp,wn} */
+
+ /* R = R - Q = floor(U/2^kk) - S^k */
+ if (perf_pow != 0)
+ {
+ mpn_sub (rp, rp, rn, qp, qn);
+ MPN_NORMALIZE_NOT_ZERO (rp, rn);
+
+ /* first multiply the remainder by 2^b */
+ MPN_LSHIFT (cy, rp + bn, rp, rn, b % GMP_NUMB_BITS);
+ rn = rn + bn;
+ if (cy != 0)
+ {
+ rp[rn] = cy;
+ rn++;
+ }
+
+ save = rp[bn];
+ /* we have to save rp[bn] up to rp[nl-1], i.e. 1 or 2 limbs */
+ if (nl - 1 > bn)
+ save2 = rp[bn + 1];
+ }
+ else
+ {
+ rn = bn;
+ save2 = save = 0;
+ }
+ /* 2: current buffers: {sp,sn}, {rp,rn}, {wp,wn} */
+
+ /* Now insert bits [kk,kk+b-1] from the input U */
+ MPN_RSHIFT (rp, up + kk / GMP_NUMB_BITS, nl, kk % GMP_NUMB_BITS);
+ /* set to zero high bits of rp[bn] */
+ rp[bn] &= (CNST_LIMB (1) << (b % GMP_NUMB_BITS)) - 1;
+ /* restore corresponding bits */
+ rp[bn] |= save;
+ if (nl - 1 > bn)
+ rp[bn + 1] = save2; /* the low b bits go in rp[0..bn] only, since
+ they start by bit 0 in rp[0], so they use
+ at most ceil(b/GMP_NUMB_BITS) limbs */
+ /* FIXME: Should we normalise {rp,rn} here ?*/
+
+ /* 3: current buffers: {sp,sn}, {rp,rn}, {wp,wn} */
+
+ /* compute {wp, wn} = k * {sp, sn}^(k-1) */
+ cy = mpn_mul_1 (wp, wp, wn, k);
+ wp[wn] = cy;
+ wn += cy != 0;
+
+ /* 6: current buffers: {sp,sn}, {qp,qn} */
+
+ /* multiply the root approximation by 2^b */
+ MPN_LSHIFT (cy, sp + b / GMP_NUMB_BITS, sp, sn, b % GMP_NUMB_BITS);
+ sn = sn + b / GMP_NUMB_BITS;
+ if (cy != 0)
+ {
+ sp[sn] = cy;
+ sn++;
+ }
+
+ save = sp[b / GMP_NUMB_BITS];
+
+ /* Number of limbs used by b bits, when least significant bit is
+ aligned to least limb */
+ bn = (b - 1) / GMP_NUMB_BITS + 1;
+
+ /* 4: current buffers: {sp,sn}, {rp,rn}, {wp,wn} */
+
+ /* now divide {rp, rn} by {wp, wn} to get the low part of the root */
+ if (UNLIKELY (rn < wn))
+ {
+ MPN_FILL (sp, bn, 0);
+ }
+ else
+ {
+ qn = rn - wn; /* expected quotient size */
+ if (qn <= bn) { /* Divide only if result is not too big. */
+ mpn_div_q (qp, rp, rn, wp, wn, scratch);
+ qn += qp[qn] != 0;
+ }
+
+ /* 5: current buffers: {sp,sn}, {qp,qn}.
+ Note: {rp,rn} is not needed any more since we'll compute it from
+ scratch at the end of the loop.
+ */
+
+ /* the quotient should be smaller than 2^b, since the previous
+ approximation was correctly rounded toward zero */
+ if (qn > bn || (qn == bn && (b % GMP_NUMB_BITS != 0) &&
+ qp[qn - 1] >= (CNST_LIMB (1) << (b % GMP_NUMB_BITS))))
+ {
+ for (qn = 1; qn < bn; ++qn)
+ sp[qn - 1] = GMP_NUMB_MAX;
+ sp[qn - 1] = GMP_NUMB_MAX >> (GMP_NUMB_BITS - 1 - ((b - 1) % GMP_NUMB_BITS));
+ }
+ else
+ {
+ /* 7: current buffers: {sp,sn}, {qp,qn} */
+
+ /* Combine sB and q to form sB + q. */
+ MPN_COPY (sp, qp, qn);
+ MPN_ZERO (sp + qn, bn - qn);
+ }
+ }
+ sp[b / GMP_NUMB_BITS] |= save;
+
+ /* 8: current buffer: {sp,sn} */
+
+ }
+
+ /* otherwise we have rn > 0, thus the return value is ok */
+ if (!approx || sp[0] <= CNST_LIMB (1))
+ {
+ for (c = 0;; c++)
+ {
+ /* Compute S^k in {qp,qn}. */
+ /* Last iteration: we don't need W anymore. */
+ /* mpn_pow_1 requires that both qp and wp have enough
+ space to store the result {sp,sn}^k + 1 limb */
+ qn = mpn_pow_1 (qp, sp, sn, k, wp);
+
+ perf_pow = 1;
+ if (qn > un || (qn == un && (perf_pow=mpn_cmp (qp, up, un)) > 0))
+ MPN_DECR_U (sp, sn, 1);
+ else
+ break;
+ };
+
+ /* sometimes two corrections are needed with logbased_root*/
+ ASSERT (c <= 1 + LOGROOT_NEEDS_TWO_CORRECTIONS);
+
+ rn = perf_pow != 0;
+ if (rn != 0 && remp != NULL)
+ {
+ mpn_sub (remp, up, un, qp, qn);
+ rn = un;
+ MPN_NORMALIZE_NOT_ZERO (remp, rn);
+ }
+ }
+
+ TMP_FREE;
+ return rn;
+}