#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
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;
/*
* Generate a k-bit random number with top and bottom bits set.
*/
- p = newbn((bits+15)/16);
+ p = bn_power_2(bits-1);
for (i = 0; i < bits; i++) {
if (i == 0 || i == bits-1)
v = 1;
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++)
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;
*/
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 ...
* Invent a random number between 1 and p-1 inclusive.
*/
while (1) {
- w = newbn((bits+15)/16);
+ w = bn_power_2(bits-1);
for (i = 0; i < bits; i++) {
if (bitsleft <= 0)
bitsleft = 8; byte = random_byte();
bitsleft--;
bignum_set_bit(w, i, v);
}
+ bn_restore_invariant(w);
if (bignum_cmp(w, p) >= 0 || bignum_cmp(w, Zero) == 0) {
freebn(w);
continue;
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;
}
* compositeness of p.
*/
freebn(p);
+ freebn(pm1);
freebn(q);
goto STARTOVER;
}
* We have a prime!
*/
freebn(q);
+ freebn(pm1);
return p;
}
projects / u / mdw / putty / blobdiff