[Gmp-commit] /var/hg/gmp: mini-gmp: Implement mpn_gcd.
mercurial at gmplib.org
mercurial at gmplib.org
Fri Oct 18 19:17:34 CEST 2024
details: /var/hg/gmp/rev/6df5dd697f5a
changeset: 18474:6df5dd697f5a
user: Niels Möller <nisse at lysator.liu.se>
date: Fri Oct 18 19:16:58 2024 +0200
description:
mini-gmp: Implement mpn_gcd.
diffstat:
mini-gmp/ChangeLog | 7 +++
mini-gmp/mini-gmp.c | 97 ++++++++++++++++++++++++++++++++-----------------
mini-gmp/mini-gmp.h | 1 +
mini-gmp/tests/t-gcd.c | 39 +++++++++++++++++++-
4 files changed, 109 insertions(+), 35 deletions(-)
diffs (197 lines):
diff -r ca451d583385 -r 6df5dd697f5a mini-gmp/ChangeLog
--- a/mini-gmp/ChangeLog Wed May 15 20:51:11 2024 +0200
+++ b/mini-gmp/ChangeLog Fri Oct 18 19:16:58 2024 +0200
@@ -1,3 +1,10 @@
+2024-10-18 Niels Möller <nisse at lysator.liu.se>
+
+ * mini-gmp.c (mpn_gcd): New function.
+ (mpz_gcd): Use mpn_gcd.
+ * mini-gmp.h: Add declaration of mpn_gcd.
+ * tests/t-gcd.c (test_one): Add testing of mpn_gcd.
+
2024-03-25 Marco Bodrato <bodrato at mail.dm.unipi.it>
* tests/t-import.c: Do not print unused value.
diff -r ca451d583385 -r 6df5dd697f5a mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c Wed May 15 20:51:11 2024 +0200
+++ b/mini-gmp/mini-gmp.c Fri Oct 18 19:16:58 2024 +0200
@@ -2706,6 +2706,66 @@
return u << shift;
}
+mp_size_t
+mpn_gcd (mp_ptr rp, mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn)
+{
+ assert (un >= vn);
+ assert (vn > 0);
+ assert (!GMP_MPN_OVERLAP_P (up, un, vp, vn));
+ assert (vp[vn-1] > 0);
+ assert ((up[0] | vp[0]) & 1);
+
+ if (un > vn)
+ mpn_div_qr (NULL, up, un, vp, vn);
+
+ un = mpn_normalized_size (up, vn);
+ if (un == 0)
+ {
+ mpn_copyi (rp, vp, vn);
+ return vn;
+ }
+
+ if (!(vp[0] & 1))
+ MPN_PTR_SWAP (up, un, vp, vn);
+
+ while (un > 1 || vn > 1)
+ {
+ int shift;
+ assert (vp[0] & 1);
+
+ while (up[0] == 0)
+ {
+ up++;
+ un--;
+ }
+ gmp_ctz (shift, up[0]);
+ if (shift > 0)
+ {
+ gmp_assert_nocarry (mpn_rshift(up, up, un, shift));
+ un -= (up[un-1] == 0);
+ }
+
+ if (un < vn)
+ MPN_PTR_SWAP (up, un, vp, vn);
+ else if (un == vn)
+ {
+ int c = mpn_cmp (up, vp, un);
+ if (c == 0)
+ {
+ mpn_copyi (rp, up, un);
+ return un;
+ }
+ else if (c < 0)
+ MP_PTR_SWAP (up, vp);
+ }
+
+ gmp_assert_nocarry (mpn_sub (up, up, un, vp, vn));
+ un = mpn_normalized_size (up, un);
+ }
+ rp[0] = mpn_gcd_11 (up[0], vp[0]);
+ return 1;
+}
+
unsigned long
mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
{
@@ -2765,42 +2825,11 @@
if (tu->_mp_size < tv->_mp_size)
mpz_swap (tu, tv);
- mpz_tdiv_r (tu, tu, tv);
- if (tu->_mp_size == 0)
- {
- mpz_swap (g, tv);
- }
- else
- for (;;)
- {
- int c;
-
- mpz_make_odd (tu);
- c = mpz_cmp (tu, tv);
- if (c == 0)
- {
- mpz_swap (g, tu);
- break;
- }
- if (c < 0)
- mpz_swap (tu, tv);
-
- if (tv->_mp_size == 1)
- {
- mp_limb_t *gp;
-
- mpz_tdiv_r (tu, tu, tv);
- gp = MPZ_REALLOC (g, 1); /* gp = mpz_limbs_modify (g, 1); */
- *gp = mpn_gcd_11 (tu->_mp_d[0], tv->_mp_d[0]);
-
- g->_mp_size = *gp != 0; /* mpz_limbs_finish (g, 1); */
- break;
- }
- mpz_sub (tu, tu, tv);
- }
+ tu->_mp_size = mpn_gcd (tu->_mp_d, tu->_mp_d, tu->_mp_size, tv->_mp_d, tv->_mp_size);
+ mpz_mul_2exp (g, tu, gz);
+
mpz_clear (tu);
mpz_clear (tv);
- mpz_mul_2exp (g, g, gz);
}
void
diff -r ca451d583385 -r 6df5dd697f5a mini-gmp/mini-gmp.h
--- a/mini-gmp/mini-gmp.h Wed May 15 20:51:11 2024 +0200
+++ b/mini-gmp/mini-gmp.h Fri Oct 18 19:16:58 2024 +0200
@@ -105,6 +105,7 @@
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_size_t mpn_gcd (mp_ptr, mp_ptr, mp_size_t, mp_ptr, 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);
diff -r ca451d583385 -r 6df5dd697f5a mini-gmp/tests/t-gcd.c
--- a/mini-gmp/tests/t-gcd.c Wed May 15 20:51:11 2024 +0200
+++ b/mini-gmp/tests/t-gcd.c Fri Oct 18 19:16:58 2024 +0200
@@ -106,7 +106,8 @@
return 1;
}
-static void test_one(const mpz_t a, const mpz_t b)
+static void
+test_one (const mpz_t a, const mpz_t b)
{
mpz_t g, s, t;
@@ -137,6 +138,42 @@
abort ();
}
+ /* Test mpn_gcd, if inputs are valid. */
+ if (mpz_sgn (a) && mpz_sgn (b) && (mpz_odd_p (a) || mpz_odd_p (b)))
+ {
+ mp_size_t an, bn, gn;
+ mp_ptr ap, bp, tp;
+ mpz_t t;
+
+ an = mpz_size (a); ap = a->_mp_d;
+ bn = mpz_size (b); bp = b->_mp_d;
+
+ if (an < bn)
+ {
+ mp_ptr sp = ap;
+ mp_size_t sn = an;
+ ap = bp; an = bn;
+ bp = sp; bn = sn;
+ }
+
+ tp = malloc ((an + bn) * sizeof (mp_limb_t));
+ if (!tp)
+ abort ();
+
+ mpn_copyi (tp, ap, an);
+ mpn_copyi (tp + an, bp, bn);
+ gn = mpn_gcd (tp, tp, an, tp + an, bn);
+ if (mpz_cmp (s, mpz_roinit_n (t, tp, gn)))
+ {
+ fprintf (stderr, "mpn_gcd failed:\n");
+ dump ("a", a);
+ dump ("b", b);
+ dump ("r", t);
+ dump ("ref", s);
+ abort ();
+ }
+ }
+
mpz_clear (g);
mpz_clear (s);
mpz_clear (t);
More information about the gmp-commit
mailing list