mirror of
https://review.haiku-os.org/buildtools
synced 2025-02-12 08:47:41 +01:00
4c74403188
Change-Id: Iae65e95cfa0d92091b8b0a424ae36d88efa76aa9 Reviewed-on: https://review.haiku-os.org/c/buildtools/+/3020 Reviewed-by: Adrien Destugues <pulkomandy@gmail.com>
556 lines
18 KiB
C
556 lines
18 KiB
C
/* 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 <stdio.h>
|
|
#include <stdlib.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)
|
|
#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 + 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 (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;
|
|
}
|