diff --git a/Recipe b/Recipe
index df040ab1b383e8931455c3fc3a921172ed98b58e..51b8c21c5101b9e0d22fb146899386ff5a07b678 100644
--- a/Recipe
+++ b/Recipe
@@ -282,7 +282,7 @@ UXSSH = SSH uxnoise uxagentc uxgss uxshare
SFTP = psftpcommon sftp sftpcommon logging cmdline
# Components of the prime-generation system.
-SSHPRIME = sshprime smallprimes primecandidate mpunsafe
+SSHPRIME = sshprime smallprimes primecandidate millerrabin mpunsafe
# Miscellaneous objects appearing in all the utilities, or all the
# network ones, or the Unix or Windows subsets of those in turn.
diff --git a/millerrabin.c b/millerrabin.c
new file mode 100644
index 0000000000000000000000000000000000000000..3358bc51a19aa39c48c98d4320f97ef7a9c351ad
--- /dev/null
+++ b/millerrabin.c
@@ -0,0 +1,214 @@
+/*
+ * millerrabin.c: Miller-Rabin probabilistic primality testing, as
+ * declared in sshkeygen.h.
+ */
+
+#include <assert.h>
+#include "ssh.h"
+#include "sshkeygen.h"
+#include "mpint.h"
+#include "mpunsafe.h"
+
+/*
+ * The Miller-Rabin primality test is an extension to the Fermat
+ * test. The Fermat test just checks that a^(p-1) == 1 mod p; this
+ * is vulnerable to Carmichael numbers. Miller-Rabin considers how
+ * that 1 is derived as well.
+ *
+ * Lemma: if a^2 == 1 (mod p), and p is prime, then either a == 1
+ * or a == -1 (mod p).
+ *
+ * Proof: p divides a^2-1, i.e. p divides (a+1)(a-1). Hence,
+ * since p is prime, either p divides (a+1) or p divides (a-1).
+ * But this is the same as saying that either a is congruent to
+ * -1 mod p or a is congruent to +1 mod p. []
+ *
+ * Comment: This fails when p is not prime. Consider p=mn, so
+ * that mn divides (a+1)(a-1). Now we could have m dividing (a+1)
+ * and n dividing (a-1), without the whole of mn dividing either.
+ * For example, consider a=10 and p=99. 99 = 9 * 11; 9 divides
+ * 10-1 and 11 divides 10+1, so a^2 is congruent to 1 mod p
+ * without a having to be congruent to either 1 or -1.
+ *
+ * So the Miller-Rabin test, as well as considering a^(p-1),
+ * considers a^((p-1)/2), a^((p-1)/4), and so on as far as it can
+ * go. In other words. we write p-1 as q * 2^k, with k as large as
+ * possible (i.e. q must be odd), and we consider the powers
+ *
+ * a^(q*2^0) a^(q*2^1) ... a^(q*2^(k-1)) a^(q*2^k)
+ * i.e. a^((n-1)/2^k) a^((n-1)/2^(k-1)) ... a^((n-1)/2) a^(n-1)
+ *
+ * If p is to be prime, the last of these must be 1. Therefore, by
+ * the above lemma, the one before it must be either 1 or -1. And
+ * _if_ it's 1, then the one before that must be either 1 or -1,
+ * and so on ... In other words, we expect to see a trailing chain
+ * of 1s preceded by a -1. (If we're unlucky, our trailing chain of
+ * 1s will be as long as the list so we'll never get to see what
+ * lies before it. This doesn't count as a test failure because it
+ * hasn't _proved_ that p is not prime.)
+ *
+ * For example, consider a=2 and p=1729. 1729 is a Carmichael
+ * number: although it's not prime, it satisfies a^(p-1) == 1 mod p
+ * for any a coprime to it. So the Fermat test wouldn't have a
+ * problem with it at all, unless we happened to stumble on an a
+ * which had a common factor.
+ *
+ * So. 1729 - 1 equals 27 * 2^6. So we look at
+ *
+ * 2^27 mod 1729 == 645
+ * 2^108 mod 1729 == 1065
+ * 2^216 mod 1729 == 1
+ * 2^432 mod 1729 == 1
+ * 2^864 mod 1729 == 1
+ * 2^1728 mod 1729 == 1
+ *
+ * We do have a trailing string of 1s, so the Fermat test would
+ * have been happy. But this trailing string of 1s is preceded by
+ * 1065; whereas if 1729 were prime, we'd expect to see it preceded
+ * by -1 (i.e. 1728.). Guards! Seize this impostor.
+ *
+ * (If we were unlucky, we might have tried a=16 instead of a=2;
+ * now 16^27 mod 1729 == 1, so we would have seen a long string of
+ * 1s and wouldn't have seen the thing _before_ the 1s. So, just
+ * like the Fermat test, for a given p there may well exist values
+ * of a which fail to show up its compositeness. So we try several,
+ * just like the Fermat test. The difference is that Miller-Rabin
+ * is not _in general_ fooled by Carmichael numbers.)
+ *
+ * Put simply, then, the Miller-Rabin test requires us to:
+ *
+ * 1. write p-1 as q * 2^k, with q odd
+ * 2. compute z = (a^q) mod p.
+ * 3. report success if z == 1 or z == -1.
+ * 4. square z at most k-1 times, and report success if it becomes
+ * -1 at any point.
+ * 5. report failure otherwise.
+ *
+ * (We expect z to become -1 after at most k-1 squarings, because
+ * if it became -1 after k squarings then a^(p-1) would fail to be
+ * 1. And we don't need to investigate what happens after we see a
+ * -1, because we _know_ that -1 squared is 1 modulo anything at
+ * all, so after we've seen a -1 we can be sure of seeing nothing
+ * but 1s.)
+ */
+
+struct MillerRabin {
+ MontyContext *mc;
+
+ size_t k;
+ mp_int *q;
+
+ mp_int *two, *pm1, *m_pm1;
+};
+
+MillerRabin *miller_rabin_new(mp_int *p)
+{
+ MillerRabin *mr = snew(MillerRabin);
+
+ assert(mp_hs_integer(p, 2));
+ assert(mp_get_bit(p, 0) == 1);
+
+ mr->k = 1;
+ while (!mp_get_bit(p, mr->k))
+ mr->k++;
+ mr->q = mp_rshift_safe(p, mr->k);
+
+ mr->two = mp_from_integer(2);
+
+ mr->pm1 = mp_unsafe_copy(p);
+ mp_sub_integer_into(mr->pm1, mr->pm1, 1);
+
+ mr->mc = monty_new(p);
+ mr->m_pm1 = monty_import(mr->mc, mr->pm1);
+
+ return mr;
+}
+
+void miller_rabin_free(MillerRabin *mr)
+{
+ mp_free(mr->q);
+ mp_free(mr->two);
+ mp_free(mr->pm1);
+ mp_free(mr->m_pm1);
+ monty_free(mr->mc);
+ smemclr(mr, sizeof(*mr));
+ sfree(mr);
+}
+
+struct mr_result {
+ bool passed;
+ bool potential_primitive_root;
+};
+
+static struct mr_result miller_rabin_test_inner(MillerRabin *mr, mp_int *w)
+{
+ /*
+ * Compute w^q mod p.
+ */
+ mp_int *wqp = monty_pow(mr->mc, w, mr->q);
+
+ /*
+ * See if this is 1, or if it is -1, or if it becomes -1
+ * when squared at most k-1 times.
+ */
+ struct mr_result result;
+ result.passed = false;
+ result.potential_primitive_root = false;
+
+ if (mp_cmp_eq(wqp, monty_identity(mr->mc))) {
+ result.passed = true;
+ } else {
+ for (size_t i = 0; i < mr->k; i++) {
+ if (mp_cmp_eq(wqp, mr->m_pm1)) {
+ result.passed = true;
+ result.potential_primitive_root = (i == mr->k - 1);
+ break;
+ }
+ if (i == mr->k - 1)
+ break;
+ monty_mul_into(mr->mc, wqp, wqp, wqp);
+ }
+ }
+
+ mp_free(wqp);
+
+ return result;
+}
+
+bool miller_rabin_test_random(MillerRabin *mr)
+{
+ mp_int *mw = mp_random_in_range(mr->two, mr->pm1);
+ struct mr_result result = miller_rabin_test_inner(mr, mw);
+ mp_free(mw);
+ return result.passed;
+}
+
+mp_int *miller_rabin_find_potential_primitive_root(MillerRabin *mr)
+{
+ while (true) {
+ mp_int *mw = mp_unsafe_shrink(mp_random_in_range(mr->two, mr->pm1));
+ struct mr_result result = miller_rabin_test_inner(mr, mw);
+
+ if (result.passed && result.potential_primitive_root) {
+ mp_int *pr = monty_export(mr->mc, mw);
+ mp_free(mw);
+ return pr;
+ }
+
+ mp_free(mw);
+
+ if (!result.passed) {
+ return NULL;
+ }
+ }
+}
+
+unsigned miller_rabin_checks_needed(unsigned bits)
+{
+ /* Table 4.4 from Handbook of Applied Cryptography */
+ return (bits >= 1300 ? 2 : bits >= 850 ? 3 : bits >= 650 ? 4 :
+ bits >= 550 ? 5 : bits >= 450 ? 6 : bits >= 400 ? 7 :
+ bits >= 350 ? 8 : bits >= 300 ? 9 : bits >= 250 ? 12 :
+ bits >= 200 ? 15 : bits >= 150 ? 18 : 27);
+}
+
diff --git a/sshkeygen.h b/sshkeygen.h
index 44952d6ebf38849399cb0dc7929357df9a2ede2f..f82d0dc975253cf5eabddbb8307446d1079df857 100644
--- a/sshkeygen.h
+++ b/sshkeygen.h
@@ -67,6 +67,31 @@ void pcs_inspect(PrimeCandidateSource *pcs, mp_int **limit_out,
/* Query functions for primegen to use */
unsigned pcs_get_bits(PrimeCandidateSource *pcs);
+/* ----------------------------------------------------------------------
+ * A system for doing Miller-Rabin probabilistic primality tests.
+ * These benefit from having set up some context beforehand, if you're
+ * going to do more than one of them on the same candidate prime, so
+ * we declare an object type here to store that context.
+ */
+
+typedef struct MillerRabin MillerRabin;
+
+/* Make and free a Miller-Rabin context. */
+MillerRabin *miller_rabin_new(mp_int *p);
+void miller_rabin_free(MillerRabin *mr);
+
+/* Perform a single Miller-Rabin test, using a random witness value. */
+bool miller_rabin_test_random(MillerRabin *mr);
+
+/* Suggest how many tests are needed to make it sufficiently unlikely
+ * that a composite number will pass them all */
+unsigned miller_rabin_checks_needed(unsigned bits);
+
+/* An extension to the M-R test, which iterates until it either finds
+ * a witness value that is potentially a primitive root, or one
+ * that proves the number to be composite. */
+mp_int *miller_rabin_find_potential_primitive_root(MillerRabin *mr);
+
/* ----------------------------------------------------------------------
* Callback API that allows key generation to report progress to its
* caller.
diff --git a/sshprime.c b/sshprime.c
index 50325295f039ff60bed5ac10aef86c4e90c3aec7..9220d4e6280166fea54c017cb0c39c798b9d3aac 100644
--- a/sshprime.c
+++ b/sshprime.c
@@ -28,98 +28,6 @@
* - go back to square one if any M-R test fails.
*/
-/*
- * The Miller-Rabin primality test is an extension to the Fermat
- * test. The Fermat test just checks that a^(p-1) == 1 mod p; this
- * is vulnerable to Carmichael numbers. Miller-Rabin considers how
- * that 1 is derived as well.
- *
- * Lemma: if a^2 == 1 (mod p), and p is prime, then either a == 1
- * or a == -1 (mod p).
- *
- * Proof: p divides a^2-1, i.e. p divides (a+1)(a-1). Hence,
- * since p is prime, either p divides (a+1) or p divides (a-1).
- * But this is the same as saying that either a is congruent to
- * -1 mod p or a is congruent to +1 mod p. []
- *
- * Comment: This fails when p is not prime. Consider p=mn, so
- * that mn divides (a+1)(a-1). Now we could have m dividing (a+1)
- * and n dividing (a-1), without the whole of mn dividing either.
- * For example, consider a=10 and p=99. 99 = 9 * 11; 9 divides
- * 10-1 and 11 divides 10+1, so a^2 is congruent to 1 mod p
- * without a having to be congruent to either 1 or -1.
- *
- * So the Miller-Rabin test, as well as considering a^(p-1),
- * considers a^((p-1)/2), a^((p-1)/4), and so on as far as it can
- * go. In other words. we write p-1 as q * 2^k, with k as large as
- * possible (i.e. q must be odd), and we consider the powers
- *
- * a^(q*2^0) a^(q*2^1) ... a^(q*2^(k-1)) a^(q*2^k)
- * i.e. a^((n-1)/2^k) a^((n-1)/2^(k-1)) ... a^((n-1)/2) a^(n-1)
- *
- * If p is to be prime, the last of these must be 1. Therefore, by
- * the above lemma, the one before it must be either 1 or -1. And
- * _if_ it's 1, then the one before that must be either 1 or -1,
- * and so on ... In other words, we expect to see a trailing chain
- * of 1s preceded by a -1. (If we're unlucky, our trailing chain of
- * 1s will be as long as the list so we'll never get to see what
- * lies before it. This doesn't count as a test failure because it
- * hasn't _proved_ that p is not prime.)
- *
- * For example, consider a=2 and p=1729. 1729 is a Carmichael
- * number: although it's not prime, it satisfies a^(p-1) == 1 mod p
- * for any a coprime to it. So the Fermat test wouldn't have a
- * problem with it at all, unless we happened to stumble on an a
- * which had a common factor.
- *
- * So. 1729 - 1 equals 27 * 2^6. So we look at
- *
- * 2^27 mod 1729 == 645
- * 2^108 mod 1729 == 1065
- * 2^216 mod 1729 == 1
- * 2^432 mod 1729 == 1
- * 2^864 mod 1729 == 1
- * 2^1728 mod 1729 == 1
- *
- * We do have a trailing string of 1s, so the Fermat test would
- * have been happy. But this trailing string of 1s is preceded by
- * 1065; whereas if 1729 were prime, we'd expect to see it preceded
- * by -1 (i.e. 1728.). Guards! Seize this impostor.
- *
- * (If we were unlucky, we might have tried a=16 instead of a=2;
- * now 16^27 mod 1729 == 1, so we would have seen a long string of
- * 1s and wouldn't have seen the thing _before_ the 1s. So, just
- * like the Fermat test, for a given p there may well exist values
- * of a which fail to show up its compositeness. So we try several,
- * just like the Fermat test. The difference is that Miller-Rabin
- * is not _in general_ fooled by Carmichael numbers.)
- *
- * Put simply, then, the Miller-Rabin test requires us to:
- *
- * 1. write p-1 as q * 2^k, with q odd
- * 2. compute z = (a^q) mod p.
- * 3. report success if z == 1 or z == -1.
- * 4. square z at most k-1 times, and report success if it becomes
- * -1 at any point.
- * 5. report failure otherwise.
- *
- * (We expect z to become -1 after at most k-1 squarings, because
- * if it became -1 after k squarings then a^(p-1) would fail to be
- * 1. And we don't need to investigate what happens after we see a
- * -1, because we _know_ that -1 squared is 1 modulo anything at
- * all, so after we've seen a -1 we can be sure of seeing nothing
- * but 1s.)
- */
-
-static unsigned miller_rabin_checks_needed(unsigned bits)
-{
- /* Table 4.4 from Handbook of Applied Cryptography */
- return (bits >= 1300 ? 2 : bits >= 850 ? 3 : bits >= 650 ? 4 :
- bits >= 550 ? 5 : bits >= 450 ? 6 : bits >= 400 ? 7 :
- bits >= 350 ? 8 : bits >= 300 ? 9 : bits >= 250 ? 12 :
- bits >= 200 ? 15 : bits >= 150 ? 18 : 27);
-}
-
ProgressPhase primegen_add_progress_phase(ProgressReceiver *prog,
unsigned bits)
{
@@ -152,95 +60,35 @@ mp_int *primegen(PrimeCandidateSource *pcs, ProgressReceiver *prog)
{
pcs_ready(pcs);
- STARTOVER:
-
- progress_report_attempt(prog);
+ while (true) {
+ progress_report_attempt(prog);
- mp_int *p = pcs_generate(pcs);
+ mp_int *p = pcs_generate(pcs);
- /*
- * Now apply the Miller-Rabin primality test a few times. First
- * work out how many checks are needed.
- */
- unsigned checks = miller_rabin_checks_needed(pcs_get_bits(pcs));
-
- /*
- * Next, write p-1 as q*2^k.
- */
- size_t k;
- for (k = 0; mp_get_bit(p, k) == !k; k++)
- continue; /* find first 1 bit in p-1 */
- mp_int *q = mp_rshift_safe(p, k);
-
- /*
- * Set up stuff for the Miller-Rabin checks.
- */
- mp_int *two = mp_from_integer(2);
- mp_int *pm1 = mp_copy(p);
- mp_sub_integer_into(pm1, pm1, 1);
- MontyContext *mc = monty_new(p);
- mp_int *m_pm1 = monty_import(mc, pm1);
-
- bool known_bad = false;
-
- /*
- * Now, for each check ...
- */
- for (unsigned check = 0; check < checks && !known_bad; check++) {
- /*
- * Invent a random number between 1 and p-1.
- */
- mp_int *w = mp_random_in_range(two, pm1);
- monty_import_into(mc, w, w);
-
- /*
- * Compute w^q mod p.
- */
- mp_int *wqp = monty_pow(mc, w, q);
- mp_free(w);
-
- /*
- * See if this is 1, or if it is -1, or if it becomes -1
- * when squared at most k-1 times.
- */
- bool passed = false;
-
- if (mp_cmp_eq(wqp, monty_identity(mc)) || mp_cmp_eq(wqp, m_pm1)) {
- passed = true;
- } else {
- for (size_t i = 0; i < k - 1; i++) {
- monty_mul_into(mc, wqp, wqp, wqp);
- if (mp_cmp_eq(wqp, m_pm1)) {
- passed = true;
- break;
- }
+ MillerRabin *mr = miller_rabin_new(p);
+ bool known_bad = false;
+ unsigned nchecks = miller_rabin_checks_needed(mp_get_nbits(p));
+ for (unsigned check = 0; check < nchecks; check++) {
+ if (!miller_rabin_test_random(mr)) {
+ known_bad = true;
+ break;
}
}
+ miller_rabin_free(mr);
+
+ if (!known_bad) {
+ /*
+ * We have a prime!
+ */
+ pcs_free(pcs);
+ return p;
+ }
- if (!passed)
- known_bad = true;
-
- mp_free(wqp);
- }
-
- mp_free(q);
- mp_free(two);
- mp_free(pm1);
- monty_free(mc);
- mp_free(m_pm1);
-
- if (known_bad) {
mp_free(p);
- goto STARTOVER;
}
-
- /*
- * We have a prime!
- */
- pcs_free(pcs);
- return p;
}
+
/* ----------------------------------------------------------------------
* Reusable null implementation of the progress-reporting API.
*/