New function mpn_gcd_22.

* mpn/generic/gcd.c (gcd_2): Moved to gcd_22.c below.
(mpn_gcd): Adapt for calling gcd_22.
* mpn/generic/gcd_22.c (mpn_gcd_22): New file and function.
* gmp-impl.h (mp_double_limb_t): New (typedef) struct.
* configure.ac (gmp_mpn_functions): Added gcd_22.

* tests/mpn/t-gcd_22.c: New test.
* tests/mpn/Makefile.am (check_PROGRAMS): Add t-gcd_22.
* tests/refmpz.c (refmpz_gcd): New function (plain binary gcd).

9 files changed, 239 insertions, 59 deletions

diff --git a/ChangeLog b/ChangeLog
index 748c83c14..1ac104a10 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+2019-08-16  Niels Möller  <nisse@lysator.liu.se>
+
+	* mpn/generic/gcd.c (gcd_2): Moved to gcd_22.c below.
+	(mpn_gcd): Adapt for calling gcd_22.
+	* mpn/generic/gcd_22.c (mpn_gcd_22): New file and function.
+	* gmp-impl.h (mp_double_limb_t): New (typedef) struct.
+	* configure.ac (gmp_mpn_functions): Added gcd_22.
+
+	* tests/mpn/t-gcd_22.c: New test.
+	* tests/mpn/Makefile.am (check_PROGRAMS): Add t-gcd_22.
+	* tests/refmpz.c (refmpz_gcd): New function (plain binary gcd).
+
 2019-08-13  Torbjörn Granlund  <tg@gmplib.org>
 
 	* mpn/x86_64: Add more gcd_11 variants of of x86_64 gcd_11.asm and
diff --git a/configure.ac b/configure.ac
index 18f1ee5d6..fc102b57e 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2970,7 +2970,7 @@ gmp_mpn_functions="$extra_functions					   \
   rootrem sqrtrem sizeinbase get_str set_str compute_powtab		   \
   scan0 scan1 popcount hamdist cmp zero_p				   \
   perfsqr perfpow strongfibo						   \
-  gcd_11 gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step			   \
+  gcd_11 gcd_22 gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step		   \
   gcdext_lehmer								   \
   div_q tdiv_qr jacbase jacobi_2 jacobi get_d				   \
   matrix22_mul matrix22_mul1_inverse_vector				   \
diff --git a/gmp-impl.h b/gmp-impl.h
index 5c72c5968..293defd8a 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -3968,6 +3968,14 @@ __GMP_DECLSPEC void __gmp_invalid_operation (void) ATTRIBUTE_NORETURN;
 #define PP_FIRST_OMITTED 3
 #endif
 
+typedef struct
+{
+  mp_limb_t d0, d1;
+} mp_double_limb_t;
+
+#define mpn_gcd_22 __MPN (gcd_22)
+__GMP_DECLSPEC mp_double_limb_t mpn_gcd_22 (mp_limb_t, mp_limb_t, mp_limb_t, mp_limb_t);
+
 /* BIT1 means a result value in bit 1 (second least significant bit), with a
    zero bit representing +1 and a one bit representing -1.  Bits other than
    bit 1 are garbage.  These are meant to be kept in "int"s, and casts are
diff --git a/mpn/generic/gcd.c b/mpn/generic/gcd.c
index b36c94892..808d9a826 100644
--- a/mpn/generic/gcd.c
+++ b/mpn/generic/gcd.c
@@ -76,61 +76,6 @@ gcd_hook (void *p, mp_srcptr gp, mp_size_t gn,
   ctx->gn = gn;
 }
 
-#if GMP_NAIL_BITS > 0
-/* Nail supports should be easy, replacing the sub_ddmmss with nails
- * logic. */
-#error Nails not supported.
-#endif
-
-/* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
-   Both U and V must be odd. */
-static inline mp_size_t
-gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
-{
-  mp_limb_t u0, u1, v0, v1;
-  mp_size_t gn;
-
-  u0 = up[0];
-  u1 = up[1];
-  v0 = vp[0];
-  v1 = vp[1];
-
-  ASSERT (u0 & 1);
-  ASSERT (v0 & 1);
-
-  /* Check for u0 != v0 needed to ensure that argument to
-   * count_trailing_zeros is non-zero. */
-  while (u1 != v1 && u0 != v0)
-    {
-      unsigned long int r;
-      if (u1 > v1)
-	{
-	  sub_ddmmss (u1, u0, u1, u0, v1, v0);
-	  count_trailing_zeros (r, u0);
-	  u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
-	  u1 >>= r;
-	}
-      else  /* u1 < v1.  */
-	{
-	  sub_ddmmss (v1, v0, v1, v0, u1, u0);
-	  count_trailing_zeros (r, v0);
-	  v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
-	  v1 >>= r;
-	}
-    }
-
-  gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);
-
-  /* If U == V == GCD, done.  Otherwise, compute GCD (V, |U - V|).  */
-  if (u1 == v1 && u0 == v0)
-    return gn;
-
-  v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
-  gp[0] = mpn_gcd_1 (gp, gn, v0);
-
-  return 1;
-}
-
 mp_size_t
 mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
 {
@@ -301,8 +246,12 @@ mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
       vp[1] >>= r;
     }
 
-  ctx.gn = gcd_2(gp, up, vp);
-
+  {
+    mp_double_limb_t g = mpn_gcd_22 (up[1], up[0], vp[1], vp[0]);
+    gp[0] = g.d0;
+    gp[1] = g.d1;
+    ctx.gn = 1 + (g.d1 > 0);
+  }
 done:
   TMP_FREE;
   return ctx.gn;
diff --git a/mpn/generic/gcd_22.c b/mpn/generic/gcd_22.c
new file mode 100644
index 000000000..9a7366c32
--- /dev/null
+++ b/mpn/generic/gcd_22.c
@@ -0,0 +1,81 @@
+/* mpn_gcd_22 -- double limb greatest common divisor.
+
+Copyright 1994, 1996, 2000, 2001, 2009, 2012, 2019 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 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 General Public License
+for more details.
+
+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/.  */
+
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#if GMP_NAIL_BITS > 0
+#error Nails not supported.
+#endif
+
+mp_double_limb_t
+mpn_gcd_22 (mp_limb_t u1, mp_limb_t u0, mp_limb_t v1, mp_limb_t v0)
+{
+  mp_double_limb_t g;
+  ASSERT (u0 & 1);
+  ASSERT (v0 & 1);
+
+  /* Check for u0 != v0 needed to ensure that argument to
+   * count_trailing_zeros is non-zero. */
+  while (u1 != v1 && u0 != v0)
+    {
+      unsigned long int r;
+      if (u1 > v1)
+	{
+	  sub_ddmmss (u1, u0, u1, u0, v1, v0);
+	  count_trailing_zeros (r, u0);
+	  u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
+	  u1 >>= r;
+	}
+      else  /* u1 < v1.  */
+	{
+	  sub_ddmmss (v1, v0, v1, v0, u1, u0);
+	  count_trailing_zeros (r, v0);
+	  v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
+	  v1 >>= r;
+	}
+    }
+
+  /* If U == V == GCD, done.  Otherwise, compute GCD (V, |U - V|).  */
+  if (u1 == v1 && u0 == v0)
+    {
+      g.d0 = u0; g.d1 = u1;
+      return g;
+    }
+  else {
+    mp_limb_t t[2];
+    t[0] = u0; t[1] = u1;
+
+    v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
+    g.d0 = mpn_gcd_1 (t, 1 + (u1 != 0), v0);
+    g.d1 = 0;
+    return g;
+  }
+}
diff --git a/tests/mpn/Makefile.am b/tests/mpn/Makefile.am
index 9f1cd8868..2fae43900 100644
--- a/tests/mpn/Makefile.am
+++ b/tests/mpn/Makefile.am
@@ -29,7 +29,7 @@ check_PROGRAMS = t-asmtype t-aors_1 t-divrem_1 t-mod_1 t-fat t-get_d	\
   t-toom2-sqr t-toom3-sqr t-toom4-sqr t-toom6-sqr t-toom8-sqr		\
   t-div t-mul t-mullo t-sqrlo t-mulmod_bnm1 t-sqrmod_bnm1 t-mulmid	\
   t-hgcd t-hgcd_appr t-matrix22 t-invert t-bdiv t-fib2m			\
-  t-broot t-brootinv t-minvert t-sizeinbase t-gcd_11
+  t-broot t-brootinv t-minvert t-sizeinbase t-gcd_11 t-gcd_22
 
 EXTRA_DIST = toom-shared.h toom-sqr-shared.h
 
diff --git a/tests/mpn/t-gcd_22.c b/tests/mpn/t-gcd_22.c
new file mode 100644
index 000000000..4761680c9
--- /dev/null
+++ b/tests/mpn/t-gcd_22.c
@@ -0,0 +1,83 @@
+/* Test mpn_gcd_22.
+
+Copyright 2019 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library test suite.
+
+The GNU MP Library test suite is free software; you can redistribute it
+and/or modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation; either version 3 of the License,
+or (at your option) any later version.
+
+The GNU MP Library test suite 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 General
+Public License for more details.
+
+You should have received a copy of the GNU General Public License along with
+the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "gmp-impl.h"
+#include "tests.h"
+
+#ifndef COUNT
+#define COUNT 10000
+#endif
+
+static void
+one_test (mpz_srcptr a, mpz_srcptr b, mpz_srcptr ref)
+{
+  mp_double_limb_t r = mpn_gcd_22 (mpz_getlimbn (a, 1), mpz_getlimbn (a, 0),
+				   mpz_getlimbn (b, 1), mpz_getlimbn (b, 0));
+  if (r.d0 != mpz_getlimbn (ref, 0) || r.d1 != mpz_getlimbn (ref, 1))
+    {
+      gmp_fprintf (stderr,
+		   "gcd_22 (0x%Zx, 0x%Zx) failed, got: g1 = 0x%Mx g0 = %Mx, ref: 0x%Zx\n",
+                   a, b, r.d1, r.d0, ref);
+      abort();
+    }
+}
+
+int
+main (int argc, char **argv)
+{
+  mpz_t a, b, ref;
+  int count = COUNT;
+  int test;
+  gmp_randstate_ptr rands;
+
+  TESTS_REPS (count, argv, argc);
+
+  tests_start ();
+  rands = RANDS;
+
+  mpz_init (a);
+  mpz_init (b);
+  mpz_init (ref);
+  for (test = 0; test < count; test++)
+    {
+      mp_bitcnt_t size = 1 + gmp_urandomm_ui(rands, 2*GMP_NUMB_BITS);
+      if (test & 1)
+	{
+	  mpz_urandomb (a, rands, size);
+	  mpz_urandomb (b, rands, size);
+	}
+      else
+	{
+	  mpz_rrandomb (a, rands, size);
+	  mpz_rrandomb (b, rands, size);
+	}
+
+      mpz_setbit (a, 0);
+      mpz_setbit (b, 0);
+      refmpz_gcd (ref, a, b);
+      one_test (a, b, ref);
+    }
+
+  mpz_clear (a);
+  mpz_clear (b);
+  mpz_clear (ref);
+}
diff --git a/tests/refmpz.c b/tests/refmpz.c
index cb280d534..167799f97 100644
--- a/tests/refmpz.c
+++ b/tests/refmpz.c
@@ -78,6 +78,52 @@ refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
   return ret;
 }
 
+void
+refmpz_gcd (mpz_ptr g, mpz_srcptr a_orig, mpz_srcptr b_orig)
+{
+  mp_bitcnt_t a_twos, b_twos, common_twos;
+  mpz_t a;
+  mpz_t b;
+  mpz_init (a);
+  mpz_init (b);
+  mpz_abs (a, a_orig);
+  mpz_abs (b, b_orig);
+
+  if (mpz_sgn (a) == 0)
+    {
+      mpz_set (g, b);
+      return;
+    }
+  if (mpz_sgn (b) == 0)
+    {
+      mpz_set (g, a);
+      return;
+    }
+  a_twos = mpz_scan1 (a, 0);
+  mpz_tdiv_q_2exp (a, a, a_twos);
+
+  b_twos = mpz_scan1 (b, 0);
+  mpz_tdiv_q_2exp (b, b, b_twos);
+
+  common_twos = MIN(a_twos, b_twos);
+  for (;;)
+    {
+      int c;
+      mp_bitcnt_t twos;
+      c = mpz_cmp (a, b);
+      if (c == 0)
+	break;
+      if (c < 0)
+	mpz_swap (a, b);
+      mpz_sub (a, a, b);
+      twos = mpz_scan1 (a, 0);
+      mpz_tdiv_q_2exp (a, a, twos);
+    }
+  mpz_mul_2exp (g, a, common_twos);
+
+  mpz_clear (a);
+  mpz_clear (b);
+}
 
 /* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
 #define JACOBI_0Z(b)  JACOBI_0LS (PTR(b)[0], SIZ(b))
diff --git a/tests/tests.h b/tests/tests.h
index 0206b2119..a44d39950 100644
--- a/tests/tests.h
+++ b/tests/tests.h
@@ -356,6 +356,7 @@ void refmpq_sub (mpq_ptr, mpq_srcptr, mpq_srcptr);
 
 void refmpz_combit (mpz_ptr, unsigned long);
 unsigned long refmpz_hamdist (mpz_srcptr, mpz_srcptr);
+void refmpz_gcd (mpz_ptr, mpz_srcptr, mpz_srcptr);
 int refmpz_kronecker (mpz_srcptr, mpz_srcptr);
 int refmpz_jacobi (mpz_srcptr, mpz_srcptr);
 int refmpz_legendre (mpz_srcptr, mpz_srcptr);