[hackers] [libzahl] Add zgcd || Mattias Andrée

[git]

commit 5810e189b55eaf51912dace3f338e1c4f68c199c
Author:     Mattias Andrée <maand...@kth.se>
AuthorDate: Wed Mar 2 10:23:14 2016 +0100
Commit:     Mattias Andrée <maand...@kth.se>
CommitDate: Wed Mar 2 10:23:14 2016 +0100

    Add zgcd
    
    Signed-off-by: Mattias Andrée <maand...@kth.se>

diff --git a/src/internals.h b/src/internals.h
index f911680..f7f9769 100644
--- a/src/internals.h
+++ b/src/internals.h
@@ -13,7 +13,9 @@
        X(libzahl_tmp_str_num)\
        X(libzahl_tmp_str_mag)\
        X(libzahl_tmp_str_div)\
-       X(libzahl_tmp_str_rem)
+       X(libzahl_tmp_str_rem)\
+       X(libzahl_tmp_gcd_u)\
+       X(libzahl_tmp_gcd_v)
 
 #define LIST_CONSTS\
        X(libzahl_const_1e19, zsetu, 10000000000000000000ULL) /* The largest 
power of 10 < 2⁶⁴. */\
diff --git a/src/zgcd.c b/src/zgcd.c
new file mode 100644
index 0000000..91216ac
--- /dev/null
+++ b/src/zgcd.c
@@ -0,0 +1,60 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#define u  libzahl_tmp_gcd_u
+#define v  libzahl_tmp_gcd_v
+
+
+void
+zgcd(z_t a, z_t b, z_t c)
+{
+       /*
+        * Binary GCD algorithm.
+        */
+
+       size_t shifts = 0, i = 0;
+       zahl_char_t uv, bit;
+       int neg;
+
+       if (!zcmp(b, c)) {
+               if (a != b)
+                       zset(a, b);
+               return;
+       }
+       if (zzero(b)) {
+               if (a != c)
+                       zset(a, c);
+               return;
+       }
+       if (zzero(c)) {
+               if (a != b)
+                       zset(a, b);
+               return;
+       }
+
+       zabs(u, b);
+       zabs(v, c);
+       neg = zsignum(b) < 0 && zsignum(c) < 0;
+
+       for (;; i++) {
+               uv = (i < u->used ? u->chars[i] : 0)
+                  | (i < v->used ? v->chars[i] : 0);
+               for (bit = 1; bit; bit <<= 1, shifts++)
+                       if (uv & bit)
+                               goto loop_done;
+       }
+loop_done:
+       zrsh(u, u, shifts);
+       zrsh(v, v, shifts);
+
+       zrsh(u, u, zlsb(u));
+       do {
+               zrsh(v, v, zlsb(v));
+               if (zcmpmag(u, v) > 0) /* Both are non-negative. */
+                       zswap(u, v);
+               zsub_unsigned(v, v, u);
+       } while (!zzero(v));
+
+       zlsh(a, u, shifts);
+       SET_SIGNUM(a, neg ? -1 : 1);
+}