Updated mini-gmp files, from gmp-6.0.0.

ChangeLog

 2014-05-09  Niels Möller  <nisse@lysator.liu.se>
 
+	* mini-gmp.c: Updated, use version from gmp-6.0.0.
+	* mini-gmp.h: Likewise.
+
 	* testsuite/Makefile.in (all): Drop dependency on $(TARGETS), to
 	delay building of test programs until make check.
 

mini-gmp.c

@@ -2,24 +2,33 @@
 
    Contributed to the GNU project by Niels Möller
 
-Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001,
-2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013
-Free Software Foundation, Inc.
+Copyright 1991-1997, 1999-2014 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 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.
+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 Lesser General Public
-License for more details.
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
 
-You should have received a copy of the GNU Lesser General Public License
-along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
+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/.  */
 
 /* NOTE: All functions in this file which are not declared in
    mini-gmp.h are internal, and are not intended to be compatible
@@ -222,11 +231,13 @@ along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
   } while (0)
 #define MPZ_SRCPTR_SWAP(x, y)						\
   do {									\
-    mpz_srcptr __mpz_srcptr_swap__tmp = (x);			\
+    mpz_srcptr __mpz_srcptr_swap__tmp = (x);				\
     (x) = (y);								\
     (y) = __mpz_srcptr_swap__tmp;					\
   } while (0)
 
+const int mp_bits_per_limb = GMP_LIMB_BITS;
+
 
 /* Memory allocation and other helper functions. */
 static void
@@ -342,12 +353,10 @@ mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
 int
 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 {
-  for (; n > 0; n--)
+  while (--n >= 0)
     {
-      if (ap[n-1] < bp[n-1])
-	return -1;
-      else if (ap[n-1] > bp[n-1])
-	return 1;
+      if (ap[n] != bp[n])
+	return ap[n] > bp[n] ? 1 : -1;
     }
   return 0;
 }
@@ -355,10 +364,8 @@ mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 static int
 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
 {
-  if (an > bn)
-    return 1;
-  else if (an < bn)
-    return -1;
+  if (an != bn)
+    return an < bn ? -1 : 1;
   else
     return mpn_cmp (ap, bp, an);
 }
@@ -373,20 +380,31 @@ mpn_normalized_size (mp_srcptr xp, mp_size_t n)
 
 #define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (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;
+}
+
 mp_limb_t
 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 {
   mp_size_t i;
 
   assert (n > 0);
-
-  for (i = 0; i < n; i++)
+  i = 0;
+  do
     {
       mp_limb_t r = ap[i] + b;
       /* Carry out */
       b = (r < b);
       rp[i] = r;
     }
+  while (++i < n);
+
   return b;
 }
 
@@ -429,7 +447,8 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 
   assert (n > 0);
 
-  for (i = 0; i < n; i++)
+  i = 0;
+  do
     {
       mp_limb_t a = ap[i];
       /* Carry out */
@@ -437,6 +456,8 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
       rp[i] = a - b;
       b = cy;
     }
+  while (++i < n);
+
   return b;
 }
 
@@ -602,7 +623,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 - 1; i != 0; i--)
+  for (i = n; --i != 0;)
     {
       low_limb = *--up;
       *--rp = high_limb | (low_limb >> tnc);
@@ -630,7 +651,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 - 1; i != 0; i--)
+  for (i = n; --i != 0;)
     {
       high_limb = *up++;
       *rp++ = low_limb | (high_limb << tnc);
@@ -641,6 +662,46 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
   return retval;
 }
 
+static mp_bitcnt_t
+mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
+		 mp_limb_t ux)
+{
+  unsigned cnt;
+
+  assert (ux == 0 || ux == GMP_LIMB_MAX);
+  assert (0 <= i && i <= un );
+
+  while (limb == 0)
+    {
+      i++;
+      if (i == un)
+	return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
+      limb = ux ^ up[i];
+    }
+  gmp_ctz (cnt, limb);
+  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
+}
+
+mp_bitcnt_t
+mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
+{
+  mp_size_t i;
+  i = bit / GMP_LIMB_BITS;
+
+  return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
+			  i, ptr, i, 0);
+}
+
+mp_bitcnt_t
+mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
+{
+  mp_size_t i;
+  i = bit / GMP_LIMB_BITS;
+
+  return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
+			  i, ptr, i, GMP_LIMB_MAX);
+}
+
 
 /* MPN division interface. */
 mp_limb_t
@@ -715,8 +776,7 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
       if (r < th)
 	{
 	  m--;
-	  if (r > u1 || (r == u1 && tl > u0))
-	    m--;
+	  m -= ((r > u1) | ((r == u1) & (tl > u0)));
 	}
     }
 
@@ -836,14 +896,20 @@ mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
   assert (d > 0);
 
   /* Special case for powers of two. */
-  if (d > 1 && (d & (d-1)) == 0)
+  if ((d & (d-1)) == 0)
     {
-      unsigned shift;
       mp_limb_t r = np[0] & (d-1);
-      gmp_ctz (shift, d);
       if (qp)
-	mpn_rshift (qp, np, nn, shift);
-
+	{
+	  if (d <= 1)
+	    mpn_copyi (qp, np, nn);
+	  else
+	    {
+	      unsigned shift;
+	      gmp_ctz (shift, d);
+	      mpn_rshift (qp, np, nn, shift);
+	    }
+	}
       return r;
     }
   else
@@ -880,7 +946,8 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
 
   r0 = np[nn - 1];
 
-  for (i = nn - 2; i >= 0; i--)
+  i = nn - 2;
+  do
     {
       mp_limb_t n0, q;
       n0 = np[i];
@@ -889,6 +956,7 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
       if (qp)
 	qp[i] = q;
     }
+  while (--i >= 0);
 
   if (shift > 0)
     {
@@ -930,18 +998,19 @@ mpn_div_qr_pi1 (mp_ptr qp,
 
   assert (dn > 2);
   assert (nn >= dn);
-  assert ((dp[dn-1] & GMP_LIMB_HIGHBIT) != 0);
 
   d1 = dp[dn - 1];
   d0 = dp[dn - 2];
 
+  assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
   /* Iteration variable is the index of the q limb.
    *
    * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
    * by            <d1,          d0,        dp[dn-3],  ..., dp[0] >
    */
 
-  for (i = nn - dn; i >= 0; i--)
+  i = nn - dn;
+  do
     {
       mp_limb_t n0 = np[dn-1+i];
 
@@ -973,6 +1042,7 @@ mpn_div_qr_pi1 (mp_ptr qp,
       if (qp)
 	qp[i] = q;
     }
+  while (--i >= 0);
 
   np[dn - 1] = n1;
 }
@@ -994,7 +1064,9 @@ mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
       mp_limb_t nh;
       unsigned shift;
 
-      assert (dp[dn-1] & GMP_LIMB_HIGHBIT);
+      assert (inv->d1 == dp[dn-1]);
+      assert (inv->d0 == dp[dn-2]);
+      assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
 
       shift = inv->shift;
       if (shift > 0)
@@ -1002,9 +1074,6 @@ mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
       else
 	nh = 0;
 
-      assert (inv->d1 == dp[dn-1]);
-      assert (inv->d0 == dp[dn-2]);
-
       mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
 
       if (shift > 0)
@@ -1238,15 +1307,14 @@ mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
 {
   mp_size_t rn;
   mp_limb_t w;
-  unsigned first;
   unsigned k;
   size_t j;
 
-  first = 1 + (sn - 1) % info->exp;
+  k = 1 + (sn - 1) % info->exp;
 
   j = 0;
   w = sp[j++];
-  for (k = 1; k < first; k++)
+  for (; --k > 0; )
     w = w * b + sp[j++];
 
   rp[0] = w;
@@ -1300,7 +1368,7 @@ mpz_init (mpz_t r)
 }
 
 /* The utility of this function is a bit limited, since many functions
-   assings the result variable using mpz_swap. */
+   assigns the result variable using mpz_swap. */
 void
 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
 {
@@ -1422,7 +1490,7 @@ mpz_fits_ulong_p (const mpz_t u)
 {
   mp_size_t us = u->_mp_size;
 
-  return us == 0 || us == 1;
+  return (us == (us > 0));
 }
 
 long int
@@ -1459,6 +1527,48 @@ mpz_getlimbn (const mpz_t u, mp_size_t n)
     return 0;
 }
 
+void
+mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
+{
+  mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
+}
+
+mp_srcptr
+mpz_limbs_read (mpz_srcptr x)
+{
+  return x->_mp_d;;
+}
+
+mp_ptr
+mpz_limbs_modify (mpz_t x, mp_size_t n)
+{
+  assert (n > 0);
+  return MPZ_REALLOC (x, n);
+}
+
+mp_ptr
+mpz_limbs_write (mpz_t x, mp_size_t n)
+{
+  return mpz_limbs_modify (x, n);
+}
+
+void
+mpz_limbs_finish (mpz_t x, mp_size_t xs)
+{
+  mp_size_t xn;
+  xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
+  x->_mp_size = xs < 0 ? -xn : xn;
+}
+
+mpz_srcptr
+mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
+{
+  x->_mp_alloc = 0;
+  x->_mp_d = (mp_ptr) xp;
+  mpz_limbs_finish (x, xs);
+  return x;
+}
+
 
 /* Conversions and comparison to double. */
 void
@@ -1473,19 +1583,15 @@ mpz_set_d (mpz_t r, double x)
 
   /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
      zero or infinity. */
-  if (x == 0.0 || x != x || x == x * 0.5)
+  if (x != x || x == x * 0.5)
     {
       r->_mp_size = 0;
       return;
     }
 
-  if (x < 0.0)
-    {
-      x = - x;
-      sign = 1;
-    }
-  else
-    sign = 0;
+  sign = x < 0.0 ;
+  if (sign)
+    x = - x;
 
   if (x < 1.0)
     {
@@ -1502,8 +1608,9 @@ mpz_set_d (mpz_t r, double x)
   f = (mp_limb_t) x;
   x -= f;
   assert (x < 1.0);
-  rp[rn-1] = f;
-  for (i = rn-1; i-- > 0; )
+  i = rn-1;
+  rp[i] = f;
+  while (--i >= 0)
     {
       x = B * x;
       f = (mp_limb_t) x;
@@ -1611,12 +1718,7 @@ mpz_sgn (const mpz_t u)
 {
   mp_size_t usize = u->_mp_size;
 
-  if (usize > 0)
-    return 1;
-  else if (usize < 0)
-    return -1;
-  else
-    return 0;
+  return (usize > 0) - (usize < 0);
 }
 
 int
@@ -1635,10 +1737,9 @@ mpz_cmp_si (const mpz_t u, long v)
       mp_limb_t ul = u->_mp_d[0];
       if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
 	return -1;
-      else 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 0;
 }
 
 int
@@ -1653,12 +1754,8 @@ mpz_cmp_ui (const mpz_t u, unsigned long v)
   else
     {
       mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
-      if (ul > v)
-	return 1;
-      else if (ul < v)
-	return -1;
+      return (ul > v) - (ul < v);
     }
-  return 0;
 }
 
 int
@@ -1667,16 +1764,12 @@ mpz_cmp (const mpz_t a, const mpz_t b)
   mp_size_t asize = a->_mp_size;
   mp_size_t bsize = b->_mp_size;
 
-  if (asize > bsize)
-    return 1;
-  else if (asize < bsize)
-    return -1;
-  else if (asize > 0)
+  if (asize != bsize)
+    return (asize < bsize) ? -1 : 1;
+  else if (asize >= 0)
     return mpn_cmp (a->_mp_d, b->_mp_d, asize);
-  else if (asize < 0)
-    return -mpn_cmp (a->_mp_d, b->_mp_d, -asize);
   else
-    return 0;
+    return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
 }
 
 int
@@ -1690,12 +1783,7 @@ mpz_cmpabs_ui (const mpz_t u, unsigned long v)
 
   ul = (un == 1) ? u->_mp_d[0] : 0;
 
-  if (ul > v)
-    return 1;
-  else if (ul < v)
-    return -1;
-  else
-    return 0;
+  return (ul > v) - (ul < v);
 }
 
 int
@@ -1753,7 +1841,7 @@ mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
 
   cy = mpn_add_1 (rp, a->_mp_d, an, b);
   rp[an] = cy;
-  an += (cy > 0);
+  an += cy;
 
   return an;
 }
@@ -1815,20 +1903,21 @@ mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
 {
   mp_size_t an = GMP_ABS (a->_mp_size);
   mp_size_t bn = GMP_ABS (b->_mp_size);
-  mp_size_t rn;
   mp_ptr rp;
   mp_limb_t cy;
 
-  rn = GMP_MAX (an, bn);
-  rp = MPZ_REALLOC (r, rn + 1);
-  if (an >= bn)
-    cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
-  else
-    cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an);
+  if (an < bn)
+    {
+      MPZ_SRCPTR_SWAP (a, b);
+      MP_SIZE_T_SWAP (an, bn);
+    }
+
+  rp = MPZ_REALLOC (r, an + 1);
+  cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
 
-  rp[rn] = cy;
+  rp[an] = cy;
 
-  return rn + (cy > 0);
+  return an + cy;
 }
 
 static mp_size_t
@@ -1899,31 +1988,26 @@ mpz_mul_si (mpz_t r, const mpz_t u, long int v)
 void
 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
 {
-  mp_size_t un;
-  mpz_t t;
+  mp_size_t un, us;
   mp_ptr tp;
   mp_limb_t cy;
 
-  un = GMP_ABS (u->_mp_size);
+  us = u->_mp_size;
 
-  if (un == 0 || v == 0)
+  if (us == 0 || v == 0)
     {
       r->_mp_size = 0;
       return;
     }
 
-  mpz_init2 (t, (un + 1) * GMP_LIMB_BITS);
+  un = GMP_ABS (us);
 
-  tp = t->_mp_d;
+  tp = MPZ_REALLOC (r, un + 1);
   cy = mpn_mul_1 (tp, u->_mp_d, un, v);
   tp[un] = cy;
 
-  t->_mp_size = un + (cy > 0);
-  if (u->_mp_size < 0)
-    t->_mp_size = - t->_mp_size;
-
-  mpz_swap (r, t);
-  mpz_clear (t);
+  un += (cy > 0);
+  r->_mp_size = (us < 0) ? - un : un;
 }
 
 void
@@ -1934,8 +2018,8 @@ mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
   mpz_t t;
   mp_ptr tp;
 
-  un = GMP_ABS (u->_mp_size);
-  vn = GMP_ABS (v->_mp_size);
+  un = u->_mp_size;
+  vn = v->_mp_size;
 
   if (un == 0 || vn == 0)
     {
@@ -1943,7 +2027,10 @@ mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
       return;
     }
 
-  sign = (u->_mp_size ^ v->_mp_size) < 0;
+  sign = (un ^ vn) < 0;
+
+  un = GMP_ABS (un);
+  vn = GMP_ABS (vn);
 
   mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
 
@@ -1996,6 +2083,46 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
   r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
 }
 
+void
+mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
+{
+  mpz_t t;
+  mpz_init (t);
+  mpz_mul_ui (t, u, v);
+  mpz_add (r, r, t);
+  mpz_clear (t);
+}
+
+void
+mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
+{
+  mpz_t t;
+  mpz_init (t);
+  mpz_mul_ui (t, u, v);
+  mpz_sub (r, r, t);
+  mpz_clear (t);
+}
+
+void
+mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
+{
+  mpz_t t;
+  mpz_init (t);
+  mpz_mul (t, u, v);
+  mpz_add (r, r, t);
+  mpz_clear (t);
+}
+
+void
+mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
+{
+  mpz_t t;
+  mpz_init (t);
+  mpz_mul (t, u, v);
+  mpz_sub (r, r, t);
+  mpz_clear (t);