commit 302edc17336dbd46e4040f77cb36c0f99b736743
Author: Mattias Andrée <maandree_AT_kth.se>
AuthorDate: Thu Mar 3 23:58:26 2016 +0100
Commit: Mattias Andrée <maandree_AT_kth.se>
CommitDate: Thu Mar 3 23:58:26 2016 +0100
Add zptest
Signed-off-by: Mattias Andrée <maandree_AT_kth.se>
diff --git a/man/zptest.3 b/man/zptest.3
index 404b66e..f8193d0 100644
--- a/man/zptest.3
+++ b/man/zptest.3
_AT_@ -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
index ba9bd44..873adbb 100644
--- a/src/internals.h
+++ b/src/internals.h
_AT_@ -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
_AT_@ -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;
+}
Received on Fri Mar 04 2016 - 09:28:30 CET
This archive was generated by hypermail 2.3.0
: Fri Mar 04 2016 - 09:36:17 CET