Add Lucas test from FIPS 186-4

This eliminates an issue identified in the paper
"Prime and Prejudice: Primality Testing Under Adversarial Conditions"
by Albrecht, Massimo, Paterson and Somorovsky

where DL_Group::verify_group with strong=false would accept a composite
q with probability 1/4096, which is exactly as the error bound is
documented, but still unfortunate.

11 files changed, 403 insertions, 123 deletions

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; }