From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001
From: Duncan Wilkie <antigravityd@gmail.com>
Date: Sat, 18 Nov 2023 06:11:09 -0600
Subject: Initial commit.

---
 gmp-6.3.0/mpz/gcd.c | 164 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 164 insertions(+)
 create mode 100644 gmp-6.3.0/mpz/gcd.c

(limited to 'gmp-6.3.0/mpz/gcd.c')

diff --git a/gmp-6.3.0/mpz/gcd.c b/gmp-6.3.0/mpz/gcd.c
new file mode 100644
index 0000000..9557155
--- /dev/null
+++ b/gmp-6.3.0/mpz/gcd.c
@@ -0,0 +1,164 @@
+/* mpz/gcd.c:   Calculate the greatest common divisor of two integers.
+
+Copyright 1991, 1993, 1994, 1996, 2000-2002, 2005, 2010, 2015, 2016
+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"
+
+
+void
+mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
+{
+  unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
+  mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
+  mp_ptr tp;
+  mp_ptr up;
+  mp_size_t usize;
+  mp_ptr vp;
+  mp_size_t vsize;
+  mp_size_t gsize;
+  TMP_DECL;
+
+  up = PTR(u);
+  usize = ABSIZ (u);
+  vp = PTR(v);
+  vsize = ABSIZ (v);
+  /* GCD(0, V) == V.  */
+  if (usize == 0)
+    {
+      SIZ (g) = vsize;
+      if (g == v)
+	return;
+      tp = MPZ_NEWALLOC (g, vsize);
+      MPN_COPY (tp, vp, vsize);
+      return;
+    }
+
+  /* GCD(U, 0) == U.  */
+  if (vsize == 0)
+    {
+      SIZ (g) = usize;
+      if (g == u)
+	return;
+      tp = MPZ_NEWALLOC (g, usize);
+      MPN_COPY (tp, up, usize);
+      return;
+    }
+
+  if (usize == 1)
+    {
+      SIZ (g) = 1;
+      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (vp, vsize, up[0]);
+      return;
+    }
+
+  if (vsize == 1)
+    {
+      SIZ(g) = 1;
+      MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (up, usize, vp[0]);
+      return;
+    }
+
+  TMP_MARK;
+
+  /*  Eliminate low zero bits from U and V and move to temporary storage.  */
+  tp = up;
+  while (*tp == 0)
+    tp++;
+  u_zero_limbs = tp - up;
+  usize -= u_zero_limbs;
+  count_trailing_zeros (u_zero_bits, *tp);
+  up = TMP_ALLOC_LIMBS (usize);
+  if (u_zero_bits != 0)
+    {
+      mpn_rshift (up, tp, usize, u_zero_bits);
+      usize -= up[usize - 1] == 0;
+    }
+  else
+    MPN_COPY (up, tp, usize);
+
+  tp = vp;
+  while (*tp == 0)
+    tp++;
+  v_zero_limbs = tp - vp;
+  vsize -= v_zero_limbs;
+  count_trailing_zeros (v_zero_bits, *tp);
+  vp = TMP_ALLOC_LIMBS (vsize);
+  if (v_zero_bits != 0)
+    {
+      mpn_rshift (vp, tp, vsize, v_zero_bits);
+      vsize -= vp[vsize - 1] == 0;
+    }
+  else
+    MPN_COPY (vp, tp, vsize);
+
+  if (u_zero_limbs > v_zero_limbs)
+    {
+      g_zero_limbs = v_zero_limbs;
+      g_zero_bits = v_zero_bits;
+    }
+  else
+    {
+      g_zero_limbs = u_zero_limbs;
+      if (u_zero_limbs < v_zero_limbs)
+	g_zero_bits = u_zero_bits;
+      else  /*  Equal.  */
+	g_zero_bits = MIN (u_zero_bits, v_zero_bits);
+    }
+
+  /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
+  vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
+    ? mpn_gcd (vp, vp, vsize, up, usize)
+    : mpn_gcd (vp, up, usize, vp, vsize);
+
+  /*  Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits).  */
+  gsize = vsize + g_zero_limbs;
+  if (g_zero_bits != 0)
+    {
+      mp_limb_t cy_limb;
+      gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
+      tp = MPZ_NEWALLOC (g, gsize);
+      MPN_ZERO (tp, g_zero_limbs);
+
+      tp = tp + g_zero_limbs;
+      cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
+      if (cy_limb != 0)
+	tp[vsize] = cy_limb;
+    }
+  else
+    {
+      tp = MPZ_NEWALLOC (g, gsize);
+      MPN_ZERO (tp, g_zero_limbs);
+      MPN_COPY (tp + g_zero_limbs, vp, vsize);
+    }
+
+  SIZ (g) = gsize;
+  TMP_FREE;
+}
-- 
cgit v1.2.3