Merge GH #1636 Add Lucas primality test

16 files changed, 538 insertions, 147 deletions

diff --git a/src/cli/speed.cpp b/src/cli/speed.cpp
index 0efb34ba1..cf35826e5 100644
--- a/src/cli/speed.cpp
+++ b/src/cli/speed.cpp
@@ -92,6 +92,7 @@
    #include <botan/pow_mod.h>
    #include <botan/reducer.h>
    #include <botan/curve_nistp.h>
+   #include <botan/internal/primality.h>
 #endif
 
 #if defined(BOTAN_HAS_ECC_GROUP)
@@ -944,6 +945,10 @@ class Speed final : public Command
 #endif
 
 #if defined(BOTAN_HAS_NUMBERTHEORY)
+            else if(algo == "primality_test")
+               {
+               bench_primality_tests(msec);
+               }
             else if(algo == "random_prime")
                {
                bench_random_prime(msec);
@@ -1700,6 +1705,38 @@ class Speed final : public Command
             }
          }
 
+      void bench_primality_tests(const std::chrono::milliseconds runtime)
+         {
+         for(size_t bits : { 256, 512, 1024 })
+            {
+            std::unique_ptr<Timer> mr_timer = make_timer("Miller-Rabin-" + std::to_string(bits));
+            std::unique_ptr<Timer> bpsw_timer = make_timer("Bailie-PSW-" + std::to_string(bits));
+            std::unique_ptr<Timer> lucas_timer = make_timer("Lucas-" + std::to_string(bits));
+
+            Botan::BigInt n = Botan::random_prime(rng(), bits);
+
+            while(lucas_timer->under(runtime))
+               {
+               Botan::Modular_Reducer mod_n(n);
+
+               mr_timer->run([&]() {
+                  return Botan::is_miller_rabin_probable_prime(n, mod_n, rng(), 2); });
+
+               bpsw_timer->run([&]() {
+                  return Botan::is_bailie_psw_probable_prime(n, mod_n); });
+
+               lucas_timer->run([&]() {
+                  return Botan::is_lucas_probable_prime(n, mod_n); });
+
+               n += 2;
+               }
+
+            record_result(mr_timer);
+            record_result(bpsw_timer);
+            record_result(lucas_timer);
+            }
+         }
+
       void bench_random_prime(const std::chrono::milliseconds runtime)
          {
          const size_t coprime = 65537; // simulates RSA key gen
diff --git a/src/lib/math/bigint/bigint.cpp b/src/lib/math/bigint/bigint.cpp
index 495907d1a..5283c893c 100644
--- a/src/lib/math/bigint/bigint.cpp
+++ b/src/lib/math/bigint/bigint.cpp
@@ -341,6 +341,21 @@ void BigInt::binary_decode(const uint8_t buf[], size_t length)
       m_reg[length / WORD_BYTES] = (m_reg[length / WORD_BYTES] << 8) | buf[i];
    }
 
+void BigInt::ct_cond_assign(bool predicate, BigInt& other)
+   {
+   const size_t t_words = size();
+   const size_t o_words = other.size();
+
+   const size_t r_words = std::max(t_words, o_words);
+
+   const word mask = CT::expand_mask<word>(predicate);
+
+   for(size_t i = 0; i != r_words; ++i)
+      {
+      this->set_word_at(i, CT::select<word>(mask, other.word_at(i), this->word_at(i)));
+      }
+   }
+
 #if defined(BOTAN_HAS_VALGRIND)
 void BigInt::const_time_poison() const
    {
diff --git a/src/lib/math/bigint/bigint.h b/src/lib/math/bigint/bigint.h
index 4a07723b7..8e09e4283 100644
--- a/src/lib/math/bigint/bigint.h
+++ b/src/lib/math/bigint/bigint.h
@@ -616,6 +616,12 @@ class BOTAN_PUBLIC_API(2,0) BigInt final
      */
      void encode_words(word out[], size_t size) const;
 
+     /**
+     * If predicate is true assign other to *this
+     * Uses a masked operation to avoid side channels
+     */
+     void ct_cond_assign(bool predicate, BigInt& other);
+
 #if defined(BOTAN_HAS_VALGRIND)
      void const_time_poison() const;
      void const_time_unpoison() const;
diff --git a/src/lib/math/numbertheory/info.txt b/src/lib/math/numbertheory/info.txt
index 8da52f9f4..84bd24434 100644
--- a/src/lib/math/numbertheory/info.txt
+++ b/src/lib/math/numbertheory/info.txt
@@ -13,6 +13,7 @@ monty.h
 </header:public>
 
 <header:internal>
+primality.h
 def_powm.h
 monty_exp.h
 </header:internal>
diff --git a/src/lib/math/numbertheory/jacobi.cpp b/src/lib/math/numbertheory/jacobi.cpp
index d3e8d7557..284fc2b20 100644
--- a/src/lib/math/numbertheory/jacobi.cpp
+++ b/src/lib/math/numbertheory/jacobi.cpp
@@ -14,12 +14,11 @@ namespace Botan {
 */
 int32_t jacobi(const BigInt& a, const BigInt& n)
    {
-   if(a.is_negative())
-      throw Invalid_Argument("jacobi: first argument must be non-negative");
    if(n.is_even() || n < 2)
       throw Invalid_Argument("jacobi: second argument must be odd and > 1");
 
-   BigInt x = a, y = n;
+   BigInt x = a % n;
+   BigInt y = n;
    int32_t J = 1;
 
    while(y > 1)
diff --git a/src/lib/math/numbertheory/monty.cpp b/src/lib/math/numbertheory/monty.cpp
index b91560fd5..61a10eae5 100644
--- a/src/lib/math/numbertheory/monty.cpp
+++ b/src/lib/math/numbertheory/monty.cpp
@@ -13,7 +13,7 @@ namespace Botan {
 Montgomery_Params::Montgomery_Params(const BigInt& p,
                                      const Modular_Reducer& mod_p)
    {
-   if(p.is_negative() || p.is_even())
+   if(p.is_even() || p < 3)
       throw Invalid_Argument("Montgomery_Params invalid modulus");
 
    m_p = p;
diff --git a/src/lib/math/numbertheory/numthry.cpp b/src/lib/math/numbertheory/numthry.cpp
index a312ba3a1..5bf649407 100644
--- a/src/lib/math/numbertheory/numthry.cpp
+++ b/src/lib/math/numbertheory/numthry.cpp
@@ -14,6 +14,7 @@
 #include <botan/internal/mp_core.h>
 #include <botan/internal/ct_utils.h>
 #include <botan/internal/monty_exp.h>
+#include <botan/internal/primality.h>
 #include <algorithm>
 
 namespace Botan {
@@ -434,78 +435,43 @@ BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
       }
    }
 
-namespace {
 
-bool mr_witness(BigInt&& y,
-                const Modular_Reducer& reducer_n,
-                const BigInt& n_minus_1, size_t s)
+BigInt is_perfect_square(const BigInt& C)
    {
-   if(y == 1 || y == n_minus_1)
-      return false;
-
-   for(size_t i = 1; i != s; ++i)
-      {
-      y = reducer_n.square(y);
-
-      if(y == 1) // found a non-trivial square root
-         return true;
-
-      /*
-      -1 is the trivial square root of unity, so ``a`` is not a
-      witness for this number - give up
-      */
-      if(y == n_minus_1)
-         return false;
-      }
-
-   return true; // is a witness
-   }
+   if(C < 1)
+      throw Invalid_Argument("is_perfect_square requires C >= 1");
+   if(C == 1)
+      return 1;
 
-size_t mr_test_iterations(size_t n_bits, size_t prob, bool random)
-   {
-   const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
+   const size_t n = C.bits();
+   const size_t m = (n + 1) / 2;
+   const BigInt B = C + BigInt::power_of_2(m);
 
-   /*
-   * If the candidate prime was maliciously constructed, we can't rely
-   * on arguments based on p being random.
-   */
-   if(random == false)
-      return base;
+   BigInt X = BigInt::power_of_2(m) - 1;
+   BigInt X2 = (X*X);
 
-   /*
-   * For randomly chosen numbers we can use the estimates from
-   * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
-   *
-   * These values are derived from the inequality for p(k,t) given on
-   * the second page.
-   */
-   if(prob <= 128)
+   for(;;)
       {
-      if(n_bits >= 1536)
-         return 4; // < 2^-133
-      if(n_bits >= 1024)
-         return 6; // < 2^-133
-      if(n_bits >= 512)
-         return 12; // < 2^-129
-      if(n_bits >= 256)
-         return 29; // < 2^-128
+      X = (X2 + C) / (2*X);
+      X2 = (X*X);
+
+      if(X2 < B)
+         break;
       }
 
-   /*
-   If the user desires a smaller error probability than we have
-   precomputed error estimates for, just fall back to using the worst
-   case error rate.
-   */
-   return base;
+   if(X2 == C)
+      return X;
+   else
+      return 0;
    }
 
-}
-
 /*
 * Test for primality using Miller-Rabin
 */
-bool is_prime(const BigInt& n, RandomNumberGenerator& rng,
-              size_t prob, bool is_random)
+bool is_prime(const BigInt& n,
+              RandomNumberGenerator& rng,
+              size_t prob,
+              bool is_random)
    {
    if(n == 2)
       return true;
@@ -520,47 +486,21 @@ bool is_prime(const BigInt& n, RandomNumberGenerator& rng,
       return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
       }
 
-   const size_t test_iterations =
-      mr_test_iterations(n.bits(), prob, is_random && rng.is_seeded());
-
-   const BigInt n_minus_1 = n - 1;
-   const size_t s = low_zero_bits(n_minus_1);
-   const BigInt nm1_s = n_minus_1 >> s;
-   const size_t n_bits = n.bits();
+   const size_t t = miller_rabin_test_iterations(n.bits(), prob, is_random);
 
-   const Modular_Reducer mod_n(n);
-   auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
+   Modular_Reducer mod_n(n);
 
-   const size_t powm_window = 4;
-
-   for(size_t i = 0; i != test_iterations; ++i)
+   if(rng.is_seeded())
       {
-      BigInt a;
-
-      if(rng.is_seeded())
-         {
-         a = BigInt::random_integer(rng, 2, n_minus_1);
-         }
-      else
-         {
-         /*
-         * If passed a null RNG just use 2,3,5, ... as bases
-         *
-         * This is not ideal but in certain circumstances we need to
-         * test for primality but have no RNG available.
-         */
-         a = PRIMES[i];
-         }
-
-      auto powm_a_n = monty_precompute(monty_n, a, powm_window);
-
-      BigInt y = monty_execute(*powm_a_n, nm1_s, n_bits);
-
-      if(mr_witness(std::move(y), mod_n, n_minus_1, s))
+      if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false)
          return false;
-      }
 
-   return true;
+      return is_lucas_probable_prime(n, mod_n);
+      }
+   else
+      {
+      return is_bailie_psw_probable_prime(n, mod_n);
+      }
    }
 
 }
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index 7097979bd..7a978cd3f 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -1,6 +1,6 @@
 /*
 * Number Theory Functions
-* (C) 1999-2007 Jack Lloyd
+* (C) 1999-2007,2018 Jack Lloyd
 *
 * Botan is released under the Simplified BSD License (see license.txt)
 */
@@ -22,8 +22,8 @@ class RandomNumberGenerator;
 * @return (a*b)+c
 */
 BigInt BOTAN_PUBLIC_API(2,0) mul_add(const BigInt& a,
-                         const BigInt& b,
-                         const BigInt& c);
+                                     const BigInt& b,
+                                     const BigInt& c);
 
 /**
 * Fused subtract-multiply
@@ -33,8 +33,8 @@ BigInt BOTAN_PUBLIC_API(2,0) mul_add(const BigInt& a,
 * @return (a-b)*c
 */
 BigInt BOTAN_PUBLIC_API(2,0) sub_mul(const BigInt& a,
-                         const BigInt& b,
-                         const BigInt& c);
+                                     const BigInt& b,
+                                     const BigInt& c);
 
 /**
 * Fused multiply-subtract
@@ -44,8 +44,8 @@ BigInt BOTAN_PUBLIC_API(2,0) sub_mul(const BigInt& a,
 * @return (a*b)-c
 */
 BigInt BOTAN_PUBLIC_API(2,0) mul_sub(const BigInt& a,
-                         const BigInt& b,
-                         const BigInt& c);
+                                     const BigInt& b,
+                                     const BigInt& c);
 
 /**
 * Return the absolute value
@@ -108,8 +108,8 @@ BigInt BOTAN_PUBLIC_API(2,0) ct_inverse_mod_odd_modulus(const BigInt& n, const B
 * Not const time
 */
 size_t BOTAN_PUBLIC_API(2,0) almost_montgomery_inverse(BigInt& result,
-                                           const BigInt& a,
-                                           const BigInt& b);
+                                                       const BigInt& a,
+                                                       const BigInt& b);
 
 /**
 * Call almost_montgomery_inverse and correct the result to a^-1 mod b
@@ -126,8 +126,7 @@ BigInt BOTAN_PUBLIC_API(2,0) normalized_montgomery_inverse(const BigInt& a, cons
 * @param n is an odd integer > 1
 * @return (n / m)
 */
-int32_t BOTAN_PUBLIC_API(2,0) jacobi(const BigInt& a,
-                        const BigInt& n);
+int32_t BOTAN_PUBLIC_API(2,0) jacobi(const BigInt& a, const BigInt& n);
 
 /**
 * Modular exponentation
@@ -137,8 +136,8 @@ int32_t BOTAN_PUBLIC_API(2,0) jacobi(const BigInt& a,
 * @return (b^x) % m
 */
 BigInt BOTAN_PUBLIC_API(2,0) power_mod(const BigInt& b,
-                           const BigInt& x,
-                           const BigInt& m);
+                                       const BigInt& x,
+                                       const BigInt& m);
 
 /**
 * Compute the square root of x modulo a prime using the
@@ -175,9 +174,18 @@ size_t BOTAN_PUBLIC_API(2,0) low_zero_bits(const BigInt& x);
 */
 bool BOTAN_PUBLIC_API(2,0) is_prime(const BigInt& n,
                                     RandomNumberGenerator& rng,
-                                    size_t prob = 56,
+                                    size_t prob = 64,
                                     bool is_random = false);
 
+/**
+* Test if the positive integer x is a perfect square ie if there
+* exists some positive integer y st y*y == x
+* See FIPS 186-4 sec C.4
+* @return 0 if the integer is not a perfect square, otherwise
+*         returns the positive y st y*y == x
+*/
+BigInt BOTAN_PUBLIC_API(2,8) is_perfect_square(const BigInt& x);
+
 inline bool quick_check_prime(const BigInt& n, RandomNumberGenerator& rng)
    { return is_prime(n, rng, 32); }
 
@@ -187,7 +195,6 @@ inline bool check_prime(const BigInt& n, RandomNumberGenerator& rng)
 inline bool verify_prime(const BigInt& n, RandomNumberGenerator& rng)
    { return is_prime(n, rng, 80); }
 
-
 /**
 * Randomly generate a prime suitable for discrete logarithm parameters
 * @param rng a random number generator
diff --git a/src/lib/math/numbertheory/primality.cpp b/src/lib/math/numbertheory/primality.cpp
new file mode 100644
index 000000000..5e683c1ff
--- /dev/null
+++ b/src/lib/math/numbertheory/primality.cpp
@@ -0,0 +1,206 @@
+/*
+* (C) 2016,2018 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include <botan/internal/primality.h>
+#include <botan/internal/monty_exp.h>
+#include <botan/bigint.h>
+#include <botan/monty.h>
+#include <botan/reducer.h>
+#include <botan/rng.h>
+#include <algorithm>
+
+namespace Botan {
+
+bool is_lucas_probable_prime(const BigInt& C, const Modular_Reducer& mod_C)
+   {
+   if(C <= 1)
+      return false;
+   else if(C == 2)
+      return true;
+   else if(C.is_even())
+      return false;
+   else if(C == 3 || C == 5 || C == 7 || C == 11 || C == 13)
+      return true;
+
+   BigInt D = 5;
+
+   for(;;)
+      {
+      int32_t j = jacobi(D, C);
+      if(j == 0)
+         return false;
+
+      if(j == -1)
+         break;
+
+      // Check 5, -7, 9, -11, 13, -15, 17, ...
+      if(D.is_negative())
+         {
+         D.flip_sign();
+         D += 2;
+         }
+      else
+         {
+         D += 2;
+         D.flip_sign();
+         }
+
+      if(D == 17 && is_perfect_square(C).is_nonzero())
+         return false;
+      }
+
+   const BigInt K = C + 1;
+   const size_t K_bits = K.bits() - 1;
+
+   BigInt U = 1;
+   BigInt V = 1;
+
+   BigInt Ut, Vt, U2, V2;
+
+   for(size_t i = 0; i != K_bits; ++i)
+      {
+      const uint8_t k_bit = K.get_bit(K_bits - 1 - i);
+
+      Ut = mod_C.multiply(U, V);
+
+      Vt = mod_C.reduce(mod_C.square(V) + mod_C.multiply(D, mod_C.square(U)));
+      if(Vt.is_odd())
+         Vt += C;
+      Vt >>= 1;
+      Vt = mod_C.reduce(Vt);
+
+      U = Ut;
+      V = Vt;
+
+      U2 = mod_C.reduce(Ut + Vt);
+      if(U2.is_odd())
+         U2 += C;
+      U2 >>= 1;
+
+      V2 = mod_C.reduce(Vt + Ut*D);
+      if(V2.is_odd())
+         V2 += C;
+      V2 >>= 1;
+
+      U.ct_cond_assign(k_bit, U2);
+      V.ct_cond_assign(k_bit, V2);
+      }
+
+   return (U == 0);
+   }
+
+bool is_bailie_psw_probable_prime(const BigInt& n, const Modular_Reducer& mod_n)
+   {
+   auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
+   return passes_miller_rabin_test(n, mod_n, monty_n, 2) && is_lucas_probable_prime(n, mod_n);
+   }
+
+bool is_bailie_psw_probable_prime(const BigInt& n)
+   {
+   Modular_Reducer mod_n(n);
+   return is_bailie_psw_probable_prime(n, mod_n);
+   }
+
+bool passes_miller_rabin_test(const BigInt& n,
+                              const Modular_Reducer& mod_n,
+                              const std::shared_ptr<Montgomery_Params>& monty_n,
+                              const BigInt& a)
+   {
+   BOTAN_ASSERT_NOMSG(n > 1);
+
+   const BigInt n_minus_1 = n - 1;
+   const size_t s = low_zero_bits(n_minus_1);
+   const BigInt nm1_s = n_minus_1 >> s;
+   const size_t n_bits = n.bits();
+
+   const size_t powm_window = 4;
+
+   auto powm_a_n = monty_precompute(monty_n, a, powm_window);
+
+   BigInt y = monty_execute(*powm_a_n, nm1_s, n_bits);
+
+   if(y == 1 || y == n_minus_1)
+      return true;
+
+   for(size_t i = 1; i != s; ++i)
+      {
+      y = mod_n.square(y);
+
+      if(y == 1) // found a non-trivial square root
+         return false;
+
+      /*
+      -1 is the trivial square root of unity, so ``a`` is not a
+      witness for this number - give up
+      */
+      if(y == n_minus_1)
+         return true;
+      }
+
+   return false;
+   }
+
+bool is_miller_rabin_probable_prime(const BigInt& n,
+                                    const Modular_Reducer& mod_n,
+                                    RandomNumberGenerator& rng,
+                                    size_t test_iterations)
+   {
+   BOTAN_ASSERT_NOMSG(n > 1);
+
+   auto monty_n = std::make_shared<Montgomery_Params>(n, mod_n);
+
+   for(size_t i = 0; i != test_iterations; ++i)
+      {
+      const BigInt a = BigInt::random_integer(rng, 2, n);
+
+      if(!passes_miller_rabin_test(n, mod_n, monty_n, a))
+         return false;
+      }
+
+   // Failed to find a counterexample
+   return true;
+   }
+
+
+size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random)
+   {
+   const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
+
+   /*
+   * If the candidate prime was maliciously constructed, we can't rely
+   * on arguments based on p being random.
+   */
+   if(random == false)
+      return base;
+
+   /*
+   * For randomly chosen numbers we can use the estimates from
+   * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
+   *
+   * These values are derived from the inequality for p(k,t) given on
+   * the second page.
+   */
+   if(prob <= 128)
+      {
+      if(n_bits >= 1536)
+         return 4; // < 2^-133
+      if(n_bits >= 1024)
+         return 6; // < 2^-133
+      if(n_bits >= 512)
+         return 12; // < 2^-129
+      if(n_bits >= 256)
+         return 29; // < 2^-128
+      }
+
+   /*
+   If the user desires a smaller error probability than we have
+   precomputed error estimates for, just fall back to using the worst
+   case error rate.
+   */
+   return base;
+   }
+
+}
diff --git a/src/lib/math/numbertheory/primality.h b/src/lib/math/numbertheory/primality.h
new file mode 100644
index 000000000..db7a76a74
--- /dev/null
+++ b/src/lib/math/numbertheory/primality.h
@@ -0,0 +1,100 @@
+/*
+* (C) 2018 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#ifndef BOTAN_PRIMALITY_TEST_H_
+#define BOTAN_PRIMALITY_TEST_H_
+
+#include <botan/types.h>
+#include <memory>
+
+namespace Botan {
+
+class BigInt;
+class Modular_Reducer;
+class Montgomery_Params;
+class RandomNumberGenerator;
+
+/**
+* Perform Lucas primality test
+* @see FIPS 186-4 C.3.3
+*
+* @warning it is possible to construct composite integers which pass
+* this test alone.
+*
+* @param n the positive integer to test
+* @param mod_n a pre-created Modular_Reducer for n
+* @return true if n seems probably prime, false if n is composite
+*/
+bool BOTAN_TEST_API is_lucas_probable_prime(const BigInt& n, const Modular_Reducer& mod_n);
+
+/**
+* Perform Bailie-PSW primality test
+*
+* This is a combination of Miller-Rabin with base 2 and a Lucas test. No known
+* composite integer passes both tests, though it is conjectured that infinitely
+* many composite counterexamples exist.
+*
+* @param n the positive integer to test
+* @param mod_n a pre-created Modular_Reducer for n
+* @return true if n seems probably prime, false if n is composite
+*/
+bool BOTAN_TEST_API is_bailie_psw_probable_prime(const BigInt& n, const Modular_Reducer& mod_n);
+
+/**
+* Perform Bailie-PSW primality test
+*
+* This is a combination of Miller-Rabin with base 2 and a Lucas test. No known
+* composite integer passes both tests, though it is conjectured that infinitely
+* many composite counterexamples exist.
+*
+* @param n the positive integer to test
+* @return true if n seems probably prime, false if n is composite
+*/
+bool is_bailie_psw_probable_prime(const BigInt& n);
+
+/**
+* Return required number of Miller-Rabin tests in order to
+* reach the specified probability of error.
+*
+* @param n_bits the bit-length of the integer being tested
+* @param prob chance of false positive is bounded by 1/2**prob
+* @param random is set if (and only if) the integer was randomly generated by us
+*        and thus cannot have been maliciously constructed.
+*/
+size_t miller_rabin_test_iterations(size_t n_bits, size_t prob, bool random);
+
+/**
+* Perform a single Miller-Rabin test with specified base
+*
+* @param n the positive integer to test
+* @param mod_n a pre-created Modular_Reducer for n
+* @param monty_n Montgomery parameters for n
+* @param a the base to check
+* @return result of primality test
+*/
+bool passes_miller_rabin_test(const BigInt& n,
+                              const Modular_Reducer& mod_n,
+                              const std::shared_ptr<Montgomery_Params>& monty_n,
+                              const BigInt& a);
+
+/**
+* Perform t iterations of a Miller-Rabin primality test with random bases
+*
+* @param n the positive integer to test
+* @param mod_n a pre-created Modular_Reducer for n
+* @param rng a random number generator
+* @param t number of tests to perform
+*
+* @return result of primality test
+*/
+bool BOTAN_TEST_API is_miller_rabin_probable_prime(const BigInt& n,
+                                                   const Modular_Reducer& mod_n,
+                                                   RandomNumberGenerator& rng,
+                                                   size_t t);
+
+}
+
+#endif
diff --git a/src/lib/math/numbertheory/reducer.cpp b/src/lib/math/numbertheory/reducer.cpp
index 98cf698ed..a5321c47c 100644
--- a/src/lib/math/numbertheory/reducer.cpp
+++ b/src/lib/math/numbertheory/reducer.cpp
@@ -32,11 +32,18 @@ Modular_Reducer::Modular_Reducer(const BigInt& mod)
       }
    }
 
-/*
-* Barrett Reduction
-*/
 BigInt Modular_Reducer::reduce(const BigInt& x) const
    {
+   BigInt r;
+   secure_vector<word> ws;
+   reduce(r, x, ws);
+   return r;
+   }
+
+void Modular_Reducer::reduce(BigInt& t1, const BigInt& x, secure_vector<word>& ws) const
+   {
+   if(&t1 == &x)
+      throw Invalid_State("Modular_Reducer arguments cannot alias");
    if(m_mod_words == 0)
       throw Invalid_State("Modular_Reducer: Never initalized");
 
@@ -45,12 +52,11 @@ BigInt Modular_Reducer::reduce(const BigInt& x) const
    if(x_sw >= (2*m_mod_words - 1) && x.cmp(m_modulus_2, false) >= 0)
       {
       // too big, fall back to normal division
-      return (x % m_modulus);
+      t1 = x % m_modulus;
+      return;
       }
 
-   secure_vector<word> ws;
-
-   BigInt t1 = x;
+   t1 = x;
    t1.set_sign(BigInt::Positive);
    t1 >>= (BOTAN_MP_WORD_BITS * (m_mod_words - 1));
 
@@ -83,8 +89,6 @@ BigInt Modular_Reducer::reduce(const BigInt& x) const
       {
       t1.rev_sub(m_modulus.data(), m_modulus.size(), ws);
       }
-
-   return t1;
    }
 
 }
diff --git a/src/lib/math/numbertheory/reducer.h b/src/lib/math/numbertheory/reducer.h
index c66c22034..5276adbbc 100644
--- a/src/lib/math/numbertheory/reducer.h
+++ b/src/lib/math/numbertheory/reducer.h
@@ -47,6 +47,8 @@ class BOTAN_PUBLIC_API(2,0) Modular_Reducer
       BigInt cube(const BigInt& x) const
          { return multiply(x, this->square(x)); }
 
+      void reduce(BigInt& out, const BigInt& x, secure_vector<word>& ws) const;
+
       bool initialized() const { return (m_mod_words != 0); }
 
       Modular_Reducer() { m_mod_words = 0; }
diff --git a/src/lib/pubkey/ec_group/ec_group.cpp b/src/lib/pubkey/ec_group/ec_group.cpp
index 586603507..30a0bb141 100644
--- a/src/lib/pubkey/ec_group/ec_group.cpp
+++ b/src/lib/pubkey/ec_group/ec_group.cpp
@@ -10,6 +10,7 @@
 
 #include <botan/ec_group.h>
 #include <botan/internal/point_mul.h>
+#include <botan/internal/primality.h>
 #include <botan/ber_dec.h>
 #include <botan/der_enc.h>
 #include <botan/oids.h>
@@ -19,12 +20,6 @@
 #include <botan/rng.h>
 #include <vector>
 
-#if defined(BOTAN_HAS_SYSTEM_RNG)
-   #include <botan/system_rng.h>
-#elif defined(BOTAN_HAS_HMAC_DRBG) && defined(BOTAN_HAS_SHA2_32)
-   #include <botan/hmac_drbg.h>
-#endif
-
 namespace Botan {
 
 class EC_Group_Data final
@@ -318,23 +313,7 @@ std::shared_ptr<EC_Group_Data> EC_Group::BER_decode_EC_group(const uint8_t bits[
          .end_cons()
          .verify_end();
 
-#if defined(BOTAN_HAS_SYSTEM_RNG)
-      System_RNG rng;
-#elif defined(BOTAN_HAS_HMAC_DRBG) && defined(BOTAN_HAS_SHA2_32)
-      /*
-      * This is not ideal because the data is attacker controlled, but
-      * it seems like it would be difficult for someone to come up
-      * with an valid ASN.1 encoding where the prime happened to pass
-      * Miller-Rabin test with exactly the values chosen when
-      * HMAC_DRBG is seeded with the overall data.
-      */
-      HMAC_DRBG rng("SHA-256");
-      rng.add_entropy(bits, len);
-#else
-      Null_RNG rng;
-#endif
-
-      if(p.bits() < 64 || p.is_negative() || (is_prime(p, rng) == false))
+      if(p.bits() < 64 || p.is_negative() || !is_bailie_psw_probable_prime(p))
          throw Decoding_Error("Invalid ECC p parameter");
 
       if(a.is_negative() || a >= p)
@@ -343,7 +322,7 @@ std::shared_ptr<EC_Group_Data> EC_Group::BER_decode_EC_group(const uint8_t bits[
       if(b <= 0 || b >= p)
          throw Decoding_Error("Invalid ECC b parameter");
 
-      if(order <= 0)
+      if(order <= 0 || !is_bailie_psw_probable_prime(order))
          throw Decoding_Error("Invalid ECC order parameter");
 
       if(cofactor <= 0 || cofactor >= 16)
diff --git a/src/tests/data/bn/isprime.vec b/src/tests/data/bn/isprime.vec
index ec6aca606..458208a1e 100644
--- a/src/tests/data/bn/isprime.vec
+++ b/src/tests/data/bn/isprime.vec
@@ -32,6 +32,18 @@ X = 0xFCEEE64D4D40D734058A51944F2B53152FFE7F15
 
 X = 0xEE23CE225FDEE2080403C2358C17A72D57C5B7CBE171D6D2BA59FE82DAABA9D3
 
+# Prime with low Hamming weight
+X = 0xFF4C000000000000000000000000000000000001
+
+# Mersenne primes
+X = 0x1ffff
+X = 0x7ffff
+X = 0x1fffffffffffffff
+X = 0x1ffffffffffffffffffffff
+X = 0x7ffffffffffffffffffffffffff
+X = 0x7fffffffffffffffffffffffffffffff
+X = 0x1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
+
 [NonPrime]
 X = 0
 X = 1
diff --git a/src/tests/data/bn/perfect_square.vec b/src/tests/data/bn/perfect_square.vec
new file mode 100644
index 000000000..b9821218d
--- /dev/null
+++ b/src/tests/data/bn/perfect_square.vec
@@ -0,0 +1,21 @@
+
+X = 2
+R = 0
+
+X = 4
+R = 2
+
+X = 8
+R = 0
+
+X = 9
+R = 3
+
+X = 16440582869613492769
+R = 4054698863
+
+X = 16440582861504095043
+R = 0
+
+X = 371930958274059827465211239444089
+R = 0
diff --git a/src/tests/test_bigint.cpp b/src/tests/test_bigint.cpp
index 8074a4f3b..146a9e8c9 100644
--- a/src/tests/test_bigint.cpp
+++ b/src/tests/test_bigint.cpp
@@ -9,6 +9,7 @@
 #if defined(BOTAN_HAS_NUMBERTHEORY)
    #include <botan/bigint.h>
    #include <botan/numthry.h>
+   #include <botan/internal/primality.h>
    #include <botan/reducer.h>
    #include <botan/pow_mod.h>
    #include <botan/parsing.h>
@@ -559,12 +560,32 @@ class BigInt_IsPrime_Test final : public Text_Based_Test
 
          Test::Result result("BigInt Test " + header);
          result.test_eq("is_prime", Botan::is_prime(value, Test::rng()), is_prime);
+
          return result;
          }
    };
 
 BOTAN_REGISTER_TEST("bn_isprime", BigInt_IsPrime_Test);
 
+class BigInt_IsSquare_Test final : public Text_Based_Test
+   {
+   public:
+      BigInt_IsSquare_Test() : Text_Based_Test("bn/perfect_square.vec", "X,R") {}
+
+      Test::Result run_one_test(const std::string&, const VarMap& vars) override
+         {
+         const BigInt value = vars.get_req_bn("X");
+         const BigInt expected = vars.get_req_bn("R");
+         const BigInt computed = Botan::is_perfect_square(value);
+
+         Test::Result result("BigInt IsSquare");
+         result.test_eq("is_perfect_square", computed, expected);
+         return result;
+         }
+   };
+
+BOTAN_REGISTER_TEST("bn_issquare", BigInt_IsSquare_Test);
+
 class BigInt_Ressol_Test final : public Text_Based_Test
    {
    public:
@@ -661,6 +682,47 @@ class BigInt_Rand_Test final : public Text_Based_Test
 
 BOTAN_REGISTER_TEST("bn_rand", BigInt_Rand_Test);
 
+class Lucas_Primality_Test final : public Test
+   {
+   public:
+
+      std::vector<Test::Result> run() override
+         {
+         const uint32_t lucas_max = (Test::run_long_tests() ? 100000 : 6000);
+
+         // OEIS A217120
+         std::set<uint32_t> lucas_pp{
+            323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179,
+            10877, 11419, 11663, 13919, 14839, 16109, 16211, 18407, 18971,
+            19043, 22499, 23407, 24569, 25199, 25877, 26069, 27323, 32759,
+            34943, 35207, 39059, 39203, 39689, 40309, 44099, 46979, 47879,
+            50183, 51983, 53663, 56279, 58519, 60377, 63881, 69509, 72389,
+            73919, 75077, 77219, 79547, 79799, 82983, 84419, 86063, 90287,
+            94667, 97019, 97439,
+         };
+
+         Test::Result result("Lucas primality test");
+
+         for(uint32_t i = 3; i <= lucas_max; i += 2)
+            {
+            Botan::Modular_Reducer mod_i(i);
+            const bool passes_lucas = Botan::is_lucas_probable_prime(i, mod_i);
+            const bool is_prime = Botan::is_prime(i, Test::rng());
+
+            const bool is_lucas_pp = (is_prime == false && passes_lucas == true);
+
+            if(is_lucas_pp)
+               result.confirm("Lucas pseudoprime is in list", lucas_pp.count(i) == 1);
+            else
+               result.confirm("Lucas non-pseudoprime is not in list", lucas_pp.count(i) == 0);
+            }
+
+         return {result};
+         }
+   };
+
+BOTAN_REGISTER_TEST("bn_lucas", Lucas_Primality_Test);
+
 class DSA_ParamGen_Test final : public Text_Based_Test
    {
    public: