X-Git-Url: https://git.distorted.org.uk/u/mdw/putty/blobdiff_plain/59600f67b6dd739bb31a2cb23e887adef13c3ac8..ded53fce5af6818acc2da88653b56c5c94377b8e:/sshprime.c

diff --git a/sshprime.c b/sshprime.c
index 4a2a660e..61488a07 100644
--- a/sshprime.c
+++ b/sshprime.c
@@ -5,6 +5,107 @@
 #include "ssh.h"
 
 /*
+ * This prime generation algorithm is pretty much cribbed from
+ * OpenSSL. The algorithm is:
+ * 
+ *  - invent a B-bit random number and ensure the top and bottom
+ *    bits are set (so it's definitely B-bit, and it's definitely
+ *    odd)
+ * 
+ *  - see if it's coprime to all primes below 2^16; increment it by
+ *    two until it is (this shouldn't take long in general)
+ * 
+ *  - perform the Miller-Rabin primality test enough times to
+ *    ensure the probability of it being composite is 2^-80 or
+ *    less
+ * 
+ *  - 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.)
+ */
+
+/*
  * The first few odd primes.
  *
  * import sys
@@ -561,7 +662,7 @@ Bignum primegen(int bits, int modulus, int residue,
                 int phase, progfn_t pfn, void *pfnparam) {
     int i, k, v, byte, bitsleft, check, checks;
     unsigned long delta, moduli[NPRIMES+1], residues[NPRIMES+1];
-    Bignum p, q, wqp, wqp2;
+    Bignum p, pm1, q, wqp, wqp2;
     int progress = 0;
 
     byte = 0; bitsleft = 0;
@@ -596,7 +697,8 @@ Bignum primegen(int bits, int modulus, int residue,
         residues[i] = bignum_mod_short(p, primes[i]);
     }
     moduli[NPRIMES] = modulus;
-    residues[NPRIMES] = bignum_mod_short(p, modulus) + modulus - residue;
+    residues[NPRIMES] = (bignum_mod_short(p, (unsigned short)modulus)
+                         + modulus - residue);
     delta = 0;
     while (1) {
         for (i = 0; i < (sizeof(moduli) / sizeof(*moduli)); i++)
@@ -617,8 +719,8 @@ Bignum primegen(int bits, int modulus, int residue,
     freebn(q);
 
     /*
-     * Now apply the Fermat primality test a few times. First work
-     * out how many checks are needed.
+     * Now apply the Miller-Rabin primality test a few times. First
+     * work out how many checks are needed.
      */
     checks = 27;
     if (bits >= 150) checks = 18;
@@ -638,6 +740,9 @@ Bignum primegen(int bits, int modulus, int residue,
      */
     for (k = 0; bignum_bit(p, k) == !k; k++);   /* find first 1 bit in p-1 */
     q = bignum_rshift(p, k);
+    /* And store p-1 itself, which we'll need. */
+    pm1 = copybn(p);
+    decbn(pm1);
 
     /*
      * Now, for each check ...
@@ -674,21 +779,21 @@ Bignum primegen(int bits, int modulus, int residue,
         freebn(w);
 
         /*
-         * See if this is 1, or if it becomes 1 when squared at
-         * most k times.
+         * See if this is 1, or if it is -1, or if it becomes -1
+         * when squared at most k-1 times.
          */
-        if (bignum_cmp(wqp, One) == 0) {
+        if (bignum_cmp(wqp, One) == 0 || bignum_cmp(wqp, pm1) == 0) {
             freebn(wqp);
             continue;
         }
-        for (i = 0; i < k; i++) {
+        for (i = 0; i < k-1; i++) {
             wqp2 = modmul(wqp, wqp, p);
             freebn(wqp);
             wqp = wqp2;
-            if (bignum_cmp(wqp, One) == 0)
+            if (bignum_cmp(wqp, pm1) == 0)
                 break;
         }
-        if (i < k) {
+        if (i < k-1) {
             freebn(wqp);
             continue;
         }
@@ -698,6 +803,7 @@ Bignum primegen(int bits, int modulus, int residue,
          * compositeness of p.
          */
         freebn(p);
+        freebn(pm1);
         freebn(q);
         goto STARTOVER;
     }
@@ -706,5 +812,6 @@ Bignum primegen(int bits, int modulus, int residue,
      * We have a prime!
      */
     freebn(q);
+    freebn(pm1);
     return p;
 }