From 901f58ce5f4603e14f4cd1bb64d4d5cf06985b89 Mon Sep 17 00:00:00 2001 From: Paul Eggert Date: Sun, 26 Jan 2020 12:52:56 -0800 Subject: [PATCH] Update mini-gmp * src/mini-gmp.c, src/mini-gmp.h: Copy from GMP 6.2.0. This incorporates: 2019-12-05 remove some sizeof(mp_limb_t) 2019-12-04 (mpn_invert_3by2): Remove special code for limb sizes 2019-12-04 (mpn_invert_3by2): Limit size of an intermediate 2019-11-20 (mpn_invert_3by2): Use xor instead of negation 2019-11-19 (mpn_invert_3by2): Move an assert earlier 2019-11-19 (mpn_invert_3by2): Add a new shortcut 2019-11-17 Prepend "unsigned" to MINI_GMP_LIMB_TYPE 2019-11-17 Enable testing with different limb sizes (types) 2019-11-20 Use already defined constants 2019-11-09 Avoid undefined behaviour with small limb sizes --- src/mini-gmp.c | 222 ++++++++++++++++++++++++++----------------------- src/mini-gmp.h | 8 +- 2 files changed, 123 insertions(+), 107 deletions(-) diff --git a/src/mini-gmp.c b/src/mini-gmp.c index bf8a6164981..09c5436be0a 100644 --- a/src/mini-gmp.c +++ b/src/mini-gmp.c @@ -94,11 +94,13 @@ see https://www.gnu.org/licenses/. */ #define gmp_clz(count, x) do { \ mp_limb_t __clz_x = (x); \ - unsigned __clz_c; \ - for (__clz_c = 0; \ - (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ - __clz_c += 8) \ - __clz_x <<= 8; \ + unsigned __clz_c = 0; \ + int LOCAL_SHIFT_BITS = 8; \ + if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) \ + for (; \ + (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ + __clz_c += 8) \ + { __clz_x <<= LOCAL_SHIFT_BITS; } \ for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \ __clz_x <<= 1; \ (count) = __clz_c; \ @@ -143,27 +145,27 @@ see https://www.gnu.org/licenses/. */ w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \ } \ else { \ - mp_limb_t __x0, __x1, __x2, __x3; \ - unsigned __ul, __vl, __uh, __vh; \ - mp_limb_t __u = (u), __v = (v); \ + mp_limb_t __x0, __x1, __x2, __x3; \ + unsigned __ul, __vl, __uh, __vh; \ + mp_limb_t __u = (u), __v = (v); \ \ - __ul = __u & GMP_LLIMB_MASK; \ - __uh = __u >> (GMP_LIMB_BITS / 2); \ - __vl = __v & GMP_LLIMB_MASK; \ - __vh = __v >> (GMP_LIMB_BITS / 2); \ + __ul = __u & GMP_LLIMB_MASK; \ + __uh = __u >> (GMP_LIMB_BITS / 2); \ + __vl = __v & GMP_LLIMB_MASK; \ + __vh = __v >> (GMP_LIMB_BITS / 2); \ \ - __x0 = (mp_limb_t) __ul * __vl; \ - __x1 = (mp_limb_t) __ul * __vh; \ - __x2 = (mp_limb_t) __uh * __vl; \ - __x3 = (mp_limb_t) __uh * __vh; \ + __x0 = (mp_limb_t) __ul * __vl; \ + __x1 = (mp_limb_t) __ul * __vh; \ + __x2 = (mp_limb_t) __uh * __vl; \ + __x3 = (mp_limb_t) __uh * __vh; \ \ - __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ - __x1 += __x2; /* but this indeed can */ \ - if (__x1 < __x2) /* did we get it? */ \ - __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ + __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ + __x1 += __x2; /* but this indeed can */ \ + if (__x1 < __x2) /* did we get it? */ \ + __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ \ - (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \ - (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \ + (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \ + (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \ } \ } while (0) @@ -768,91 +770,81 @@ mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n) mp_limb_t mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) { - int GMP_LIMB_BITS_MUL_3 = GMP_LIMB_BITS * 3; - if (sizeof (unsigned) * CHAR_BIT > GMP_LIMB_BITS * 3) - { - return (((unsigned) 1 << GMP_LIMB_BITS_MUL_3) - 1) / - (((unsigned) u1 << GMP_LIMB_BITS_MUL_3 / 3) + u0); - } - else if (GMP_ULONG_BITS > GMP_LIMB_BITS * 3) - { - return (((unsigned long) 1 << GMP_LIMB_BITS_MUL_3) - 1) / - (((unsigned long) u1 << GMP_LIMB_BITS_MUL_3 / 3) + u0); - } - else { - mp_limb_t r, p, m, ql; - unsigned ul, uh, qh; + mp_limb_t r, m; - assert (u1 >= GMP_LIMB_HIGHBIT); + { + mp_limb_t p, ql; + unsigned ul, uh, qh; - /* 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); + /* 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; + /* 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 ^ GMP_LIMB_MAX) / uh; - /* Adjust to get a half-limb 3/2 inverse, i.e., we want + /* 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), + qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u + = floor( (b (~u) + b-1) / u), - and the remainder + and the remainder - r = b (~u) + b-1 - qh (b uh + ul) + 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. - */ + 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; + r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; - p = (mp_limb_t) qh * ul; - /* Adjustment steps taken from udiv_qrnnd_c */ - if (r < p) - { - qh--; - r += u1; - if (r >= u1) /* i.e. we didn't get carry when adding to r */ - if (r < p) - { - qh--; - r += u1; - } - } - r -= p; + p = (mp_limb_t) qh * ul; + /* Adjustment steps taken from udiv_qrnnd_c */ + if (r < p) + { + qh--; + r += u1; + if (r >= u1) /* i.e. we didn't get carry when adding to r */ + if (r < p) + { + qh--; + r += u1; + } + } + r -= p; - /* Low half of the quotient is + /* 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. */ + 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; + 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 trick, we don't need the high half limb. */ - r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; + /* 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))) - { - ql--; - r += u1; - } - m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; - if (r >= u1) - { - m++; - r -= u1; - } + if (r >= (GMP_LIMB_MAX & (p << (GMP_LIMB_BITS / 2)))) + { + ql--; + r += u1; + } + m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; + if (r >= u1) + { + m++; + r -= u1; + } + } /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a 3/2 inverse. */ @@ -881,7 +873,6 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) } return m; - } } struct gmp_div_inverse @@ -3332,7 +3323,7 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k) mpz_fac_ui (t, k); for (; k > 0; --k) - mpz_mul_ui (r, r, n--); + mpz_mul_ui (r, r, n--); mpz_divexact (r, r, t); mpz_clear (t); @@ -3427,7 +3418,7 @@ gmp_lucas_mod (mpz_t V, mpz_t Qk, long Q, gmp_lucas_step_k_2k (V, Qk, n); /* A step k->k+1 is performed if the bit in $n$ is 1 */ - /* mpz_tstbit(n,bs) or the bit is 0 in $n$ but */ + /* mpz_tstbit(n,bs) or the the bit is 0 in $n$ but */ /* should be 1 in $n+1$ (bs == b0) */ if (b0 == bs || mpz_tstbit (n, bs)) { @@ -3990,13 +3981,18 @@ gmp_popcount_limb (mp_limb_t x) unsigned c; /* Do 16 bits at a time, to avoid limb-sized constants. */ - for (c = 0; x > 0; x >>= 16) + int LOCAL_SHIFT_BITS = 16; + for (c = 0; x > 0;) { unsigned w = x - ((x >> 1) & 0x5555); w = ((w >> 2) & 0x3333) + (w & 0x3333); w = (w >> 4) + w; w = ((w >> 8) & 0x000f) + (w & 0x000f); c += w; + if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) + x >>= LOCAL_SHIFT_BITS; + else + x = 0; } return c; } @@ -4492,7 +4488,7 @@ mpz_export (void *r, size_t *countp, int order, size_t size, int endian, ptrdiff_t word_step; /* The current (partial) limb. */ mp_limb_t limb; - /* The number of bytes left to do in this limb. */ + /* The number of bytes left to to in this limb. */ size_t bytes; /* The index where the limb was read. */ mp_size_t i; @@ -4503,10 +4499,15 @@ mpz_export (void *r, size_t *countp, int order, size_t size, int endian, limb = u->_mp_d[un-1]; assert (limb != 0); - k = 0; - do { - k++; limb >>= CHAR_BIT; - } while (limb != 0); + k = (GMP_LIMB_BITS <= CHAR_BIT); + if (!k) + { + do { + int LOCAL_CHAR_BIT = CHAR_BIT; + k++; limb >>= LOCAL_CHAR_BIT; + } while (limb != 0); + } + /* else limb = 0; */ count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size; @@ -4535,17 +4536,28 @@ mpz_export (void *r, size_t *countp, int order, size_t size, int endian, for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step) { size_t j; - for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + for (j = 0; j < size; ++j, p -= (ptrdiff_t) endian) { - if (bytes == 0) + if (sizeof (mp_limb_t) == 1) { if (i < un) - limb = u->_mp_d[i++]; - bytes = sizeof (mp_limb_t); + *p = u->_mp_d[i++]; + else + *p = 0; + } + else + { + int LOCAL_CHAR_BIT = CHAR_BIT; + if (bytes == 0) + { + if (i < un) + limb = u->_mp_d[i++]; + bytes = sizeof (mp_limb_t); + } + *p = limb; + limb >>= LOCAL_CHAR_BIT; + bytes--; } - *p = limb; - limb >>= CHAR_BIT; - bytes--; } } assert (i == un); diff --git a/src/mini-gmp.h b/src/mini-gmp.h index 27e0c0671a2..7cce3f7a328 100644 --- a/src/mini-gmp.h +++ b/src/mini-gmp.h @@ -1,6 +1,6 @@ /* mini-gmp, a minimalistic implementation of a GNU GMP subset. -Copyright 2011-2015, 2017 Free Software Foundation, Inc. +Copyright 2011-2015, 2017, 2019 Free Software Foundation, Inc. This file is part of the GNU MP Library. @@ -53,7 +53,11 @@ void mp_get_memory_functions (void *(**) (size_t), void *(**) (void *, size_t, size_t), void (**) (void *, size_t)); -typedef unsigned long mp_limb_t; +#ifndef MINI_GMP_LIMB_TYPE +#define MINI_GMP_LIMB_TYPE long +#endif + +typedef unsigned MINI_GMP_LIMB_TYPE mp_limb_t; typedef long mp_size_t; typedef unsigned long mp_bitcnt_t; -- 2.39.2