From ff8105752ae9d29c50d18e35bc766935642b7476 Mon Sep 17 00:00:00 2001 From: Harsha Sharma Date: Sat, 31 Mar 2018 20:19:41 +0530 Subject: src: Updates for mini-gmp.{c,h} updates from latest stable release of libgmp to get in sync with them Signed-off-by: Harsha Sharma Signed-off-by: Pablo Neira Ayuso --- include/mini-gmp.h | 8 +- src/mini-gmp.c | 436 +++++++++++++++++++++++++++++------------------------ 2 files changed, 245 insertions(+), 199 deletions(-) diff --git a/include/mini-gmp.h b/include/mini-gmp.h index c043ca7e..27e0c067 100644 --- a/include/mini-gmp.h +++ b/include/mini-gmp.h @@ -1,6 +1,6 @@ /* mini-gmp, a minimalistic implementation of a GNU GMP subset. -Copyright 2011-2014 Free Software Foundation, Inc. +Copyright 2011-2015, 2017 Free Software Foundation, Inc. This file is part of the GNU MP Library. @@ -82,6 +82,7 @@ void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t); void mpn_zero (mp_ptr, mp_size_t); int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t); +int mpn_zero_p (mp_srcptr, mp_size_t); mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); @@ -107,6 +108,9 @@ mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t); mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t); +void mpn_com (mp_ptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_neg (mp_ptr, mp_srcptr, mp_size_t); + mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t); mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t); @@ -213,6 +217,8 @@ void mpz_rootrem (mpz_t, mpz_t, const mpz_t, unsigned long); int mpz_root (mpz_t, const mpz_t, unsigned long); void mpz_fac_ui (mpz_t, unsigned long); +void mpz_2fac_ui (mpz_t, unsigned long); +void mpz_mfac_uiui (mpz_t, unsigned long, unsigned long); void mpz_bin_uiui (mpz_t, unsigned long, unsigned long); int mpz_probab_prime_p (const mpz_t, int); diff --git a/src/mini-gmp.c b/src/mini-gmp.c index acbe1bec..360bfa94 100644 --- a/src/mini-gmp.c +++ b/src/mini-gmp.c @@ -2,7 +2,7 @@ Contributed to the GNU project by Niels Möller -Copyright 1991-1997, 1999-2014 Free Software Foundation, Inc. +Copyright 1991-1997, 1999-2018 Free Software Foundation, Inc. This file is part of the GNU MP Library. @@ -69,8 +69,16 @@ see https://www.gnu.org/licenses/. */ #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b)) #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b)) +#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b))) + +/* Return non-zero if xp,xsize and yp,ysize overlap. + If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no + overlap. If both these are false, there's an overlap. */ +#define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize) \ + ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp)) + #define gmp_assert_nocarry(x) do { \ - mp_limb_t __cy = x; \ + mp_limb_t __cy = (x); \ assert (__cy == 0); \ } while (0) @@ -264,12 +272,12 @@ gmp_default_alloc (size_t size) static void * gmp_default_realloc (void *old, size_t old_size, size_t new_size) { - mp_ptr p; + void * p; p = realloc (old, new_size); if (!p) - gmp_die("gmp_default_realoc: Virtual memory exhausted."); + gmp_die("gmp_default_realloc: Virtual memory exhausted."); return p; } @@ -322,14 +330,14 @@ mp_set_memory_functions (void *(*alloc_func) (size_t), static mp_ptr gmp_xalloc_limbs (mp_size_t size) { - return gmp_xalloc (size * sizeof (mp_limb_t)); + return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t)); } static mp_ptr gmp_xrealloc_limbs (mp_ptr old, mp_size_t size) { assert (size > 0); - return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t)); + return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t)); } @@ -346,7 +354,7 @@ mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n) void mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n) { - while (n-- > 0) + while (--n >= 0) d[n] = s[n]; } @@ -373,20 +381,22 @@ mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) static mp_size_t mpn_normalized_size (mp_srcptr xp, mp_size_t n) { - for (; n > 0 && xp[n-1] == 0; n--) - ; + while (n > 0 && xp[n-1] == 0) + --n; return n; } -#define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0) +int +mpn_zero_p(mp_srcptr rp, mp_size_t n) +{ + return mpn_normalized_size (rp, n) == 0; +} void mpn_zero (mp_ptr rp, mp_size_t n) { - mp_size_t i; - - for (i = 0; i < n; i++) - rp[i] = 0; + while (--n >= 0) + rp[n] = 0; } mp_limb_t @@ -452,7 +462,7 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) { mp_limb_t a = ap[i]; /* Carry out */ - mp_limb_t cy = a < b;; + mp_limb_t cy = a < b; rp[i] = a - b; b = cy; } @@ -572,23 +582,24 @@ mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn) { assert (un >= vn); assert (vn >= 1); + assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un)); + assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn)); /* We first multiply by the low order limb. This result can be stored, not added, to rp. We also avoid a loop for zeroing this way. */ rp[un] = mpn_mul_1 (rp, up, un, vp[0]); - rp += 1, vp += 1, vn -= 1; /* Now accumulate the product of up[] and the next higher limb from vp[]. */ - while (vn >= 1) + while (--vn >= 1) { + rp += 1, vp += 1; rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); - rp += 1, vp += 1, vn -= 1; } - return rp[un - 1]; + return rp[un]; } void @@ -608,7 +619,6 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) { mp_limb_t high_limb, low_limb; unsigned int tnc; - mp_size_t i; mp_limb_t retval; assert (n >= 1); @@ -623,7 +633,7 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) retval = low_limb >> tnc; high_limb = (low_limb << cnt); - for (i = n; --i != 0;) + while (--n != 0) { low_limb = *--up; *--rp = high_limb | (low_limb >> tnc); @@ -639,7 +649,6 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) { mp_limb_t high_limb, low_limb; unsigned int tnc; - mp_size_t i; mp_limb_t retval; assert (n >= 1); @@ -651,7 +660,7 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) retval = (high_limb << tnc); low_limb = high_limb >> cnt; - for (i = n; --i != 0;) + while (--n != 0) { high_limb = *up++; *rp++ = low_limb | (high_limb << tnc); @@ -702,24 +711,68 @@ mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit) i, ptr, i, GMP_LIMB_MAX); } +void +mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (--n >= 0) + *rp++ = ~ *up++; +} + +mp_limb_t +mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (*up == 0) + { + *rp = 0; + if (!--n) + return 0; + ++up; ++rp; + } + *rp = - *up; + mpn_com (++rp, ++up, --n); + return 1; +} + /* MPN division interface. */ + +/* The 3/2 inverse is defined as + + m = floor( (B^3-1) / (B u1 + u0)) - B +*/ mp_limb_t mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) { - mp_limb_t r, p, m; - unsigned ul, uh; - unsigned ql, qh; + mp_limb_t r, p, m, ql; + unsigned ul, uh, qh; - /* First, do a 2/1 inverse. */ - /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 < - * B^2 - (B + m) u1 <= u1 */ assert (u1 >= GMP_LIMB_HIGHBIT); + /* For notation, let b denote the half-limb base, so that B = b^2. + Split u1 = b uh + ul. */ ul = u1 & GMP_LLIMB_MASK; uh = u1 >> (GMP_LIMB_BITS / 2); + /* Approximation of the high half of quotient. Differs from the 2/1 + inverse of the half limb uh, since we have already subtracted + u0. */ qh = ~u1 / uh; + + /* Adjust to get a half-limb 3/2 inverse, i.e., we want + + qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u + = floor( (b (~u) + b-1) / u), + + and the remainder + + r = b (~u) + b-1 - qh (b uh + ul) + = b (~u - qh uh) + b-1 - qh ul + + Subtraction of qh ul may underflow, which implies adjustments. + But by normalization, 2 u >= B > qh ul, so we need to adjust by + at most 2. + */ + r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; p = (mp_limb_t) qh * ul; @@ -737,11 +790,19 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) } r -= p; - /* Do a 3/2 division (with half limb size) */ + /* Low half of the quotient is + + ql = floor ( (b r + b-1) / u1). + + This is a 3/2 division (on half-limbs), for which qh is a + suitable inverse. */ + p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; + /* Unlike full-limb 3/2, we can add 1 without overflow. For this to + work, it is essential that ql is a full mp_limb_t. */ ql = (p >> (GMP_LIMB_BITS / 2)) + 1; - /* By the 3/2 method, we don't need the high half limb. */ + /* By the 3/2 trick, we don't need the high half limb. */ r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; if (r >= (p << (GMP_LIMB_BITS / 2))) @@ -756,6 +817,8 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) r -= u1; } + /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a + 3/2 inverse. */ if (u0 > 0) { mp_limb_t th, tl; @@ -867,7 +930,8 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, if (inv->shift > 0) { - tp = gmp_xalloc_limbs (nn); + /* Shift, reusing qp area if possible. In-place shift if qp == np. */ + tp = qp ? qp : gmp_xalloc_limbs (nn); r = mpn_lshift (tp, np, nn, inv->shift); np = tp; } @@ -876,7 +940,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, d = inv->d1; di = inv->di; - while (nn-- > 0) + while (--nn >= 0) { mp_limb_t q; @@ -884,7 +948,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, if (qp) qp[nn] = q; } - if (inv->shift > 0) + if ((inv->shift > 0) && (tp != qp)) gmp_free (tp); return r >> inv->shift; @@ -921,13 +985,12 @@ mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d) } static void -mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, +mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, const struct gmp_div_inverse *inv) { unsigned shift; mp_size_t i; mp_limb_t d1, d0, di, r1, r0; - mp_ptr tp; assert (nn >= 2); shift = inv->shift; @@ -936,11 +999,7 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, di = inv->di; if (shift > 0) - { - tp = gmp_xalloc_limbs (nn); - r1 = mpn_lshift (tp, np, nn, shift); - np = tp; - } + r1 = mpn_lshift (np, np, nn, shift); else r1 = 0; @@ -963,27 +1022,12 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, assert ((r0 << (GMP_LIMB_BITS - shift)) == 0); r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift)); r1 >>= shift; - - gmp_free (tp); } - rp[1] = r1; - rp[0] = r0; + np[1] = r1; + np[0] = r0; } -#if 0 -static void -mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, - mp_limb_t d1, mp_limb_t d0) -{ - struct gmp_div_inverse inv; - assert (nn >= 2); - - mpn_div_qr_2_invert (&inv, d1, d0); - mpn_div_qr_2_preinv (qp, rp, np, nn, &inv); -} -#endif - static void mpn_div_qr_pi1 (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_limb_t n1, @@ -1058,7 +1102,7 @@ mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, if (dn == 1) np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv); else if (dn == 2) - mpn_div_qr_2_preinv (qp, np, np, nn, inv); + mpn_div_qr_2_preinv (qp, np, nn, inv); else { mp_limb_t nh; @@ -1160,7 +1204,7 @@ mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un) unsigned char mask; size_t sn, j; mp_size_t i; - int shift; + unsigned shift; sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]) + bits - 1) / bits; @@ -1301,6 +1345,8 @@ mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn, return rn; } +/* Result is usually normalized, except for all-zero input, in which + case a single zero limb is written at *RP, and 1 is returned. */ static mp_size_t mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn, mp_limb_t b, const struct mpn_base_info *info) @@ -1310,16 +1356,18 @@ mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn, unsigned k; size_t j; + assert (sn > 0); + k = 1 + (sn - 1) % info->exp; j = 0; w = sp[j++]; - for (; --k > 0; ) + while (--k != 0) w = w * b + sp[j++]; rp[0] = w; - for (rn = (w > 0); j < sn;) + for (rn = 1; j < sn;) { mp_limb_t cy; @@ -1362,9 +1410,11 @@ mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base) void mpz_init (mpz_t r) { - r->_mp_alloc = 1; + static const mp_limb_t dummy_limb = 0xc1a0; + + r->_mp_alloc = 0; r->_mp_size = 0; - r->_mp_d = gmp_xalloc_limbs (1); + r->_mp_d = (mp_ptr) &dummy_limb; } /* The utility of this function is a bit limited, since many functions @@ -1385,15 +1435,19 @@ mpz_init2 (mpz_t r, mp_bitcnt_t bits) void mpz_clear (mpz_t r) { - gmp_free (r->_mp_d); + if (r->_mp_alloc) + gmp_free (r->_mp_d); } -static void * +static mp_ptr mpz_realloc (mpz_t r, mp_size_t size) { size = GMP_MAX (size, 1); - r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size); + if (r->_mp_alloc) + r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size); + else + r->_mp_d = gmp_xalloc_limbs (size); r->_mp_alloc = size; if (GMP_ABS (r->_mp_size) > size) @@ -1416,7 +1470,7 @@ mpz_set_si (mpz_t r, signed long int x) else /* (x < 0) */ { r->_mp_size = -1; - r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x); + MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x); } } @@ -1426,7 +1480,7 @@ mpz_set_ui (mpz_t r, unsigned long int x) if (x > 0) { r->_mp_size = 1; - r->_mp_d[0] = x; + MPZ_REALLOC (r, 1)[0] = x; } else r->_mp_size = 0; @@ -1475,14 +1529,12 @@ mpz_fits_slong_p (const mpz_t u) { mp_size_t us = u->_mp_size; - if (us == 0) - return 1; - else if (us == 1) + if (us == 1) return u->_mp_d[0] < GMP_LIMB_HIGHBIT; else if (us == -1) return u->_mp_d[0] <= GMP_LIMB_HIGHBIT; else - return 0; + return (us == 0); } int @@ -1496,14 +1548,11 @@ mpz_fits_ulong_p (const mpz_t u) long int mpz_get_si (const mpz_t u) { - mp_size_t us = u->_mp_size; - - if (us > 0) - return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT); - else if (us < 0) - return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT); + if (u->_mp_size < 0) + /* This expression is necessary to properly handle 0x80000000 */ + return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT); else - return 0; + return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT); } unsigned long int @@ -1536,7 +1585,7 @@ mpz_realloc2 (mpz_t x, mp_bitcnt_t n) mp_srcptr mpz_limbs_read (mpz_srcptr x) { - return x->_mp_d;; + return x->_mp_d; } mp_ptr @@ -1560,11 +1609,19 @@ mpz_limbs_finish (mpz_t x, mp_size_t xs) x->_mp_size = xs < 0 ? -xn : xn; } -mpz_srcptr -mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs) +static mpz_srcptr +mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs) { x->_mp_alloc = 0; x->_mp_d = (mp_ptr) xp; + x->_mp_size = xs; + return x; +} + +mpz_srcptr +mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs) +{ + mpz_roinit_normal_n (x, xp, xs); mpz_limbs_finish (x, xs); return x; } @@ -1716,9 +1773,7 @@ mpz_cmp_d (const mpz_t x, double d) int mpz_sgn (const mpz_t u) { - mp_size_t usize = u->_mp_size; - - return (usize > 0) - (usize < 0); + return GMP_CMP (u->_mp_size, 0); } int @@ -1733,13 +1788,7 @@ mpz_cmp_si (const mpz_t u, long v) else if (usize >= 0) return 1; else /* usize == -1 */ - { - mp_limb_t ul = u->_mp_d[0]; - if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul) - return -1; - else - return (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul; - } + return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), u->_mp_d[0]); } int @@ -1752,10 +1801,7 @@ mpz_cmp_ui (const mpz_t u, unsigned long v) else if (usize < 0) return -1; else - { - mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0; - return (ul > v) - (ul < v); - } + return GMP_CMP (mpz_get_ui (u), v); } int @@ -1775,15 +1821,10 @@ mpz_cmp (const mpz_t a, const mpz_t b) int mpz_cmpabs_ui (const mpz_t u, unsigned long v) { - mp_size_t un = GMP_ABS (u->_mp_size); - mp_limb_t ul; - - if (un > 1) + if (GMP_ABS (u->_mp_size) > 1) return 1; - - ul = (un == 1) ? u->_mp_d[0] : 0; - - return (ul > v) - (ul < v); + else + return GMP_CMP (mpz_get_ui (u), v); } int @@ -1796,18 +1837,14 @@ mpz_cmpabs (const mpz_t u, const mpz_t v) void mpz_abs (mpz_t r, const mpz_t u) { - if (r != u) - mpz_set (r, u); - + mpz_set (r, u); r->_mp_size = GMP_ABS (r->_mp_size); } void mpz_neg (mpz_t r, const mpz_t u) { - if (r != u) - mpz_set (r, u); - + mpz_set (r, u); r->_mp_size = -r->_mp_size; } @@ -1833,7 +1870,7 @@ mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b) an = GMP_ABS (a->_mp_size); if (an == 0) { - r->_mp_d[0] = b; + MPZ_REALLOC (r, 1)[0] = b; return b > 0; } @@ -1852,14 +1889,15 @@ static mp_size_t mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b) { mp_size_t an = GMP_ABS (a->_mp_size); - mp_ptr rp = MPZ_REALLOC (r, an); + mp_ptr rp; if (an == 0) { - rp[0] = b; + MPZ_REALLOC (r, 1)[0] = b; return -(b > 0); } - else if (an == 1 && a->_mp_d[0] < b) + rp = MPZ_REALLOC (r, an); + if (an == 1 && a->_mp_d[0] < b) { rp[0] = b - a->_mp_d[0]; return -1; @@ -2077,8 +2115,7 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits) else mpn_copyd (rp + limbs, u->_mp_d, un); - while (limbs > 0) - rp[--limbs] = 0; + mpn_zero (rp, limbs); r->_mp_size = (u->_mp_size < 0) ? - rn : rn; } @@ -2331,7 +2368,6 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index, if (qn <= 0) qn = 0; - else { qp = MPZ_REALLOC (q, qn); @@ -2385,16 +2421,9 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, { /* Have to negate and sign extend. */ mp_size_t i; - mp_limb_t cy; - for (cy = 1, i = 0; i < un; i++) - { - mp_limb_t s = ~u->_mp_d[i] + cy; - cy = s < cy; - rp[i] = s; - } - assert (cy == 0); - for (; i < rn - 1; i++) + gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un)); + for (i = un; i < rn - 1; i++) rp[i] = GMP_LIMB_MAX; rp[rn-1] = mask; @@ -2419,23 +2448,13 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ { /* If r != 0, compute 2^{bit_count} - r. */ - mp_size_t i; - - for (i = 0; i < rn && rp[i] == 0; i++) - ; - if (i < rn) - { - /* r > 0, need to flip sign. */ - rp[i] = ~rp[i] + 1; - while (++i < rn) - rp[i] = ~rp[i]; + mpn_neg (rp, rp, rn); - rp[rn-1] &= mask; + rp[rn-1] &= mask; - /* us is not used for anything else, so we can modify it - here to indicate flipped sign. */ - us = -us; - } + /* us is not used for anything else, so we can modify it + here to indicate flipped sign. */ + us = -us; } } rn = mpn_normalized_size (rp, rn); @@ -2550,7 +2569,7 @@ mpz_div_qr_ui (mpz_t q, mpz_t r, if (r) { - r->_mp_d[0] = rl; + MPZ_REALLOC (r, 1)[0] = rl; r->_mp_size = rs; } if (q) @@ -3074,9 +3093,7 @@ void mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e) { mpz_t b; - mpz_init_set_ui (b, blimb); - mpz_pow_ui (r, b, e); - mpz_clear (b); + mpz_pow_ui (r, mpz_roinit_normal_n (b, &blimb, blimb != 0), e); } void @@ -3148,7 +3165,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m) } mpz_init_set_ui (tr, 1); - while (en-- > 0) + while (--en >= 0) { mp_limb_t w = e->_mp_d[en]; mp_limb_t bit; @@ -3188,9 +3205,7 @@ void mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m) { mpz_t e; - mpz_init_set_ui (e, elimb); - mpz_powm (r, b, e, m); - mpz_clear (e); + mpz_powm (r, b, mpz_roinit_normal_n (e, &elimb, elimb != 0), m); } /* x=trunc(y^(1/z)), r=y-x^z */ @@ -3215,12 +3230,8 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z) } mpz_init (u); - { - mp_bitcnt_t tb; - tb = mpz_sizeinbase (y, 2) / z + 1; - mpz_init2 (t, tb); - mpz_setbit (t, tb); - } + mpz_init (t); + mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1); if (z == 2) /* simplify sqrt loop: z-1 == 1 */ do { @@ -3301,7 +3312,7 @@ mpn_perfect_square_p (mp_srcptr p, mp_size_t n) assert (n > 0); assert (p [n-1] != 0); - return mpz_root (NULL, mpz_roinit_n (t, p, n), 2); + return mpz_root (NULL, mpz_roinit_normal_n (t, p, n), 2); } mp_size_t @@ -3315,7 +3326,7 @@ mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n) mpz_init (r); mpz_init (s); - mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2); + mpz_rootrem (s, r, mpz_roinit_normal_n (u, p, n), 2); assert (s->_mp_size == (n+1)/2); mpn_copyd (sp, s->_mp_d, s->_mp_size); @@ -3330,11 +3341,24 @@ mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n) /* Combinatorics */ void -mpz_fac_ui (mpz_t x, unsigned long n) +mpz_mfac_uiui (mpz_t x, unsigned long n, unsigned long m) { mpz_set_ui (x, n + (n == 0)); - for (;n > 2;) - mpz_mul_ui (x, x, --n); + if (m + 1 < 2) return; + while (n > m + 1) + mpz_mul_ui (x, x, n -= m); +} + +void +mpz_2fac_ui (mpz_t x, unsigned long n) +{ + mpz_mfac_uiui (x, n, 2); +} + +void +mpz_fac_ui (mpz_t x, unsigned long n) +{ + mpz_mfac_uiui (x, n, 1); } void @@ -3350,7 +3374,7 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k) mpz_init (t); mpz_fac_ui (t, k); - for (; k > 0; k--) + for (; k > 0; --k) mpz_mul_ui (r, r, n--); mpz_divexact (r, r, t); @@ -3504,7 +3528,7 @@ mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index) must be complemented. */ if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0) return bit ^ 1; - while (limb_index-- > 0) + while (--limb_index >= 0) if (d->_mp_d[limb_index] > 0) return bit ^ 1; } @@ -3647,7 +3671,7 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v) /* If the smaller input is positive, higher limbs don't matter. */ rn = vx ? un : vn; - rp = MPZ_REALLOC (r, rn + rc); + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); up = u->_mp_d; vp = v->_mp_d; @@ -3720,7 +3744,7 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v) don't matter. */ rn = vx ? vn : un; - rp = MPZ_REALLOC (r, rn + rc); + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); up = u->_mp_d; vp = v->_mp_d; @@ -3789,7 +3813,7 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v) vx = -vc; rx = -rc; - rp = MPZ_REALLOC (r, un + rc); + rp = MPZ_REALLOC (r, un + (mp_size_t) rc); up = u->_mp_d; vp = v->_mp_d; @@ -3835,10 +3859,10 @@ gmp_popcount_limb (mp_limb_t x) /* Do 16 bits at a time, to avoid limb-sized constants. */ for (c = 0; x > 0; x >>= 16) { - unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555); + unsigned w = x - ((x >> 1) & 0x5555); w = ((w >> 2) & 0x3333) + (w & 0x3333); - w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f); - w = (w >> 8) + (w & 0x00ff); + w = (w >> 4) + w; + w = ((w >> 8) & 0x000f) + (w & 0x000f); c += w; } return c; @@ -3999,7 +4023,7 @@ mpz_sizeinbase (const mpz_t u, int base) size_t ndigits; assert (base >= 2); - assert (base <= 36); + assert (base <= 62); un = GMP_ABS (u->_mp_size); if (un == 0) @@ -4049,23 +4073,26 @@ mpz_get_str (char *sp, int base, const mpz_t u) mp_size_t un; size_t i, sn; - if (base >= 0) + digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"; + if (base > 1) { - digits = "0123456789abcdefghijklmnopqrstuvwxyz"; + if (base <= 36) + digits = "0123456789abcdefghijklmnopqrstuvwxyz"; + else if (base > 62) + return NULL; } + else if (base >= -1) + base = 10; else { base = -base; - digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + if (base > 36) + return NULL; } - if (base <= 1) - base = 10; - if (base > 36) - return NULL; sn = 1 + mpz_sizeinbase (u, base); if (!sp) - sp = gmp_xalloc (1 + sn); + sp = (char *) gmp_xalloc (1 + sn); un = GMP_ABS (u->_mp_size); @@ -4109,14 +4136,14 @@ mpz_get_str (char *sp, int base, const mpz_t u) int mpz_set_str (mpz_t r, const char *sp, int base) { - unsigned bits; + unsigned bits, value_of_a; mp_size_t rn, alloc; mp_ptr rp; - size_t sn; + size_t dn; int sign; unsigned char *dp; - assert (base == 0 || (base >= 2 && base <= 36)); + assert (base == 0 || (base >= 2 && base <= 62)); while (isspace( (unsigned char) *sp)) sp++; @@ -4126,18 +4153,17 @@ mpz_set_str (mpz_t r, const char *sp, int base) if (base == 0) { - if (*sp == '0') + if (sp[0] == '0') { - sp++; - if (*sp == 'x' || *sp == 'X') + if (sp[1] == 'x' || sp[1] == 'X') { base = 16; - sp++; + sp += 2; } - else if (*sp == 'b' || *sp == 'B') + else if (sp[1] == 'b' || sp[1] == 'B') { base = 2; - sp++; + sp += 2; } else base = 8; @@ -4146,49 +4172,63 @@ mpz_set_str (mpz_t r, const char *sp, int base) base = 10; } - sn = strlen (sp); - dp = gmp_xalloc (sn + (sn == 0)); + if (!*sp) + { + r->_mp_size = 0; + return -1; + } + dp = (unsigned char *) gmp_xalloc (strlen (sp)); - for (sn = 0; *sp; sp++) + value_of_a = (base > 36) ? 36 : 10; + for (dn = 0; *sp; sp++) { unsigned digit; if (isspace ((unsigned char) *sp)) continue; - if (*sp >= '0' && *sp <= '9') + else if (*sp >= '0' && *sp <= '9') digit = *sp - '0'; else if (*sp >= 'a' && *sp <= 'z') - digit = *sp - 'a' + 10; + digit = *sp - 'a' + value_of_a; else if (*sp >= 'A' && *sp <= 'Z') digit = *sp - 'A' + 10; else digit = base; /* fail */ - if (digit >= base) + if (digit >= (unsigned) base) { gmp_free (dp); r->_mp_size = 0; return -1; } - dp[sn++] = digit; + dp[dn++] = digit; } + if (!dn) + { + gmp_free (dp); + r->_mp_size = 0; + return -1; + } bits = mpn_base_power_of_two_p (base); if (bits > 0) { - alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; rp = MPZ_REALLOC (r, alloc); - rn = mpn_set_str_bits (rp, dp, sn, bits); + rn = mpn_set_str_bits (rp, dp, dn, bits); } else { struct mpn_base_info info; mpn_get_base_info (&info, base); - alloc = (sn + info.exp - 1) / info.exp; + alloc = (dn + info.exp - 1) / info.exp; rp = MPZ_REALLOC (r, alloc); - rn = mpn_set_str_other (rp, dp, sn, base, &info); + rn = mpn_set_str_other (rp, dp, dn, base, &info); + /* Normalization, needed for all-zero input. */ + assert (rn > 0); + rn -= rp[rn-1] == 0; } assert (rn <= alloc); gmp_free (dp); -- cgit v1.2.3