mirror
/
openwrt-packages
mirror of https://git.openwrt.org/feed/packages.git
1
0
Fork 0
Code Releases Activity
openwrt-packages/libs/nettle/patches/0001-Updated-mini-gmp.patch

1739 lines
42 KiB
Diff
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

From 57700c26d73cf7fa6f5cfaec1145eccf388acab9 Mon Sep 17 00:00:00 2001
From: Nikos Mavrogiannopoulos <nmav@gnutls.org>
Date: Sun, 9 Mar 2014 11:27:42 +0100
Subject: [PATCH 1/5] Updated mini-gmp
---
mini-gmp.c | 890 +++++++++++++++++++++++++++++++++++++++----------------------
mini-gmp.h | 55 +++-
2 files changed, 618 insertions(+), 327 deletions(-)
diff --git a/mini-gmp.c b/mini-gmp.c
index 8b6f070..766df30 100644
--- a/mini-gmp.c
+++ b/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[rn] = cy;
+ rp = MPZ_REALLOC (r, an + 1);
+ cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
+
+ 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);
+}
+
/* MPZ division */
enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
@@ -2060,8 +2187,7 @@ mpz_div_qr (mpz_t q, mpz_t r,
mp_size_t qn, rn;
mpz_t tq, tr;
- mpz_init (tr);
- mpz_set (tr, n);
+ mpz_init_set (tr, n);
np = tr->_mp_d;
qn = nn - dn + 1;
@@ -2171,10 +2297,7 @@ mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
void
mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
{
- if (d->_mp_size >= 0)
- mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
- else
- mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
+ mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
}
static void
@@ -2184,7 +2307,7 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
mp_size_t un, qn;
mp_size_t limb_cnt;
mp_ptr qp;
- mp_limb_t adjust;
+ int adjust;
un = u->_mp_size;
if (un == 0)
@@ -2226,7 +2349,8 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
q->_mp_size = qn;
- mpz_add_ui (q, q, adjust);
+ if (adjust)
+ mpz_add_ui (q, q, 1);
if (un < 0)
mpz_neg (q, q);
}
@@ -2303,7 +2427,7 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
{
/* r > 0, need to flip sign. */
rp[i] = ~rp[i] + 1;
- for (i++; i < rn; i++)
+ while (++i < rn)
rp[i] = ~rp[i];
rp[rn-1] &= mask;
@@ -2366,6 +2490,24 @@ mpz_divisible_p (const mpz_t n, const mpz_t d)
return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
}
+int
+mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
+{
+ mpz_t t;
+ int res;
+
+ /* a == b (mod 0) iff a == b */
+ if (mpz_sgn (m) == 0)
+ return (mpz_cmp (a, b) == 0);
+
+ mpz_init (t);
+ mpz_sub (t, a, b);
+ res = mpz_divisible_p (t, m);
+ mpz_clear (t);
+
+ return res;
+}
+
static unsigned long
mpz_div_qr_ui (mpz_t q, mpz_t r,
const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
@@ -2579,32 +2721,16 @@ mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
}
static mp_bitcnt_t
-mpz_make_odd (mpz_t r, const mpz_t u)
+mpz_make_odd (mpz_t r)
{
- mp_size_t un, rn, i;
- mp_ptr rp;
- unsigned shift;
-
- un = GMP_ABS (u->_mp_size);
- assert (un > 0);
+ mp_bitcnt_t shift;
- for (i = 0; u->_mp_d[i] == 0; i++)
- ;
-
- gmp_ctz (shift, u->_mp_d[i]);
-
- rn = un - i;
- rp = MPZ_REALLOC (r, rn);
- if (shift > 0)
- {
- mpn_rshift (rp, u->_mp_d + i, rn, shift);
- rn -= (rp[rn-1] == 0);
- }
- else
- mpn_copyi (rp, u->_mp_d + i, rn);
+ assert (r->_mp_size > 0);
+ /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
+ shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
+ mpz_tdiv_q_2exp (r, r, shift);
- r->_mp_size = rn;
- return i * GMP_LIMB_BITS + shift;
+ return shift;
}
void
@@ -2627,8 +2753,10 @@ mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
mpz_init (tu);
mpz_init (tv);
- uz = mpz_make_odd (tu, u);
- vz = mpz_make_odd (tv, v);
+ mpz_abs (tu, u);
+ uz = mpz_make_odd (tu);
+ mpz_abs (tv, v);
+ vz = mpz_make_odd (tv);
gz = GMP_MIN (uz, vz);
if (tu->_mp_size < tv->_mp_size)
@@ -2644,7 +2772,7 @@ mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
{
int c;
- mpz_make_odd (tu, tu);
+ mpz_make_odd (tu);
c = mpz_cmp (tu, tv);
if (c == 0)
{
@@ -2706,8 +2834,10 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
mpz_init (t0);
mpz_init (t1);
- uz = mpz_make_odd (tu, u);
- vz = mpz_make_odd (tv, v);
+ mpz_abs (tu, u);
+ uz = mpz_make_odd (tu);
+ mpz_abs (tv, v);
+ vz = mpz_make_odd (tv);
gz = GMP_MIN (uz, vz);
uz -= gz;
@@ -2755,7 +2885,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
if (tu->_mp_size > 0)
{
mp_bitcnt_t shift;
- shift = mpz_make_odd (tu, tu);
+ shift = mpz_make_odd (tu);
mpz_mul_2exp (t0, t0, shift);
mpz_mul_2exp (s0, s0, shift);
power += shift;
@@ -2778,7 +2908,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
mpz_add (t0, t0, t1);
mpz_add (s0, s0, s1);
- shift = mpz_make_odd (tv, tv);
+ shift = mpz_make_odd (tv);
mpz_mul_2exp (t1, t1, shift);
mpz_mul_2exp (s1, s1, shift);
}
@@ -2788,7 +2918,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
mpz_add (t1, t0, t1);
mpz_add (s1, s0, s1);
- shift = mpz_make_odd (tu, tu);
+ shift = mpz_make_odd (tu);
mpz_mul_2exp (t0, t0, shift);
mpz_mul_2exp (s0, s0, shift);
}
@@ -2926,12 +3056,16 @@ mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
mpz_t tr;
mpz_init_set_ui (tr, 1);
- for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1)
+ bit = GMP_ULONG_HIGHBIT;
+ do
{
mpz_mul (tr, tr, tr);
if (e & bit)
mpz_mul (tr, tr, b);
+ bit >>= 1;
}
+ while (bit > 0);
+
mpz_swap (r, tr);
mpz_clear (tr);
}
@@ -2987,7 +3121,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
if (e->_mp_size < 0)
{
if (!mpz_invert (base, b, m))
- gmp_die ("mpz_powm: Negative exponent and non-invertibe base.");
+ gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
}
else
{
@@ -3019,7 +3153,8 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
mp_limb_t w = e->_mp_d[en];
mp_limb_t bit;
- for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1)
+ bit = GMP_LIMB_HIGHBIT;
+ do
{
mpz_mul (tr, tr, tr);
if (w & bit)
@@ -3029,7 +3164,9 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
}
+ bit >>= 1;
}
+ while (bit > 0);
}
/* Final reduction */
@@ -3064,21 +3201,26 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
mpz_t t, u;
sgn = y->_mp_size < 0;
- if (sgn && (z & 1) == 0)
+ if ((~z & sgn) != 0)
gmp_die ("mpz_rootrem: Negative argument, with even root.");
if (z == 0)
gmp_die ("mpz_rootrem: Zeroth root.");
if (mpz_cmpabs_ui (y, 1) <= 0) {
- mpz_set (x, y);
+ if (x)
+ mpz_set (x, y);
if (r)
r->_mp_size = 0;
return;
}
- mpz_init (t);
mpz_init (u);
- mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
+ {
+ mp_bitcnt_t tb;
+ tb = mpz_sizeinbase (y, 2) / z + 1;
+ mpz_init2 (t, tb);
+ mpz_setbit (t, tb);
+ }
if (z == 2) /* simplify sqrt loop: z-1 == 1 */
do {
@@ -3110,7 +3252,8 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
mpz_pow_ui (t, u, z);
mpz_sub (r, y, t);
}
- mpz_swap (x, u);
+ if (x)
+ mpz_swap (x, u);
mpz_clear (u);
mpz_clear (t);
}
@@ -3142,19 +3285,56 @@ mpz_sqrt (mpz_t s, const mpz_t u)
mpz_rootrem (s, NULL, u, 2);
}
+int
+mpz_perfect_square_p (const mpz_t u)
+{
+ if (u->_mp_size <= 0)
+ return (u->_mp_size == 0);
+ else
+ return mpz_root (NULL, u, 2);
+}
+
+int
+mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
+{
+ mpz_t t;
+
+ assert (n > 0);
+ assert (p [n-1] != 0);
+ return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
+}
+
+mp_size_t
+mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
+{
+ mpz_t s, r, u;
+ mp_size_t res;
+
+ assert (n > 0);
+ assert (p [n-1] != 0);
+
+ mpz_init (r);
+ mpz_init (s);
+ mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
+
+ assert (s->_mp_size == (n+1)/2);
+ mpn_copyd (sp, s->_mp_d, s->_mp_size);
+ mpz_clear (s);
+ res = r->_mp_size;
+ if (rp)
+ mpn_copyd (rp, r->_mp_d, res);
+ mpz_clear (r);
+ return res;
+}
/* Combinatorics */
void
mpz_fac_ui (mpz_t x, unsigned long n)
{
- if (n < 2) {
- mpz_set_ui (x, 1);
- return;
- }
- mpz_set_ui (x, n);
- for (;--n > 1;)
- mpz_mul_ui (x, x, n);
+ mpz_set_ui (x, n + (n == 0));
+ for (;n > 2;)
+ mpz_mul_ui (x, x, --n);
}
void
@@ -3162,25 +3342,120 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
{
mpz_t t;
- if (k > n) {
- r->_mp_size = 0;
- return;
- }
- mpz_fac_ui (r, n);
+ mpz_set_ui (r, k <= n);
+
+ if (k > (n >> 1))
+ k = (k <= n) ? n - k : 0;
+
mpz_init (t);
mpz_fac_ui (t, k);
- mpz_divexact (r, r, t);
- mpz_fac_ui (t, n - k);
+
+ for (; k > 0; k--)
+ mpz_mul_ui (r, r, n--);
+
mpz_divexact (r, r, t);
mpz_clear (t);
}
+/* Primality testing */
+static int
+gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
+ const mpz_t q, mp_bitcnt_t k)
+{
+ mp_bitcnt_t i;
+
+ /* Caller must initialize y to the base. */
+ mpz_powm (y, y, q, n);
+
+ if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
+ return 1;
+
+ for (i = 1; i < k; i++)
+ {
+ mpz_powm_ui (y, y, 2, n);
+ if (mpz_cmp (y, nm1) == 0)
+ return 1;
+ if (mpz_cmp_ui (y, 1) == 0)
+ return 0;
+ }
+ return 0;
+}
+
+/* This product is 0xc0cfd797, and fits in 32 bits. */
+#define GMP_PRIME_PRODUCT \
+ (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
+
+/* Bit (p+1)/2 is set, for each odd prime <= 61 */
+#define GMP_PRIME_MASK 0xc96996dcUL
+
+int
+mpz_probab_prime_p (const mpz_t n, int reps)
+{
+ mpz_t nm1;
+ mpz_t q;
+ mpz_t y;
+ mp_bitcnt_t k;
+ int is_prime;
+ int j;
+
+ /* Note that we use the absolute value of n only, for compatibility
+ with the real GMP. */
+ if (mpz_even_p (n))
+ return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
+
+ /* Above test excludes n == 0 */
+ assert (n->_mp_size != 0);
+
+ if (mpz_cmpabs_ui (n, 64) < 0)
+ return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
+
+ if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
+ return 0;
+
+ /* All prime factors are >= 31. */
+ if (mpz_cmpabs_ui (n, 31*31) < 0)
+ return 2;
+
+ /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
+ j^2 + j + 41 using Euler's polynomial. We potentially stop early,
+ if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
+ 30 (a[30] == 971 > 31*31 == 961). */
+
+ mpz_init (nm1);
+ mpz_init (q);
+ mpz_init (y);
+
+ /* Find q and k, where q is odd and n = 1 + 2**k * q. */
+ mpz_abs (nm1, n);
+ mpz_sub_ui (nm1, nm1, 1);
+ k = mpz_scan1 (nm1, 0);
+ mpz_tdiv_q_2exp (q, nm1, k);
+
+ for (j = 0, is_prime = 1; is_prime && j < reps; j++)
+ {
+ mpz_set_ui (y, (unsigned long) j*j+j+41);
+ if (mpz_cmp (y, nm1) >= 0)
+ {
+ /* Don't try any further bases. */
+ assert (j >= 30);
+ break;
+ }
+ is_prime &= gmp_millerrabin (n, nm1, y, q, k);
+ }
+ mpz_clear (nm1);
+ mpz_clear (q);
+ mpz_clear (y);
+
+ return is_prime;
+}
+
+
/* Logical operations and bit manipulation. */
/* Numbers are treated as if represented in two's complement (and
infinitely sign extended). For a negative values we get the two's
- complement from -x = ~x + 1, where ~ is bitwise complementt.
+ complement from -x = ~x + 1, where ~ is bitwise complement.
Negation transforms
xxxx10...0
@@ -3374,7 +3649,8 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
up = u->_mp_d;
vp = v->_mp_d;
- for (i = 0; i < vn; i++)
+ i = 0;
+ do
{
ul = (up[i] ^ ux) + uc;
uc = ul < uc;
@@ -3386,6 +3662,7 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
rc = rl < rc;
rp[i] = rl;
}
+ while (++i < vn);
assert (vc == 0);
for (; i < rn; i++)
@@ -3445,7 +3722,8 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
up = u->_mp_d;
vp = v->_mp_d;
- for (i = 0; i < vn; i++)
+ i = 0;
+ do
{
ul = (up[i] ^ ux) + uc;
uc = ul < uc;
@@ -3457,6 +3735,7 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
rc = rl < rc;
rp[i] = rl;
}
+ while (++i < vn);
assert (vc == 0);
for (; i < rn; i++)
@@ -3512,7 +3791,8 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
up = u->_mp_d;
vp = v->_mp_d;
- for (i = 0; i < vn; i++)
+ i = 0;
+ do
{
ul = (up[i] ^ ux) + uc;
uc = ul < uc;
@@ -3524,6 +3804,7 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
rc = rl < rc;
rp[i] = rl;
}
+ while (++i < vn);
assert (vc == 0);
for (; i < un; i++)
@@ -3561,20 +3842,28 @@ gmp_popcount_limb (mp_limb_t x)
}
mp_bitcnt_t
-mpz_popcount (const mpz_t u)
+mpn_popcount (mp_srcptr p, mp_size_t n)
{
- mp_size_t un, i;
+ mp_size_t i;
mp_bitcnt_t c;
+ for (c = 0, i = 0; i < n; i++)
+ c += gmp_popcount_limb (p[i]);
+
+ return c;
+}
+
+mp_bitcnt_t
+mpz_popcount (const mpz_t u)
+{
+ mp_size_t un;
+
un = u->_mp_size;
if (un < 0)
return ~(mp_bitcnt_t) 0;
- for (c = 0, i = 0; i < un; i++)
- c += gmp_popcount_limb (u->_mp_d[i]);
-
- return c;
+ return mpn_popcount (u->_mp_d, un);
}
mp_bitcnt_t
@@ -3591,16 +3880,13 @@ mpz_hamdist (const mpz_t u, const mpz_t v)
if ( (un ^ vn) < 0)
return ~(mp_bitcnt_t) 0;
- if (un < 0)
+ comp = - (uc = vc = (un < 0));
+ if (uc)
{
assert (vn < 0);
un = -un;
vn = -vn;
- uc = vc = 1;
- comp = - (mp_limb_t) 1;
}
- else
- uc = vc = comp = 0;
up = u->_mp_d;
vp = v->_mp_d;
@@ -3636,10 +3922,8 @@ mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
{
mp_ptr up;
mp_size_t us, un, i;
- mp_limb_t limb, ux, uc;
- unsigned cnt;
+ mp_limb_t limb, ux;
- up = u->_mp_d;
us = u->_mp_size;
un = GMP_ABS (us);
i = starting_bit / GMP_LIMB_BITS;
@@ -3649,36 +3933,24 @@ mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
if (i >= un)
return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
- if (us < 0)
- {
- ux = GMP_LIMB_MAX;
- uc = mpn_zero_p (up, i);
- }
- else
- ux = uc = 0;
-
- limb = (ux ^ up[i]) + uc;
- uc = limb < uc;
-
- /* Mask to 0 all bits before starting_bit, thus ignoring them. */
- limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
+ up = u->_mp_d;
+ ux = 0;
+ limb = up[i];
- while (limb == 0)
+ if (starting_bit != 0)
{
- i++;
- if (i == un)
+ if (us < 0)
{
- assert (uc == 0);
- /* For the u > 0 case, this can happen only for the first
- masked limb. For the u < 0 case, it happens when the
- highest limbs of the absolute value are all ones. */
- return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
+ ux = mpn_zero_p (up, i);
+ limb = ~ limb + ux;
+ ux = - (mp_limb_t) (limb >= ux);
}
- limb = (ux ^ up[i]) + uc;
- uc = limb < uc;
+
+ /* Mask to 0 all bits before starting_bit, thus ignoring them. */
+ limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
}
- gmp_ctz (cnt, limb);
- return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
+
+ return mpn_common_scan (limb, i, up, un, ux);
}
mp_bitcnt_t
@@ -3686,46 +3958,28 @@ mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
{
mp_ptr up;
mp_size_t us, un, i;
- mp_limb_t limb, ux, uc;
- unsigned cnt;
+ mp_limb_t limb, ux;
- up = u->_mp_d;
us = u->_mp_size;
+ ux = - (mp_limb_t) (us >= 0);
un = GMP_ABS (us);
i = starting_bit / GMP_LIMB_BITS;
/* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
u<0. Notice this test picks up all cases of u==0 too. */
if (i >= un)
- return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0);
+ return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
- if (us < 0)
- {
- ux = GMP_LIMB_MAX;
- uc = mpn_zero_p (up, i);
- }
- else
- ux = uc = 0;
+ up = u->_mp_d;
+ limb = up[i] ^ ux;
- limb = (ux ^ up[i]) + uc;
- uc = limb < uc;
+ if (ux == 0)
+ limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
- /* Mask to 1 all bits before starting_bit, thus ignoring them. */
- limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1;
+ /* Mask all bits before starting_bit, thus ignoring them. */
+ limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
- while (limb == GMP_LIMB_MAX)
- {
- i++;
- if (i == un)
- {
- assert (uc == 0);
- return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0);
- }
- limb = (ux ^ up[i]) + uc;
- uc = limb < uc;
- }
- gmp_ctz (cnt, ~limb);
- return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
+ return mpn_common_scan (limb, i, up, un, ux);
}
@@ -3771,11 +4025,15 @@ mpz_sizeinbase (const mpz_t u, int base)
mpn_copyi (tp, up, un);
mpn_div_qr_1_invert (&bi, base);
- for (ndigits = 0; un > 0; ndigits++)
+ ndigits = 0;
+ do
{
+ ndigits++;
mpn_div_qr_1_preinv (tp, tp, un, &bi);
un -= (tp[un-1] == 0);
}
+ while (un > 0);
+
gmp_free (tp);
return ndigits;
}
@@ -3852,7 +4110,6 @@ mpz_set_str (mpz_t r, const char *sp, int base)
mp_size_t rn, alloc;
mp_ptr rp;
size_t sn;
- size_t dn;
int sign;
unsigned char *dp;
@@ -3861,13 +4118,8 @@ mpz_set_str (mpz_t r, const char *sp, int base)
while (isspace( (unsigned char) *sp))
sp++;
- if (*sp == '-')
- {
- sign = 1;
- sp++;
- }
- else
- sign = 0;
+ sign = (*sp == '-');
+ sp += sign;
if (base == 0)
{
@@ -3894,7 +4146,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
sn = strlen (sp);
dp = gmp_xalloc (sn + (sn == 0));
- for (dn = 0; *sp; sp++)
+ for (sn = 0; *sp; sp++)
{
unsigned digit;
@@ -3916,7 +4168,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
return -1;
}
- dp[dn++] = digit;
+ dp[sn++] = digit;
}
bits = mpn_base_power_of_two_p (base);
@@ -3925,7 +4177,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
{
alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
rp = MPZ_REALLOC (r, alloc);
- rn = mpn_set_str_bits (rp, dp, dn, bits);
+ rn = mpn_set_str_bits (rp, dp, sn, bits);
}
else
{
@@ -3933,7 +4185,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
mpn_get_base_info (&info, base);
alloc = (sn + info.exp - 1) / info.exp;
rp = MPZ_REALLOC (r, alloc);
- rn = mpn_set_str_other (rp, dp, dn, base, &info);
+ rn = mpn_set_str_other (rp, dp, sn, base, &info);
}
assert (rn <= alloc);
gmp_free (dp);
@@ -3967,14 +4219,9 @@ mpz_out_str (FILE *stream, int base, const mpz_t x)
static int
gmp_detect_endian (void)
{
- static const int i = 1;
+ static const int i = 2;
const unsigned char *p = (const unsigned char *) &i;
- if (*p == 1)
- /* Little endian */
- return -1;
- else
- /* Big endian */
- return 1;
+ return 1 - *p;
}
/* Import and export. Does not support nails. */
@@ -4037,29 +4284,22 @@ mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
}
}
}
- if (bytes > 0)
+ assert (i + (bytes > 0) == rn);
+ if (limb != 0)
rp[i++] = limb;
- assert (i == rn);
+ else
+ i = mpn_normalized_size (rp, i);
- r->_mp_size = mpn_normalized_size (rp, i);
+ r->_mp_size = i;
}
void *
mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
size_t nails, const mpz_t u)
{
- unsigned char *p;
- ptrdiff_t word_step;
- size_t count, k;
+ size_t count;
mp_size_t un;
- /* The current (partial) limb. */
- mp_limb_t 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;
-
if (nails != 0)
gmp_die ("mpz_import: Nails not supported.");
@@ -4067,62 +4307,74 @@ mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
assert (endian >= -1 && endian <= 1);
assert (size > 0 || u->_mp_size == 0);
- un = GMP_ABS (u->_mp_size);
- if (un == 0)
- {
- if (countp)
- *countp = 0;
- return r;
- }
+ un = u->_mp_size;
+ count = 0;
+ if (un != 0)
+ {
+ size_t k;
+ unsigned char *p;
+ ptrdiff_t word_step;
+ /* The current (partial) limb. */
+ mp_limb_t 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;
- /* Count bytes in top limb. */
- for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT)
- ;
+ un = GMP_ABS (un);
- assert (k > 0);
+ /* Count bytes in top limb. */
+ limb = u->_mp_d[un-1];
+ assert (limb != 0);
- count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
+ k = 0;
+ do {
+ k++; limb >>= CHAR_BIT;
+ } while (limb != 0);
- if (!r)
- r = gmp_xalloc (count * size);
+ count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
- if (endian == 0)
- endian = gmp_detect_endian ();
+ if (!r)
+ r = gmp_xalloc (count * size);
- p = (unsigned char *) r;
+ if (endian == 0)
+ endian = gmp_detect_endian ();
- word_step = (order != endian) ? 2 * size : 0;
+ p = (unsigned char *) r;
- /* Process bytes from the least significant end, so point p at the
- least significant word. */
- if (order == 1)
- {
- p += size * (count - 1);
- word_step = - word_step;
- }
+ word_step = (order != endian) ? 2 * size : 0;
- /* And at least significant byte of that word. */
- if (endian == 1)
- p += (size - 1);
+ /* Process bytes from the least significant end, so point p at the
+ least significant word. */
+ if (order == 1)
+ {
+ p += size * (count - 1);
+ word_step = - word_step;
+ }
- 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)
- {
- if (bytes == 0)
- {
- if (i < un)
- limb = u->_mp_d[i++];
- bytes = sizeof (mp_limb_t);
- }
- *p = limb;
- limb >>= CHAR_BIT;
- bytes--;
- }
- }
- assert (i == un);
- assert (k == count);
+ /* And at least significant byte of that word. */
+ if (endian == 1)
+ p += (size - 1);
+
+ 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)
+ {
+ if (bytes == 0)
+ {
+ if (i < un)
+ limb = u->_mp_d[i++];
+ bytes = sizeof (mp_limb_t);
+ }
+ *p = limb;
+ limb >>= CHAR_BIT;
+ bytes--;
+ }
+ }
+ assert (i == un);
+ assert (k == count);
+ }
if (countp)
*countp = count;
diff --git a/mini-gmp.h b/mini-gmp.h
index 8c94ca2..d8f691f 100644
--- a/mini-gmp.h
+++ b/mini-gmp.h
@@ -1,21 +1,32 @@
/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
-Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
+Copyright 2011-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/. */
/* About mini-gmp: This is a minimal implementation of a subset of the
GMP interface. It is intended for inclusion into applications which
@@ -64,8 +75,11 @@ typedef __mpz_struct mpz_t[1];
typedef __mpz_struct *mpz_ptr;
typedef const __mpz_struct *mpz_srcptr;
+extern const int mp_bits_per_limb;
+
void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t);
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);
@@ -84,10 +98,17 @@ mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t);
+int mpn_perfect_square_p (mp_srcptr, mp_size_t);
+mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t);
mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
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);
+
+mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t);
+
mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t);
#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0)
@@ -124,6 +145,10 @@ void mpz_mul_si (mpz_t, const mpz_t, long int);
void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int);
void mpz_mul (mpz_t, const mpz_t, const mpz_t);
void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
+void mpz_addmul_ui (mpz_t, const mpz_t, unsigned long int);
+void mpz_addmul (mpz_t, const mpz_t, const mpz_t);
+void mpz_submul_ui (mpz_t, const mpz_t, unsigned long int);
+void mpz_submul (mpz_t, const mpz_t, const mpz_t);
void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
@@ -147,6 +172,7 @@ void mpz_mod (mpz_t, const mpz_t, const mpz_t);
void mpz_divexact (mpz_t, const mpz_t, const mpz_t);
int mpz_divisible_p (const mpz_t, const mpz_t);
+int mpz_congruent_p (const mpz_t, const mpz_t, const mpz_t);
unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
@@ -176,6 +202,7 @@ int mpz_invert (mpz_t, const mpz_t, const mpz_t);
void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t);
void mpz_sqrt (mpz_t, const mpz_t);
+int mpz_perfect_square_p (const mpz_t);
void mpz_pow_ui (mpz_t, const mpz_t, unsigned long);
void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long);
@@ -188,6 +215,9 @@ int mpz_root (mpz_t, const mpz_t, unsigned long);
void mpz_fac_ui (mpz_t, unsigned long);
void mpz_bin_uiui (mpz_t, unsigned long, unsigned long);
+int
+mpz_probab_prime_p (const mpz_t, int);
+
int mpz_tstbit (const mpz_t, mp_bitcnt_t);
void mpz_setbit (mpz_t, mp_bitcnt_t);
void mpz_clrbit (mpz_t, mp_bitcnt_t);
@@ -211,6 +241,15 @@ double mpz_get_d (const mpz_t);
size_t mpz_size (const mpz_t);
mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t);
+void mpz_realloc2 (mpz_t, mp_bitcnt_t);
+mp_srcptr mpz_limbs_read (mpz_srcptr);
+mp_ptr mpz_limbs_modify (mpz_t, mp_size_t);
+mp_ptr mpz_limbs_write (mpz_t, mp_size_t);
+void mpz_limbs_finish (mpz_t, mp_size_t);
+mpz_srcptr mpz_roinit_n (mpz_t, mp_srcptr, mp_size_t);
+
+#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }}
+
void mpz_set_si (mpz_t, signed long int);
void mpz_set_ui (mpz_t, unsigned long int);
void mpz_set (mpz_t, const mpz_t);
--
1.9.2
View Git Blame Copy Permalink