From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/mpn/generic/sqrtrem.c | 555 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 555 insertions(+) create mode 100644 gmp-6.3.0/mpn/generic/sqrtrem.c (limited to 'gmp-6.3.0/mpn/generic/sqrtrem.c') diff --git a/gmp-6.3.0/mpn/generic/sqrtrem.c b/gmp-6.3.0/mpn/generic/sqrtrem.c new file mode 100644 index 0000000..cc6dd9c --- /dev/null +++ b/gmp-6.3.0/mpn/generic/sqrtrem.c @@ -0,0 +1,555 @@ +/* mpn_sqrtrem -- square root and remainder + + Contributed to the GNU project by Paul Zimmermann (most code), + Torbjorn Granlund (mpn_sqrtrem1) and Marco Bodrato (mpn_dc_sqrt). + + THE FUNCTIONS IN THIS FILE EXCEPT mpn_sqrtrem ARE INTERNAL WITH MUTABLE + INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. + IN FACT, IT IS ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A + FUTURE GMP RELEASE. + +Copyright 1999-2002, 2004, 2005, 2008, 2010, 2012, 2015, 2017 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/. */ + + +/* See "Karatsuba Square Root", reference in gmp.texi. */ + + +#include +#include + +#include "gmp-impl.h" +#include "longlong.h" +#define USE_DIVAPPR_Q 1 +#define TRACE(x) + +static const unsigned char invsqrttab[384] = /* The common 0x100 was removed */ +{ + 0xff,0xfd,0xfb,0xf9,0xf7,0xf5,0xf3,0xf2, /* sqrt(1/80)..sqrt(1/87) */ + 0xf0,0xee,0xec,0xea,0xe9,0xe7,0xe5,0xe4, /* sqrt(1/88)..sqrt(1/8f) */ + 0xe2,0xe0,0xdf,0xdd,0xdb,0xda,0xd8,0xd7, /* sqrt(1/90)..sqrt(1/97) */ + 0xd5,0xd4,0xd2,0xd1,0xcf,0xce,0xcc,0xcb, /* sqrt(1/98)..sqrt(1/9f) */ + 0xc9,0xc8,0xc6,0xc5,0xc4,0xc2,0xc1,0xc0, /* sqrt(1/a0)..sqrt(1/a7) */ + 0xbe,0xbd,0xbc,0xba,0xb9,0xb8,0xb7,0xb5, /* sqrt(1/a8)..sqrt(1/af) */ + 0xb4,0xb3,0xb2,0xb0,0xaf,0xae,0xad,0xac, /* sqrt(1/b0)..sqrt(1/b7) */ + 0xaa,0xa9,0xa8,0xa7,0xa6,0xa5,0xa4,0xa3, /* sqrt(1/b8)..sqrt(1/bf) */ + 0xa2,0xa0,0x9f,0x9e,0x9d,0x9c,0x9b,0x9a, /* sqrt(1/c0)..sqrt(1/c7) */ + 0x99,0x98,0x97,0x96,0x95,0x94,0x93,0x92, /* sqrt(1/c8)..sqrt(1/cf) */ + 0x91,0x90,0x8f,0x8e,0x8d,0x8c,0x8c,0x8b, /* sqrt(1/d0)..sqrt(1/d7) */ + 0x8a,0x89,0x88,0x87,0x86,0x85,0x84,0x83, /* sqrt(1/d8)..sqrt(1/df) */ + 0x83,0x82,0x81,0x80,0x7f,0x7e,0x7e,0x7d, /* sqrt(1/e0)..sqrt(1/e7) */ + 0x7c,0x7b,0x7a,0x79,0x79,0x78,0x77,0x76, /* sqrt(1/e8)..sqrt(1/ef) */ + 0x76,0x75,0x74,0x73,0x72,0x72,0x71,0x70, /* sqrt(1/f0)..sqrt(1/f7) */ + 0x6f,0x6f,0x6e,0x6d,0x6d,0x6c,0x6b,0x6a, /* sqrt(1/f8)..sqrt(1/ff) */ + 0x6a,0x69,0x68,0x68,0x67,0x66,0x66,0x65, /* sqrt(1/100)..sqrt(1/107) */ + 0x64,0x64,0x63,0x62,0x62,0x61,0x60,0x60, /* sqrt(1/108)..sqrt(1/10f) */ + 0x5f,0x5e,0x5e,0x5d,0x5c,0x5c,0x5b,0x5a, /* sqrt(1/110)..sqrt(1/117) */ + 0x5a,0x59,0x59,0x58,0x57,0x57,0x56,0x56, /* sqrt(1/118)..sqrt(1/11f) */ + 0x55,0x54,0x54,0x53,0x53,0x52,0x52,0x51, /* sqrt(1/120)..sqrt(1/127) */ + 0x50,0x50,0x4f,0x4f,0x4e,0x4e,0x4d,0x4d, /* sqrt(1/128)..sqrt(1/12f) */ + 0x4c,0x4b,0x4b,0x4a,0x4a,0x49,0x49,0x48, /* sqrt(1/130)..sqrt(1/137) */ + 0x48,0x47,0x47,0x46,0x46,0x45,0x45,0x44, /* sqrt(1/138)..sqrt(1/13f) */ + 0x44,0x43,0x43,0x42,0x42,0x41,0x41,0x40, /* sqrt(1/140)..sqrt(1/147) */ + 0x40,0x3f,0x3f,0x3e,0x3e,0x3d,0x3d,0x3c, /* sqrt(1/148)..sqrt(1/14f) */ + 0x3c,0x3b,0x3b,0x3a,0x3a,0x39,0x39,0x39, /* sqrt(1/150)..sqrt(1/157) */ + 0x38,0x38,0x37,0x37,0x36,0x36,0x35,0x35, /* sqrt(1/158)..sqrt(1/15f) */ + 0x35,0x34,0x34,0x33,0x33,0x32,0x32,0x32, /* sqrt(1/160)..sqrt(1/167) */ + 0x31,0x31,0x30,0x30,0x2f,0x2f,0x2f,0x2e, /* sqrt(1/168)..sqrt(1/16f) */ + 0x2e,0x2d,0x2d,0x2d,0x2c,0x2c,0x2b,0x2b, /* sqrt(1/170)..sqrt(1/177) */ + 0x2b,0x2a,0x2a,0x29,0x29,0x29,0x28,0x28, /* sqrt(1/178)..sqrt(1/17f) */ + 0x27,0x27,0x27,0x26,0x26,0x26,0x25,0x25, /* sqrt(1/180)..sqrt(1/187) */ + 0x24,0x24,0x24,0x23,0x23,0x23,0x22,0x22, /* sqrt(1/188)..sqrt(1/18f) */ + 0x21,0x21,0x21,0x20,0x20,0x20,0x1f,0x1f, /* sqrt(1/190)..sqrt(1/197) */ + 0x1f,0x1e,0x1e,0x1e,0x1d,0x1d,0x1d,0x1c, /* sqrt(1/198)..sqrt(1/19f) */ + 0x1c,0x1b,0x1b,0x1b,0x1a,0x1a,0x1a,0x19, /* sqrt(1/1a0)..sqrt(1/1a7) */ + 0x19,0x19,0x18,0x18,0x18,0x18,0x17,0x17, /* sqrt(1/1a8)..sqrt(1/1af) */ + 0x17,0x16,0x16,0x16,0x15,0x15,0x15,0x14, /* sqrt(1/1b0)..sqrt(1/1b7) */ + 0x14,0x14,0x13,0x13,0x13,0x12,0x12,0x12, /* sqrt(1/1b8)..sqrt(1/1bf) */ + 0x12,0x11,0x11,0x11,0x10,0x10,0x10,0x0f, /* sqrt(1/1c0)..sqrt(1/1c7) */ + 0x0f,0x0f,0x0f,0x0e,0x0e,0x0e,0x0d,0x0d, /* sqrt(1/1c8)..sqrt(1/1cf) */ + 0x0d,0x0c,0x0c,0x0c,0x0c,0x0b,0x0b,0x0b, /* sqrt(1/1d0)..sqrt(1/1d7) */ + 0x0a,0x0a,0x0a,0x0a,0x09,0x09,0x09,0x09, /* sqrt(1/1d8)..sqrt(1/1df) */ + 0x08,0x08,0x08,0x07,0x07,0x07,0x07,0x06, /* sqrt(1/1e0)..sqrt(1/1e7) */ + 0x06,0x06,0x06,0x05,0x05,0x05,0x04,0x04, /* sqrt(1/1e8)..sqrt(1/1ef) */ + 0x04,0x04,0x03,0x03,0x03,0x03,0x02,0x02, /* sqrt(1/1f0)..sqrt(1/1f7) */ + 0x02,0x02,0x01,0x01,0x01,0x01,0x00,0x00 /* sqrt(1/1f8)..sqrt(1/1ff) */ +}; + +/* Compute s = floor(sqrt(a0)), and *rp = a0 - s^2. */ + +#if GMP_NUMB_BITS > 32 +#define MAGIC CNST_LIMB(0x10000000000) /* 0xffe7debbfc < MAGIC < 0x232b1850f410 */ +#else +#define MAGIC CNST_LIMB(0x100000) /* 0xfee6f < MAGIC < 0x29cbc8 */ +#endif + +static mp_limb_t +mpn_sqrtrem1 (mp_ptr rp, mp_limb_t a0) +{ +#if GMP_NUMB_BITS > 32 + mp_limb_t a1; +#endif + mp_limb_t x0, t2, t, x2; + unsigned abits; + + ASSERT_ALWAYS (GMP_NAIL_BITS == 0); + ASSERT_ALWAYS (GMP_LIMB_BITS == 32 || GMP_LIMB_BITS == 64); + ASSERT (a0 >= GMP_NUMB_HIGHBIT / 2); + + /* Use Newton iterations for approximating 1/sqrt(a) instead of sqrt(a), + since we can do the former without division. As part of the last + iteration convert from 1/sqrt(a) to sqrt(a). */ + + abits = a0 >> (GMP_LIMB_BITS - 1 - 8); /* extract bits for table lookup */ + x0 = 0x100 | invsqrttab[abits - 0x80]; /* initial 1/sqrt(a) */ + + /* x0 is now an 8 bits approximation of 1/sqrt(a0) */ + +#if GMP_NUMB_BITS > 32 + a1 = a0 >> (GMP_LIMB_BITS - 1 - 32); + t = (mp_limb_signed_t) (CNST_LIMB(0x2000000000000) - 0x30000 - a1 * x0 * x0) >> 16; + x0 = (x0 << 16) + ((mp_limb_signed_t) (x0 * t) >> (16+2)); + + /* x0 is now a 16 bits approximation of 1/sqrt(a0) */ + + t2 = x0 * (a0 >> (32-8)); + t = t2 >> 25; + t = ((mp_limb_signed_t) ((a0 << 14) - t * t - MAGIC) >> (32-8)); + x0 = t2 + ((mp_limb_signed_t) (x0 * t) >> 15); + x0 >>= 32; +#else + t2 = x0 * (a0 >> (16-8)); + t = t2 >> 13; + t = ((mp_limb_signed_t) ((a0 << 6) - t * t - MAGIC) >> (16-8)); + x0 = t2 + ((mp_limb_signed_t) (x0 * t) >> 7); + x0 >>= 16; +#endif + + /* x0 is now a full limb approximation of sqrt(a0) */ + + x2 = x0 * x0; + if (x2 + 2*x0 <= a0 - 1) + { + x2 += 2*x0 + 1; + x0++; + } + + *rp = a0 - x2; + return x0; +} + + +#define Prec (GMP_NUMB_BITS >> 1) +#if ! defined(SQRTREM2_INPLACE) +#define SQRTREM2_INPLACE 0 +#endif + +/* same as mpn_sqrtrem, but for size=2 and {np, 2} normalized + return cc such that {np, 2} = sp[0]^2 + cc*2^GMP_NUMB_BITS + rp[0] */ +#if SQRTREM2_INPLACE +#define CALL_SQRTREM2_INPLACE(sp,rp) mpn_sqrtrem2 (sp, rp) +static mp_limb_t +mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp) +{ + mp_srcptr np = rp; +#else +#define CALL_SQRTREM2_INPLACE(sp,rp) mpn_sqrtrem2 (sp, rp, rp) +static mp_limb_t +mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np) +{ +#endif + mp_limb_t q, u, np0, sp0, rp0, q2; + int cc; + + ASSERT (np[1] >= GMP_NUMB_HIGHBIT / 2); + + np0 = np[0]; + sp0 = mpn_sqrtrem1 (rp, np[1]); + rp0 = rp[0]; + /* rp0 <= 2*sp0 < 2^(Prec + 1) */ + rp0 = (rp0 << (Prec - 1)) + (np0 >> (Prec + 1)); + q = rp0 / sp0; + /* q <= 2^Prec, if q = 2^Prec, reduce the overestimate. */ + q -= q >> Prec; + /* now we have q < 2^Prec */ + u = rp0 - q * sp0; + /* now we have (rp[0]<>Prec)/2 = q * sp0 + u */ + sp0 = (sp0 << Prec) | q; + cc = u >> (Prec - 1); + rp0 = ((u << (Prec + 1)) & GMP_NUMB_MASK) + (np0 & ((CNST_LIMB (1) << (Prec + 1)) - 1)); + /* subtract q * q from rp */ + q2 = q * q; + cc -= rp0 < q2; + rp0 -= q2; + if (cc < 0) + { + rp0 += sp0; + cc += rp0 < sp0; + --sp0; + rp0 += sp0; + cc += rp0 < sp0; + } + + rp[0] = rp0; + sp[0] = sp0; + return cc; +} + +/* writes in {sp, n} the square root (rounded towards zero) of {np, 2n}, + and in {np, n} the low n limbs of the remainder, returns the high + limb of the remainder (which is 0 or 1). + Assumes {np, 2n} is normalized, i.e. np[2n-1] >= B/4 + where B=2^GMP_NUMB_BITS. + Needs a scratch of n/2+1 limbs. */ +static mp_limb_t +mpn_dc_sqrtrem (mp_ptr sp, mp_ptr np, mp_size_t n, mp_limb_t approx, mp_ptr scratch) +{ + mp_limb_t q; /* carry out of {sp, n} */ + int c, b; /* carry out of remainder */ + mp_size_t l, h; + + ASSERT (n > 1); + ASSERT (np[2 * n - 1] >= GMP_NUMB_HIGHBIT / 2); + + l = n / 2; + h = n - l; + if (h == 1) + q = CALL_SQRTREM2_INPLACE (sp + l, np + 2 * l); + else + q = mpn_dc_sqrtrem (sp + l, np + 2 * l, h, 0, scratch); + if (q != 0) + ASSERT_CARRY (mpn_sub_n (np + 2 * l, np + 2 * l, sp + l, h)); + TRACE(printf("tdiv_qr(,,,,%u,,%u) -> %u\n", (unsigned) n, (unsigned) h, (unsigned) (n - h + 1))); + mpn_tdiv_qr (scratch, np + l, 0, np + l, n, sp + l, h); + q += scratch[l]; + c = scratch[0] & 1; + mpn_rshift (sp, scratch, l, 1); + sp[l - 1] |= (q << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK; + if (UNLIKELY ((sp[0] & approx) != 0)) /* (sp[0] & mask) > 1 */ + return 1; /* Remainder is non-zero */ + q >>= 1; + if (c != 0) + c = mpn_add_n (np + l, np + l, sp + l, h); + TRACE(printf("sqr(,,%u)\n", (unsigned) l)); + mpn_sqr (np + n, sp, l); + b = q + mpn_sub_n (np, np, np + n, 2 * l); + c -= (l == h) ? b : mpn_sub_1 (np + 2 * l, np + 2 * l, 1, (mp_limb_t) b); + + if (c < 0) + { + q = mpn_add_1 (sp + l, sp + l, h, q); +#if HAVE_NATIVE_mpn_addlsh1_n_ip1 || HAVE_NATIVE_mpn_addlsh1_n + c += mpn_addlsh1_n_ip1 (np, sp, n) + 2 * q; +#else + c += mpn_addmul_1 (np, sp, n, CNST_LIMB(2)) + 2 * q; +#endif + c -= mpn_sub_1 (np, np, n, CNST_LIMB(1)); + q -= mpn_sub_1 (sp, sp, n, CNST_LIMB(1)); + } + + return c; +} + +#if USE_DIVAPPR_Q +static void +mpn_divappr_q (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, mp_ptr scratch) +{ + gmp_pi1_t inv; + mp_limb_t qh; + ASSERT (dn > 2); + ASSERT (nn >= dn); + ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0); + + MPN_COPY (scratch, np, nn); + invert_pi1 (inv, dp[dn-1], dp[dn-2]); + if (BELOW_THRESHOLD (dn, DC_DIVAPPR_Q_THRESHOLD)) + qh = mpn_sbpi1_divappr_q (qp, scratch, nn, dp, dn, inv.inv32); + else if (BELOW_THRESHOLD (dn, MU_DIVAPPR_Q_THRESHOLD)) + qh = mpn_dcpi1_divappr_q (qp, scratch, nn, dp, dn, &inv); + else + { + mp_size_t itch = mpn_mu_divappr_q_itch (nn, dn, 0); + TMP_DECL; + TMP_MARK; + /* Sadly, scratch is too small. */ + qh = mpn_mu_divappr_q (qp, np, nn, dp, dn, TMP_ALLOC_LIMBS (itch)); + TMP_FREE; + } + qp [nn - dn] = qh; +} +#endif + +/* writes in {sp, n} the square root (rounded towards zero) of {np, 2n-odd}, + returns zero if the operand was a perfect square, one otherwise. + Assumes {np, 2n-odd}*4^nsh is normalized, i.e. B > np[2n-1-odd]*4^nsh >= B/4 + where B=2^GMP_NUMB_BITS. + THINK: In the odd case, three more (dummy) limbs are taken into account, + when nsh is maximal, two limbs are discarded from the result of the + division. Too much? Is a single dummy limb enough? */ +static int +mpn_dc_sqrt (mp_ptr sp, mp_srcptr np, mp_size_t n, unsigned nsh, unsigned odd) +{ + mp_limb_t q; /* carry out of {sp, n} */ + int c; /* carry out of remainder */ + mp_size_t l, h; + mp_ptr qp, tp, scratch; + TMP_DECL; + TMP_MARK; + + ASSERT (np[2 * n - 1 - odd] != 0); + ASSERT (n > 4); + ASSERT (nsh < GMP_NUMB_BITS / 2); + + l = (n - 1) / 2; + h = n - l; + ASSERT (n >= l + 2 && l + 2 >= h && h > l && l >= 1 + odd); + scratch = TMP_ALLOC_LIMBS (l + 2 * n + 5 - USE_DIVAPPR_Q); /* n + 2-USE_DIVAPPR_Q */ + tp = scratch + n + 2 - USE_DIVAPPR_Q; /* n + h + 1, but tp [-1] is writable */ + if (nsh != 0) + { + /* o is used to exactly set the lowest bits of the dividend, is it needed? */ + int o = l > (1 + odd); + ASSERT_NOCARRY (mpn_lshift (tp - o, np + l - 1 - o - odd, n + h + 1 + o, 2 * nsh)); + } + else + MPN_COPY (tp, np + l - 1 - odd, n + h + 1); + q = mpn_dc_sqrtrem (sp + l, tp + l + 1, h, 0, scratch); + if (q != 0) + ASSERT_CARRY (mpn_sub_n (tp + l + 1, tp + l + 1, sp + l, h)); + qp = tp + n + 1; /* l + 2 */ + TRACE(printf("div(appr)_q(,,%u,,%u) -> %u \n", (unsigned) n+1, (unsigned) h, (unsigned) (n + 1 - h + 1))); +#if USE_DIVAPPR_Q + mpn_divappr_q (qp, tp, n + 1, sp + l, h, scratch); +#else + mpn_div_q (qp, tp, n + 1, sp + l, h, scratch); +#endif + q += qp [l + 1]; + c = 1; + if (q > 1) + { + /* FIXME: if s!=0 we will shift later, a noop on this area. */ + MPN_FILL (sp, l, GMP_NUMB_MAX); + } + else + { + /* FIXME: if s!=0 we will shift again later, shift just once. */ + mpn_rshift (sp, qp + 1, l, 1); + sp[l - 1] |= q << (GMP_NUMB_BITS - 1); + if (((qp[0] >> (2 + USE_DIVAPPR_Q)) | /* < 3 + 4*USE_DIVAPPR_Q */ + (qp[1] & (GMP_NUMB_MASK >> ((GMP_NUMB_BITS >> odd)- nsh - 1)))) == 0) + { + mp_limb_t cy; + /* Approximation is not good enough, the extra limb(+ nsh bits) + is smaller than needed to absorb the possible error. */ + /* {qp + 1, l + 1} equals 2*{sp, l} */ + /* FIXME: use mullo or wrap-around, or directly evaluate + remainder with a single sqrmod_bnm1. */ + TRACE(printf("mul(,,%u,,%u)\n", (unsigned) h, (unsigned) (l+1))); + ASSERT_NOCARRY (mpn_mul (scratch, sp + l, h, qp + 1, l + 1)); + /* Compute the remainder of the previous mpn_div(appr)_q. */ + cy = mpn_sub_n (tp + 1, tp + 1, scratch, h); +#if USE_DIVAPPR_Q || WANT_ASSERT + MPN_DECR_U (tp + 1 + h, l, cy); +#if USE_DIVAPPR_Q + ASSERT (mpn_cmp (tp + 1 + h, scratch + h, l) <= 0); + if (mpn_cmp (tp + 1 + h, scratch + h, l) < 0) + { + /* May happen only if div result was not exact. */ +#if HAVE_NATIVE_mpn_addlsh1_n_ip1 || HAVE_NATIVE_mpn_addlsh1_n + cy = mpn_addlsh1_n_ip1 (tp + 1, sp + l, h); +#else + cy = mpn_addmul_1 (tp + 1, sp + l, h, CNST_LIMB(2)); +#endif + ASSERT_NOCARRY (mpn_add_1 (tp + 1 + h, tp + 1 + h, l, cy)); + MPN_DECR_U (sp, l, 1); + } + /* Can the root be exact when a correction was needed? We + did not find an example, but it depends on divappr + internals, and we can not assume it true in general...*/ + /* else */ +#else /* WANT_ASSERT */ + ASSERT (mpn_cmp (tp + 1 + h, scratch + h, l) == 0); +#endif +#endif + if (mpn_zero_p (tp + l + 1, h - l)) + { + TRACE(printf("sqr(,,%u)\n", (unsigned) l)); + mpn_sqr (scratch, sp, l); + c = mpn_cmp (tp + 1, scratch + l, l); + if (c == 0) + { + if (nsh != 0) + { + mpn_lshift (tp, np, l, 2 * nsh); + np = tp; + } + c = mpn_cmp (np, scratch + odd, l - odd); + } + if (c < 0) + { + MPN_DECR_U (sp, l, 1); + c = 1; + } + } + } + } + TMP_FREE; + + if ((odd | nsh) != 0) + mpn_rshift (sp, sp, n, nsh + (odd ? GMP_NUMB_BITS / 2 : 0)); + return c; +} + + +mp_size_t +mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn) +{ + mp_limb_t cc, high, rl; + int c; + mp_size_t rn, tn; + TMP_DECL; + + ASSERT (nn > 0); + ASSERT_MPN (np, nn); + + ASSERT (np[nn - 1] != 0); + ASSERT (rp == NULL || MPN_SAME_OR_SEPARATE_P (np, rp, nn)); + ASSERT (rp == NULL || ! MPN_OVERLAP_P (sp, (nn + 1) / 2, rp, nn)); + ASSERT (! MPN_OVERLAP_P (sp, (nn + 1) / 2, np, nn)); + + high = np[nn - 1]; + if (high & (GMP_NUMB_HIGHBIT | (GMP_NUMB_HIGHBIT / 2))) + c = 0; + else + { + count_leading_zeros (c, high); + c -= GMP_NAIL_BITS; + + c = c / 2; /* we have to shift left by 2c bits to normalize {np, nn} */ + } + if (nn == 1) { + if (c == 0) + { + sp[0] = mpn_sqrtrem1 (&rl, high); + if (rp != NULL) + rp[0] = rl; + } + else + { + cc = mpn_sqrtrem1 (&rl, high << (2*c)) >> c; + sp[0] = cc; + if (rp != NULL) + rp[0] = rl = high - cc*cc; + } + return rl != 0; + } + if (nn == 2) { + mp_limb_t tp [2]; + if (rp == NULL) rp = tp; + if (c == 0) + { +#if SQRTREM2_INPLACE + rp[1] = high; + rp[0] = np[0]; + cc = CALL_SQRTREM2_INPLACE (sp, rp); +#else + cc = mpn_sqrtrem2 (sp, rp, np); +#endif + rp[1] = cc; + return ((rp[0] | cc) != 0) + cc; + } + else + { + rl = np[0]; + rp[1] = (high << (2*c)) | (rl >> (GMP_NUMB_BITS - 2*c)); + rp[0] = rl << (2*c); + CALL_SQRTREM2_INPLACE (sp, rp); + cc = sp[0] >>= c; /* c != 0, the highest bit of the root cc is 0. */ + rp[0] = rl -= cc*cc; /* Computed modulo 2^GMP_LIMB_BITS, because it's smaller. */ + return rl != 0; + } + } + tn = (nn + 1) / 2; /* 2*tn is the smallest even integer >= nn */ + + if ((rp == NULL) && (nn > 8)) + return mpn_dc_sqrt (sp, np, tn, c, nn & 1); + TMP_MARK; + if (((nn & 1) | c) != 0) + { + mp_limb_t s0[1], mask; + mp_ptr tp, scratch; + TMP_ALLOC_LIMBS_2 (tp, 2 * tn, scratch, tn / 2 + 1); + tp[0] = 0; /* needed only when 2*tn > nn, but saves a test */ + if (c != 0) + mpn_lshift (tp + (nn & 1), np, nn, 2 * c); + else + MPN_COPY (tp + (nn & 1), np, nn); + c += (nn & 1) ? GMP_NUMB_BITS / 2 : 0; /* c now represents k */ + mask = (CNST_LIMB (1) << c) - 1; + rl = mpn_dc_sqrtrem (sp, tp, tn, (rp == NULL) ? mask - 1 : 0, scratch); + /* We have 2^(2k)*N = S^2 + R where k = c + (2tn-nn)*GMP_NUMB_BITS/2, + thus 2^(2k)*N = (S-s0)^2 + 2*S*s0 - s0^2 + R where s0=S mod 2^k */ + s0[0] = sp[0] & mask; /* S mod 2^k */ + rl += mpn_addmul_1 (tp, sp, tn, 2 * s0[0]); /* R = R + 2*s0*S */ + cc = mpn_submul_1 (tp, s0, 1, s0[0]); + rl -= (tn > 1) ? mpn_sub_1 (tp + 1, tp + 1, tn - 1, cc) : cc; + mpn_rshift (sp, sp, tn, c); + tp[tn] = rl; + if (rp == NULL) + rp = tp; + c = c << 1; + if (c < GMP_NUMB_BITS) + tn++; + else + { + tp++; + c -= GMP_NUMB_BITS; + } + if (c != 0) + mpn_rshift (rp, tp, tn, c); + else + MPN_COPY_INCR (rp, tp, tn); + rn = tn; + } + else + { + if (rp != np) + { + if (rp == NULL) /* nn <= 8 */ + rp = TMP_SALLOC_LIMBS (nn); + MPN_COPY (rp, np, nn); + } + rn = tn + (rp[tn] = mpn_dc_sqrtrem (sp, rp, tn, 0, TMP_ALLOC_LIMBS(tn / 2 + 1))); + } + + MPN_NORMALIZE (rp, rn); + + TMP_FREE; + return rn; +} -- cgit v1.2.3