2 files changed, 87 insertions, 0 deletions
diff --git a/gnulib b/gnulib
deleted file mode 160000
-Subproject 443bc5ffcf7429e557f4a371b0661abe98ddbc1
diff --git a/gnulib/lib/gcd.c b/gnulib/lib/gcd.c
new file mode 100644
index 0000000..50c60df
--- /dev/null
+++ b/gnulib/lib/gcd.c
@@ -0,0 +1,87 @@
+/* Arithmetic.
+   Copyright (C) 2001-2002, 2006, 2009-2011 Free Software Foundation, Inc.
+   Written by Bruno Haible <bruno@clisp.org>, 2001.
+
+   This program 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.
+
+   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+#include <config.h>
+
+/* This file can also be used to define gcd functions for other unsigned
+   types, such as 'unsigned long long' or 'uintmax_t'.  */
+#ifndef WORD_T
+/* Specification.  */
+# include "gcd.h"
+# define WORD_T unsigned long
+# define GCD gcd
+#endif
+
+#include <stdlib.h>
+
+/* Return the greatest common divisor of a > 0 and b > 0.  */
+WORD_T
+GCD (WORD_T a, WORD_T b)
+{
+  /* Why no division, as in Euclid's algorithm? Because in Euclid's algorithm
+     the division result floor(a/b) or floor(b/a) is very often = 1 or = 2,
+     and nearly always < 8.  A sequence of a few subtractions and tests is
+     faster than a division.  */
+  /* Why not Euclid's algorithm? Because the two integers can be shifted by 1
+     bit in a single instruction, and the algorithm uses fewer variables than
+     Euclid's algorithm.  */
+
+  WORD_T c = a | b;
+  c = c ^ (c - 1);
+  /* c = largest power of 2 that divides a and b.  */
+
+  if (a & c)
+    {
+      if (b & c)
+        goto odd_odd;
+      else
+        goto odd_even;
+    }
+  else
+    {
+      if (b & c)
+        goto even_odd;
+      else
+        abort ();
+    }
+
+  for (;;)
+    {
+    odd_odd: /* a/c and b/c both odd */
+      if (a == b)
+        break;
+      if (a > b)
+        {
+          a = a - b;
+        even_odd: /* a/c even, b/c odd */
+          do
+            a = a >> 1;
+          while ((a & c) == 0);
+        }
+      else
+        {
+          b = b - a;
+        odd_even: /* a/c odd, b/c even */
+          do
+            b = b >> 1;
+          while ((b & c) == 0);
+        }
+    }
+
+  /* a = b */
+  return a;
+}