libzahl

commit 5810e189b55eaf51912dace3f338e1c4f68c199c
parent a626bacf8e45af60727882250f9d3abeb15ff3cd
Author: Mattias Andrée <maandree@kth.se>
Date:   Wed,  2 Mar 2016 10:23:14 +0100

Add zgcd

Signed-off-by: Mattias Andrée <maandree@kth.se>

Diffstat:
M
src/internals.h
 | 
4
+++-
A
src/zgcd.c
 | 
60
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 63 insertions(+), 1 deletion(-)
diff --git 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
@@ -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);
+}