Hi,
the subset of mpn functions in mini-gmp is rather small, but can be
expanded as needed. Would it make sense to add mpn_gcd? Patch below, net
increase about 35 lines.
It's a bit annoying that we need both mpn_make_odd and mpz_make_odd.
One question is on appropriate tests, is it sufficient that the code is
tested only via mpz_gcd? We don't get coverage for the case of even u,
odd v that way.
Regards,
/Niels
diff -r ca451d583385 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 Wed Oct 16 19:38:56 2024 +0200
@@ -2706,6 +2706,75 @@
return u << shift;
}
+/* X must be non-zero. */
+static mp_bitcnt_t
+mpn_make_odd (mp_ptr xp, mp_size_t xn)
+{
+ mp_bitcnt_t bits, shift;
+ mp_size_t limb_cnt;
+
+ /* Count trailing zeros, equivalent to mpn_scan1, because we know that there
is a 1 */
+ bits = mpn_scan1 (xp, 0);
+ limb_cnt = bits / GMP_LIMB_BITS;
+ shift = bits % GMP_LIMB_BITS;
+
+ if (shift > 0)
+ mpn_rshift (xp, xp + limb_cnt, xn - limb_cnt, shift);
+ else
+ mpn_copyi (xp, xp + limb_cnt, xn - limb_cnt);
+
+ mpn_zero (xp + xn - limb_cnt, limb_cnt);
+ return bits;
+}
+
+mp_size_t
+mpn_gcd (mp_ptr rp, mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t n)
+{
+ assert (un >= n);
+ assert (n > 0);
+ assert (!GMP_MPN_OVERLAP_P (up, un, vp, n));
+ assert (vp[n-1]);
+ assert ((up[0] | vp[0]) & 1);
+
+ if (un > n)
+ mpn_div_qr (NULL, up, un, vp, n);
+
+ if (mpn_zero_p (up, n))
+ {
+ mpn_copyi (rp, vp, n);
+ return n;
+ }
+
+ if (!(vp[0] & 1))
+ MP_PTR_SWAP (up, vp);
+
+ while (n > 1)
+ {
+ int c;
+
+ assert (vp[0] & 1);
+
+ mpn_make_odd (up, n);
+ c = mpn_cmp (up, vp, n);
+ if (c == 0)
+ {
+ n = mpn_normalized_size (vp, n);
+ mpn_copyi (rp, vp, n);
+ assert (vp[n-1] > 0);
+ return n;
+ }
+ if (c < 0)
+ MP_PTR_SWAP (up, vp);
+
+ mpn_sub_n (up, up, vp, n);
+
+ while ( (vp[n-1] | up[n-1]) == 0)
+ n--;
+ }
+ 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)
{
@@ -2726,14 +2795,13 @@
static mp_bitcnt_t
mpz_make_odd (mpz_t r)
{
- mp_bitcnt_t shift;
+ mp_bitcnt_t bits;
assert (r->_mp_size > 0);
- /* Count trailing zeros, equivalent to mpn_scan1, because we know that there
is a 1 */
- shift = mpn_scan1 (r->_mp_d, 0);
- mpz_tdiv_q_2exp (r, r, shift);
-
- return shift;
+ bits = mpn_make_odd (r->_mp_d, r->_mp_size);
+ r->_mp_size = mpn_normalized_size (r->_mp_d, r->_mp_size);
+
+ return bits;
}
void
@@ -2765,42 +2833,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 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 Wed Oct 16 19:38:56 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);
--
Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677.
Internet email is subject to wholesale government surveillance.
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
