Add zptest

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

2 files changed, 77 insertions, 2 deletions

diff --git a/src/internals.h b/src/internals.h
index ba9bd44..873adbb 100644
--- a/src/internals.h
+++ b/src/internals.h
@@ -29,12 +29,19 @@
 	X(libzahl_tmp_modsqr)\
 	X(libzahl_tmp_divmod_a)\
 	X(libzahl_tmp_divmod_b)\
-	X(libzahl_tmp_divmod_d)
+	X(libzahl_tmp_divmod_d)\
+	X(libzahl_tmp_ptest_x)\
+	X(libzahl_tmp_ptest_a)\
+	X(libzahl_tmp_ptest_d)\
+	X(libzahl_tmp_ptest_n1)\
+	X(libzahl_tmp_ptest_n4)
 
 #define LIST_CONSTS\
 	X(libzahl_const_1e19, zsetu, 10000000000000000000ULL) /* The largest power of 10 < 2⁶⁴. */\
 	X(libzahl_const_1e9,  zsetu, 1000000000ULL)           /* The largest power of 10 < 2³². */\
-	X(libzahl_const_1,    zsetu, 1)
+	X(libzahl_const_1,    zsetu, 1)\
+	X(libzahl_const_2,    zsetu, 2)\
+	X(libzahl_const_4,    zsetu, 4)
 
 #define X(x)  extern z_t x;
 LIST_TEMPS
diff --git a/src/zptest.c b/src/zptest.c
new file mode 100644
index 0000000..db40527
--- /dev/null
+++ b/src/zptest.c
@@ -0,0 +1,68 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#define x   libzahl_tmp_ptest_x
+#define a   libzahl_tmp_ptest_a
+#define d   libzahl_tmp_ptest_d
+#define n1  libzahl_tmp_ptest_n1
+#define n2  libzahl_tmp_ptest_n4
+
+
+enum zprimality
+zptest(z_t witness, z_t n, int t)
+{
+	/*
+	 * Miller–Rabin primarlity test.
+	 */
+
+	size_t i, r;
+
+	if (zcmpu(n, 3) <= 0) {
+		if (zcmpu(n, 1) <= 0) {
+			if (witness)
+				SET(witness, n);
+			return NONPRIME;
+		} else {
+			return PRIME;
+		}
+	}
+	if (zeven(n)) {
+		if (witness)
+			SET(witness, n);
+		return NONPRIME;
+	}
+
+	zsub_unsigned(n1, n, libzahl_const_1);
+	zsub_unsigned(n4, n, libzahl_const_4);
+
+	r = zlsb(n1);
+	zrsh(d, n1, r);
+
+	while (t--) {
+		zrand(a, n4, FAST_RANDOM, UNIFORM);
+		zadd_unsigned(a, a, libzahl_const_2);
+		zmodpow(x, a, d, tn);
+
+		if (!zcmp(x, libzahl_const_1) || !zcmp(x, n1))
+			continue;
+
+		for (i = 1; i < r; i++) {
+			zsqr(x, x);
+			zmod(x, x, tn);
+			if (!zcmp(x, libzahl_const_1)) {
+				if (witness)
+					zswap(witness, a);
+				return NONPRIME;
+			}
+			if (!zcmp(x, n1))
+				break;
+		}
+		if (i == r) {
+			if (witness)
+				zswap(witness, a);
+			return NONPRIME;
+		}
+	}
+
+	return PROBABLY_PRIME;
+}