libzahl

commit 302edc17336dbd46e4040f77cb36c0f99b736743
parent 3e4e851bfa30869e48cee9e15edc6723bee684f5
Author: Mattias Andrée <maandree@kth.se>
Date:   Thu,  3 Mar 2016 23:58:26 +0100

Add zptest

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

Diffstat:
M
man/zptest.3
 | 
5
+++++
M
src/internals.h
 | 
11
+++++++++--
A
src/zptest.c
 | 
68
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

3 files changed, 82 insertions(+), 2 deletions(-)
diff --git a/man/zptest.3 b/man/zptest.3
@@ -35,6 +35,11 @@ This factor can be either composite, prime, or 1.
 The risk that a composite number is determined to be
 a probably prime is
 .IR (1-4↑-tries) .
+.P
+It is safe to call
+.B zptest
+with non-unique parameters, and with
+.IR "(witness==0)" .
 .SH RETURN VALUE
 This function will either return:
 .TP
diff --git 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
@@ -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;
+}