PrimeCandidateSource: add Sophie Germain filtering.

A Sophie Germain prime is a prime p such that 2p+1 is also prime. The
larger prime of the pair 2p+1 is also known as a 'safe prime', and is
the preferred kind of modulus for conventional Diffie-Hellman.

Generating these is harder work than normal prime generation. There's
not really much of a technique except to just keep generating
candidate primes p and then testing 2p+1. But what you _can_ do to
speed things up is to get the prime-candidate generator to help a bit:
it's already enforcing that no small prime divides p, and it's easy to
get it to also enforce that no small prime divides 2p+1. That check
can filter out a lot of bad candidates early, before you waste time on
the more expensive checks, so you have a better chance of success with
each number that gets as far as Miller-Rabin.

Here I add an extra setup function for PrimeCandidateSource which
enables those extra checks. After you call pcs_try_sophie_germain(),
the PCS will only deliver you numbers for which both p and 2p+1 are
free of small factors.

primecandidate.c

@@ -15,7 +15,7 @@ struct avoid {
 
 struct PrimeCandidateSource {
     unsigned bits;
-    bool ready;
+    bool ready, try_sophie_germain;
 
     /* We'll start by making up a random number strictly less than this ... */
     mp_int *limit;
@@ -47,6 +47,7 @@ PrimeCandidateSource *pcs_new_with_firstbits(unsigned bits,
 
     s->bits = bits;
     s->ready = false;
+    s->try_sophie_germain = false;
 
     s->kps = NULL;
     s->nkps = s->kpsize = 0;
@@ -97,6 +98,11 @@ void pcs_free(PrimeCandidateSource *s)
     sfree(s);
 }
 
+void pcs_try_sophie_germain(PrimeCandidateSource *s)
+{
+    s->try_sophie_germain = true;
+}
+
 static void pcs_require_residue_inner(PrimeCandidateSource *s,
                                       mp_int *mod, mp_int *res)
 {
@@ -270,6 +276,29 @@ void pcs_ready(PrimeCandidateSource *s)
     for (size_t i = 0; i < NSMALLPRIMES && smallprimes[i] < limit; i++)
         ADD_AVOID(smallprimes[i], 0);
 
+    if (s->try_sophie_germain) {
+        /*
+         * If we're aiming to generate a Sophie Germain prime (i.e. p
+         * such that 2p+1 is also prime), then we also want to ensure
+         * 2p+1 is not congruent to 0 mod any small prime, because if
+         * it is, we'll waste a lot of time generating a p for which
+         * 2p+1 can't possibly work. So we have to avoid an extra
+         * residue mod each odd q.
+         *
+         * We can simplify:        2p+1 ==  0      (mod q)
+         *                    =>     2p == -1      (mod q)
+         *                    =>      p == -2^{-1} (mod q)
+         *
+         * There's no need to do Euclid's algorithm to compute those
+         * inverses, because for any odd q, the modular inverse of -2
+         * mod q is just (q-1)/2. (Proof: multiplying it by -2 gives
+         * 1-q, which is congruent to 1 mod q.)
+         */
+        for (size_t i = 0; i < NSMALLPRIMES && smallprimes[i] < limit; i++)
+            if (smallprimes[i] != 2)
+                ADD_AVOID(smallprimes[i], (smallprimes[i] - 1) / 2);
+    }
+
     /*
      * Finally, if there's a particular modulus and residue we've been
      * told to avoid, put it on the list.

sshkeygen.h

@@ -51,6 +51,9 @@ void pcs_require_residue_1_mod_prime(PrimeCandidateSource *s, mp_int *mod);
 void pcs_avoid_residue_small(PrimeCandidateSource *s,
                              unsigned mod, unsigned res);
 
+/* Exclude any prime that has no chance of being a Sophie Germain prime. */
+void pcs_try_sophie_germain(PrimeCandidateSource *s);
+
 /* Prepare a PrimeCandidateSource to actually generate numbers. This
  * function does last-minute computation that has to be delayed until
  * all constraints have been input. */

testcrypt.h

@@ -276,6 +276,7 @@ FUNC3(void, pcs_require_residue, val_pcs, val_mpint, val_mpint)
 FUNC2(void, pcs_require_residue_1, val_pcs, val_mpint)
 FUNC2(void, pcs_require_residue_1_mod_prime, val_pcs, val_mpint)
 FUNC3(void, pcs_avoid_residue_small, val_pcs, uint, uint)
+FUNC1(void, pcs_try_sophie_germain, val_pcs)
 FUNC1(void, pcs_ready, val_pcs)
 FUNC4(void, pcs_inspect, val_pcs, out_val_mpint, out_val_mpint, out_val_mpint)
 FUNC1(val_mpint, pcs_generate, val_pcs)