1 files changed, 68 insertions, 0 deletions
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;
+}