mini-gmp mpn_gcd?

Niels Möller nisse at lysator.liu.se
Wed Oct 16 20:12:25 CEST 2024


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.



More information about the gmp-devel mailing list