Improve speed of prime generation especially safe primes

First, correct a bug in the sieve code. It would break early if
a value did not match up with the sieve. However in that case,
the sieve values would be out of sync with the value of p, and
would be returning effectively random results. This caused
prime generation to be slower than it should be, both because the
sieve was incorrectly rejecting values that were not multiples of
any small prime and was allowing values that were multiples of
small primes to move on to the Miller-Rabin test.

In the sieve, also sieve so that 2*q+1 is also not a multiple of
the small primes. This speeds up safe prime generation.

GH #1411

3 files changed, 73 insertions, 37 deletions

diff --git a/src/lib/math/mp/mp_core.cpp b/src/lib/math/mp/mp_core.cpp
index 4ed61ac38..ffa5b31a8 100644
--- a/src/lib/math/mp/mp_core.cpp
+++ b/src/lib/math/mp/mp_core.cpp
@@ -436,10 +436,14 @@ word bigint_divop(word n1, word n0, word d)
 */
 word bigint_modop(word n1, word n0, word d)
    {
+#if defined(BOTAN_HAS_MP_DWORD)
+   return ((static_cast<dword>(n1) << MP_WORD_BITS) | n0) % d;
+#else
    word z = bigint_divop(n1, n0, d);
    word dummy = 0;
    z = word_madd2(z, d, &dummy);
    return (n0-z);
+#endif
    }
 
 }
diff --git a/src/lib/math/numbertheory/make_prm.cpp b/src/lib/math/numbertheory/make_prm.cpp
index f06f1978e..419252ed6 100644
--- a/src/lib/math/numbertheory/make_prm.cpp
+++ b/src/lib/math/numbertheory/make_prm.cpp
@@ -16,20 +16,22 @@ namespace Botan {
 */
 BigInt random_prime(RandomNumberGenerator& rng,
                     size_t bits, const BigInt& coprime,
-                    size_t equiv, size_t modulo)
+                    size_t equiv, size_t modulo,
+                    size_t prob)
    {
-   if(coprime <= 0)
+   if(coprime.is_negative())
       {
-      throw Invalid_Argument("random_prime: coprime must be > 0");
+      throw Invalid_Argument("random_prime: coprime must be >= 0");
       }
-   if(modulo % 2 == 1 || modulo == 0)
+   if(modulo == 0)
       {
       throw Invalid_Argument("random_prime: Invalid modulo value");
       }
-   if(equiv >= modulo || equiv % 2 == 0)
-      {
-      throw Invalid_Argument("random_prime: equiv must be < modulo, and odd");
-      }
+
+   equiv %= modulo;
+
+   if(equiv == 0)
+      throw Invalid_Argument("random_prime Invalid value for equiv/modulo");
 
    // Handle small values:
    if(bits <= 1)
@@ -50,6 +52,9 @@ BigInt random_prime(RandomNumberGenerator& rng,
       return ((rng.next_byte() % 2) ? 11 : 13);
       }
 
+   secure_vector<uint16_t> sieve(PRIME_TABLE_SIZE);
+   const size_t MAX_ATTEMPTS = 32*1024;
+
    while(true)
       {
       BigInt p(rng, bits);
@@ -59,21 +64,18 @@ BigInt random_prime(RandomNumberGenerator& rng,
       p.set_bit(bits - 2);
       p.set_bit(0);
 
-      if(p % modulo != equiv)
-         p += (modulo - p % modulo) + equiv;
+      // Force p to be equal to equiv mod modulo
+      p += (modulo - (p % modulo)) + equiv;
 
-      const size_t sieve_size = std::min(bits / 2, PRIME_TABLE_SIZE);
-      secure_vector<uint16_t> sieve(sieve_size);
-
-      for(size_t j = 0; j != sieve.size(); ++j)
-         sieve[j] = static_cast<uint16_t>(p % PRIMES[j]);
+      for(size_t i = 0; i != sieve.size(); ++i)
+         sieve[i] = static_cast<uint16_t>(p % PRIMES[i]);
 
       size_t counter = 0;
       while(true)
          {
          ++counter;
 
-         if(counter >= 4096)
+         if(counter > MAX_ATTEMPTS)
             {
             break; // don't try forever, choose a new starting point
             }
@@ -84,26 +86,39 @@ BigInt random_prime(RandomNumberGenerator& rng,
             break;
 
          bool passes_sieve = true;
-         for(size_t j = 0; j != sieve.size(); ++j)
+         for(size_t i = 0; i != sieve.size(); ++i)
             {
-            sieve[j] = (sieve[j] + modulo) % PRIMES[j];
-            if(sieve[j] == 0)
-               {
+            sieve[i] = (sieve[i] + modulo) % PRIMES[i];
+
+            /*
+            In this case, p is a multiple of PRIMES[i]
+            */
+            if(sieve[i] == 0)
+               passes_sieve = false;
+
+            /*
+            In this case, 2*p+1 will be a multiple of PRIMES[i]
+
+            So if generating a safe prime, we want to avoid this value
+            because 2*p+1 will not be useful. Since the check is cheap to
+            do and doesn't seem to affect the overall distribution of the
+            generated primes overmuch it's used in all cases.
+
+            See "Safe Prime Generation with a Combined Sieve" M. Wiener
+            https://eprint.iacr.org/2003/186.pdf
+            */
+            if(sieve[i] == PRIMES[i] / 2)
                passes_sieve = false;
-               break;
-               }
             }
 
          if(!passes_sieve)
             continue;
 
-         if(gcd(p - 1, coprime) != 1)
+         if(coprime > 0 && gcd(p - 1, coprime) != 1)
             continue;
 
-         if(is_prime(p, rng, 128, true))
-            {
+         if(is_prime(p, rng, prob, true))
             return p;
-            }
          }
       }
    }
@@ -117,12 +132,25 @@ BigInt random_safe_prime(RandomNumberGenerator& rng, size_t bits)
       throw Invalid_Argument("random_safe_prime: Can't make a prime of " +
                              std::to_string(bits) + " bits");
 
-   BigInt p;
-   do
-      p = (random_prime(rng, bits - 1) << 1) + 1;
-   while(!is_prime(p, rng, 128, true));
+   BigInt q, p;
+   for(;;)
+      {
+      /*
+      Generate q == 2 (mod 3)
+
+      Otherwise [q == 1 (mod 3) case], 2*q+1 == 3 (mod 3) and not prime.
+      */
+      q = random_prime(rng, bits - 1, 1, 2, 3, 8);
+      p = (q << 1) + 1;
 
-   return p;
+      if(is_prime(p, rng, 128, true))
+         {
+         // We did only a weak check before, go back and verify q before returning
+         if(is_prime(q, rng, 128, true))
+            return p;
+         }
+
+      }
    }
 
 }
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index 7d37ea1b7..b002ace91 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -164,9 +164,9 @@ size_t BOTAN_PUBLIC_API(2,0) low_zero_bits(const BigInt& x);
 * @return true if all primality tests passed, otherwise false
 */
 bool BOTAN_PUBLIC_API(2,0) is_prime(const BigInt& n,
-                        RandomNumberGenerator& rng,
-                        size_t prob = 56,
-                        bool is_random = false);
+                                    RandomNumberGenerator& rng,
+                                    size_t prob = 56,
+                                    bool is_random = false);
 
 inline bool quick_check_prime(const BigInt& n, RandomNumberGenerator& rng)
    { return is_prime(n, rng, 32); }
@@ -186,11 +186,15 @@ inline bool verify_prime(const BigInt& n, RandomNumberGenerator& rng)
 * @param equiv a non-negative number that the result should be
                equivalent to modulo equiv_mod
 * @param equiv_mod the modulus equiv should be checked against
+* @param prob use test so false positive is bounded by 1/2**prob
 * @return random prime with the specified criteria
 */
 BigInt BOTAN_PUBLIC_API(2,0) random_prime(RandomNumberGenerator& rng,
-                              size_t bits, const BigInt& coprime = 1,
-                              size_t equiv = 1, size_t equiv_mod = 2);
+                                          size_t bits,
+                                          const BigInt& coprime = 0,
+                                          size_t equiv = 1,
+                                          size_t equiv_mod = 2,
+                                          size_t prob = 128);
 
 /**
 * Return a 'safe' prime, of the form p=2*q+1 with q prime
@@ -199,7 +203,7 @@ BigInt BOTAN_PUBLIC_API(2,0) random_prime(RandomNumberGenerator& rng,
 * @return prime randomly chosen from safe primes of length bits
 */
 BigInt BOTAN_PUBLIC_API(2,0) random_safe_prime(RandomNumberGenerator& rng,
-                                   size_t bits);
+                                               size_t bits);
 
 /**
 * Generate DSA parameters using the FIPS 186 kosherizer