/* 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 A
MUTABLE INTERFACE. 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 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 <stdio.h>
#include <stdlib.h>
#include "gmp.h"
#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)
/* 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] */
static mp_limb_t
mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
{
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 + np0>>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 (np[2 * n - 1] >= GMP_NUMB_HIGHBIT / 2);
if (n == 1)
c = mpn_sqrtrem2 (sp, np, np);
else
{
l = n / 2;
h = n - l;
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 *tp, s0[1], 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;
}
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 mask;
mp_ptr 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;
}