[hackers] [libzahl] Add zptest || Mattias Andrée

[git]

commit 302edc17336dbd46e4040f77cb36c0f99b736743
Author:     Mattias Andrée <maand...@kth.se>
AuthorDate: Thu Mar 3 23:58:26 2016 +0100
Commit:     Mattias Andrée <maand...@kth.se>
CommitDate: Thu Mar 3 23:58:26 2016 +0100

    Add zptest
    
    Signed-off-by: Mattias Andrée <maand...@kth.se>

diff --git a/man/zptest.3 b/man/zptest.3
index 404b66e..f8193d0 100644
--- 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
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;
+}