[Gmp-commit] /var/hg/gmp-6.1: Merge mini-gmp-related changes from main-line
mercurial at gmplib.org
mercurial at gmplib.org
Thu Dec 1 05:21:33 UTC 2016
details: /var/hg/gmp-6.1/rev/af8fcdbda0a6
changeset: 16957:af8fcdbda0a6
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Thu Dec 01 06:21:09 2016 +0100
description:
Merge mini-gmp-related changes from main-line
diffstat:
Makefile.am | 5 +-
mini-gmp/mini-gmp.c | 258 ++++++++++++++++++++++++------------------
mini-gmp/mini-gmp.h | 3 +
mini-gmp/tests/Makefile | 15 +-
mini-gmp/tests/hex-random.c | 71 ++++++++++-
mini-gmp/tests/mini-random.h | 2 +-
mini-gmp/tests/run-tests | 22 ++-
mini-gmp/tests/t-cmp_d.c | 6 +-
mini-gmp/tests/t-cong.c | 2 +-
mini-gmp/tests/t-invert.c | 127 ++++++++++++++-------
mini-gmp/tests/t-limbs.c | 15 +-
mini-gmp/tests/t-logops.c | 2 +-
mini-gmp/tests/t-pprime_p.c | 11 +-
mini-gmp/tests/t-reuse.c | 2 +-
mini-gmp/tests/t-signed.c | 203 ++++++++++++++++++++++-----------
mini-gmp/tests/t-str.c | 18 ++-
mini-gmp/tests/testutils.c | 4 +-
17 files changed, 496 insertions(+), 270 deletions(-)
diffs (truncated from 1465 to 300 lines):
diff -r 53f948a9ded6 -r af8fcdbda0a6 Makefile.am
--- a/Makefile.am Wed Nov 30 10:17:19 2016 +0100
+++ b/Makefile.am Thu Dec 01 06:21:09 2016 +0100
@@ -426,15 +426,14 @@
abs_srcdir="`cd $(srcdir) && pwd`" ; \
$(MKDIR_P) mini-gmp/tests \
&& cd mini-gmp/tests \
- && LD_LIBRARY_PATH="../../.libs:$$LD_LIBRARY_PATH" \
- DYLD_LIBRARY_PATH="../../.libs:$$DYLD_LIBRARY_PATH" \
+ && TEST_LIBRARY_PATH="../../.libs" \
$(MAKE) -f "$$abs_srcdir/mini-gmp/tests/Makefile" \
VPATH="$$abs_srcdir/mini-gmp/tests" \
srcdir="$$abs_srcdir/mini-gmp/tests" \
MINI_GMP_DIR="$$abs_srcdir/mini-gmp" \
LDFLAGS="-L../../.libs" \
LIBS="-lgmp -lm" \
- CC="$(CC_FOR_BUILD)" EXTRA_CFLAGS="-g -I../.." check
+ CC="$(CC)" CFLAGS="$(CFLAGS)" CPPFLAGS="-I../.." check
clean-mini-gmp:
if [ -d mini-gmp/tests ] ; then \
diff -r 53f948a9ded6 -r af8fcdbda0a6 mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c Wed Nov 30 10:17:19 2016 +0100
+++ b/mini-gmp/mini-gmp.c Thu Dec 01 06:21:09 2016 +0100
@@ -2,7 +2,7 @@
Contributed to the GNU project by Niels Möller
-Copyright 1991-1997, 1999-2015 Free Software Foundation, Inc.
+Copyright 1991-1997, 1999-2016 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -69,8 +69,10 @@
#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)))
+
#define gmp_assert_nocarry(x) do { \
- mp_limb_t __cy = x; \
+ mp_limb_t __cy = (x); \
assert (__cy == 0); \
} while (0)
@@ -454,7 +456,7 @@
{
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;
}
@@ -701,24 +703,68 @@
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;
-
- /* 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 */
+ mp_limb_t r, p, m, ql;
+ unsigned ul, uh, qh;
+
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;
@@ -736,11 +782,19 @@
}
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)))
@@ -755,6 +809,8 @@
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;
@@ -1159,7 +1215,7 @@
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;
@@ -1300,6 +1356,8 @@
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)
@@ -1309,6 +1367,8 @@
unsigned k;
size_t j;
+ assert (sn > 0);
+
k = 1 + (sn - 1) % info->exp;
j = 0;
@@ -1318,7 +1378,7 @@
rp[0] = w;
- for (rn = (w > 0); j < sn;)
+ for (rn = 1; j < sn;)
{
mp_limb_t cy;
@@ -1361,9 +1421,11 @@
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
@@ -1384,7 +1446,8 @@
void
mpz_clear (mpz_t r)
{
- gmp_free (r->_mp_d);
+ if (r->_mp_alloc)
+ gmp_free (r->_mp_d);
}
static mp_ptr
@@ -1392,7 +1455,10 @@
{
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)
@@ -1415,7 +1481,7 @@
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);
}
}
@@ -1425,7 +1491,7 @@
if (x > 0)
{
r->_mp_size = 1;
- r->_mp_d[0] = x;
+ MPZ_REALLOC (r, 1)[0] = x;
}
else
r->_mp_size = 0;
@@ -1493,14 +1559,11 @@
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
@@ -1533,7 +1596,7 @@
mp_srcptr
mpz_limbs_read (mpz_srcptr x)
{
- return x->_mp_d;;
+ return x->_mp_d;
}
mp_ptr
@@ -1713,9 +1776,7 @@
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
@@ -1730,13 +1791,7 @@
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;
More information about the gmp-commit
mailing list