Cleanup in number theory

Remove or hide deprecated functions

Consolidate some source files

17 files changed, 231 insertions, 939 deletions

diff --git a/src/lib/math/numbertheory/dsa_gen.cpp b/src/lib/math/numbertheory/dsa_gen.cpp
index a5efbc266..8a86d0cfd 100644
--- a/src/lib/math/numbertheory/dsa_gen.cpp
+++ b/src/lib/math/numbertheory/dsa_gen.cpp
@@ -5,6 +5,7 @@
 * Botan is released under the Simplified BSD License (see license.txt)
 */
 
+#include <botan/internal/primality.h>
 #include <botan/numthry.h>
 #include <botan/hash.h>
 #include <botan/reducer.h>
diff --git a/src/lib/math/numbertheory/info.txt b/src/lib/math/numbertheory/info.txt
index a5b914a30..3d91db8d7 100644
--- a/src/lib/math/numbertheory/info.txt
+++ b/src/lib/math/numbertheory/info.txt
@@ -1,5 +1,5 @@
 <defines>
-NUMBERTHEORY -> 20131128
+NUMBERTHEORY -> 20201108
 </defines>
 
 <header:public>
@@ -11,7 +11,6 @@ reducer.h
 curve_nistp.h
 monty.h
 monty_exp.h
-pow_mod.h
 primality.h
 </header:internal>
 
diff --git a/src/lib/math/numbertheory/jacobi.cpp b/src/lib/math/numbertheory/jacobi.cpp
deleted file mode 100644
index 284fc2b20..000000000
--- a/src/lib/math/numbertheory/jacobi.cpp
+++ /dev/null
@@ -1,52 +0,0 @@
-/*
-* Jacobi Function
-* (C) 1999-2007 Jack Lloyd
-*
-* Botan is released under the Simplified BSD License (see license.txt)
-*/
-
-#include <botan/numthry.h>
-
-namespace Botan {
-
-/*
-* Calculate the Jacobi symbol
-*/
-int32_t jacobi(const BigInt& a, const BigInt& n)
-   {
-   if(n.is_even() || n < 2)
-      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
-
-   BigInt x = a % n;
-   BigInt y = n;
-   int32_t J = 1;
-
-   while(y > 1)
-      {
-      x %= y;
-      if(x > y / 2)
-         {
-         x = y - x;
-         if(y % 4 == 3)
-            J = -J;
-         }
-      if(x.is_zero())
-         return 0;
-
-      size_t shifts = low_zero_bits(x);
-      x >>= shifts;
-      if(shifts % 2)
-         {
-         word y_mod_8 = y % 8;
-         if(y_mod_8 == 3 || y_mod_8 == 5)
-            J = -J;
-         }
-
-      if(x % 4 == 3 && y % 4 == 3)
-         J = -J;
-      std::swap(x, y);
-      }
-   return J;
-   }
-
-}
diff --git a/src/lib/math/numbertheory/make_prm.cpp b/src/lib/math/numbertheory/make_prm.cpp
index bac15c707..30f6ae91d 100644
--- a/src/lib/math/numbertheory/make_prm.cpp
+++ b/src/lib/math/numbertheory/make_prm.cpp
@@ -5,12 +5,12 @@
 * Botan is released under the Simplified BSD License (see license.txt)
 */
 
+#include <botan/internal/primality.h>
 #include <botan/numthry.h>
 #include <botan/rng.h>
 #include <botan/internal/bit_ops.h>
 #include <botan/internal/loadstor.h>
 #include <botan/reducer.h>
-#include <botan/internal/primality.h>
 #include <algorithm>
 
 namespace Botan {
diff --git a/src/lib/math/numbertheory/mod_inv.cpp b/src/lib/math/numbertheory/mod_inv.cpp
index b7c6f2b1e..f6f6c6c17 100644
--- a/src/lib/math/numbertheory/mod_inv.cpp
+++ b/src/lib/math/numbertheory/mod_inv.cpp
@@ -12,81 +12,6 @@
 
 namespace Botan {
 
-/*
-Sets result to a^-1 * 2^k mod a
-with n <= k <= 2n
-Returns k
-
-"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
-https://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.75.8377
-
-A const time implementation of this algorithm is described in
-"Constant Time Modular Inversion" Joppe W. Bos
-http://www.joppebos.com/files/CTInversion.pdf
-*/
-size_t almost_montgomery_inverse(BigInt& result,
-                                 const BigInt& a,
-                                 const BigInt& p)
-   {
-   size_t k = 0;
-
-   BigInt u = p, v = a, r = 0, s = 1;
-
-   while(v > 0)
-      {
-      if(u.is_even())
-         {
-         u >>= 1;
-         s <<= 1;
-         }
-      else if(v.is_even())
-         {
-         v >>= 1;
-         r <<= 1;
-         }
-      else if(u > v)
-         {
-         u -= v;
-         u >>= 1;
-         r += s;
-         s <<= 1;
-         }
-      else
-         {
-         v -= u;
-         v >>= 1;
-         s += r;
-         r <<= 1;
-         }
-
-      ++k;
-      }
-
-   if(r >= p)
-      {
-      r -= p;
-      }
-
-   result = p - r;
-
-   return k;
-   }
-
-BigInt normalized_montgomery_inverse(const BigInt& a, const BigInt& p)
-   {
-   BigInt r;
-   size_t k = almost_montgomery_inverse(r, a, p);
-
-   for(size_t i = 0; i != k; ++i)
-      {
-      if(r.is_odd())
-         r += p;
-      r >>= 1;
-      }
-
-   return r;
-   }
-
 namespace {
 
 BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
@@ -259,8 +184,8 @@ BigInt inverse_mod(const BigInt& n, const BigInt& mod)
    if(mod.is_odd())
       {
       /*
-      Fastpath for common case. This leaks information if n > mod
-      but we don't guarantee const time behavior in that case.
+      Fastpath for common case. This leaks if n is greater than mod or
+      not, but we don't guarantee const time behavior in that case.
       */
       if(n < mod)
          return inverse_mod_odd_modulus(n, mod);
@@ -313,44 +238,4 @@ BigInt inverse_mod(const BigInt& n, const BigInt& mod)
    return h;
    }
 
-// Deprecated forwarding functions:
-BigInt inverse_euclid(const BigInt& x, const BigInt& modulus)
-   {
-   return inverse_mod(x, modulus);
-   }
-
-BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
-   {
-   return inverse_mod_odd_modulus(n, mod);
-   }
-
-word monty_inverse(word a)
-   {
-   if(a % 2 == 0)
-      throw Invalid_Argument("monty_inverse only valid for odd integers");
-
-   /*
-   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
-   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
-   */
-
-   word b = 1;
-   word r = 0;
-
-   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
-      {
-      const word bi = b % 2;
-      r >>= 1;
-      r += bi << (BOTAN_MP_WORD_BITS - 1);
-
-      b -= a * bi;
-      b >>= 1;
-      }
-
-   // Now invert in addition space
-   r = (MP_WORD_MAX - r) + 1;
-
-   return r;
-   }
-
 }
diff --git a/src/lib/math/numbertheory/monty.cpp b/src/lib/math/numbertheory/monty.cpp
index 40cee2986..709f1ebe5 100644
--- a/src/lib/math/numbertheory/monty.cpp
+++ b/src/lib/math/numbertheory/monty.cpp
@@ -10,6 +10,35 @@
 
 namespace Botan {
 
+word monty_inverse(word a)
+   {
+   if(a % 2 == 0)
+      throw Invalid_Argument("monty_inverse only valid for odd integers");
+
+   /*
+   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
+   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
+   */
+
+   word b = 1;
+   word r = 0;
+
+   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
+      {
+      const word bi = b % 2;
+      r >>= 1;
+      r += bi << (BOTAN_MP_WORD_BITS - 1);
+
+      b -= a * bi;
+      b >>= 1;
+      }
+
+   // Now invert in addition space
+   r = (MP_WORD_MAX - r) + 1;
+
+   return r;
+   }
+
 Montgomery_Params::Montgomery_Params(const BigInt& p,
                                      const Modular_Reducer& mod_p)
    {
diff --git a/src/lib/math/numbertheory/monty.h b/src/lib/math/numbertheory/monty.h
index 8e0cd342f..b4c90397e 100644
--- a/src/lib/math/numbertheory/monty.h
+++ b/src/lib/math/numbertheory/monty.h
@@ -8,7 +8,6 @@
 #define BOTAN_MONTY_INT_H_
 
 #include <botan/bigint.h>
-BOTAN_FUTURE_INTERNAL_HEADER(monty.h)
 
 namespace Botan {
 
@@ -16,6 +15,13 @@ class Modular_Reducer;
 
 class Montgomery_Params;
 
+/*
+* Compute -input^-1 mod 2^MP_WORD_BITS. Throws an exception if input
+* is even. If input is odd, then input and 2^n are relatively prime
+* and an inverse exists.
+*/
+word monty_inverse(word input);
+
 /**
 * The Montgomery representation of an integer
 */
diff --git a/src/lib/math/numbertheory/mp_numth.cpp b/src/lib/math/numbertheory/mp_numth.cpp
deleted file mode 100644
index eef641996..000000000
--- a/src/lib/math/numbertheory/mp_numth.cpp
+++ /dev/null
@@ -1,84 +0,0 @@
-/*
-* Fused and Important MP Algorithms
-* (C) 1999-2007 Jack Lloyd
-*     2016 Matthias Gierlings
-*
-* Botan is released under the Simplified BSD License (see license.txt)
-*/
-
-#include <botan/numthry.h>
-#include <botan/internal/mp_core.h>
-#include <botan/internal/rounding.h>
-#include <algorithm>
-
-namespace Botan {
-
-/*
-* Square a BigInt
-*/
-BigInt square(const BigInt& x)
-   {
-   BigInt z = x;
-   secure_vector<word> ws;
-   z.square(ws);
-   return z;
-   }
-
-/*
-* Multiply-Add Operation
-*/
-BigInt mul_add(const BigInt& a, const BigInt& b, const BigInt& c)
-   {
-   if(c.is_negative())
-      throw Invalid_Argument("mul_add: Third argument must be > 0");
-
-   BigInt::Sign sign = BigInt::Positive;
-   if(a.sign() != b.sign())
-      sign = BigInt::Negative;
-
-   const size_t a_sw = a.sig_words();
-   const size_t b_sw = b.sig_words();
-   const size_t c_sw = c.sig_words();
-
-   BigInt r(sign, std::max(a_sw + b_sw, c_sw) + 1);
-   secure_vector<word> workspace(r.size());
-
-   bigint_mul(r.mutable_data(), r.size(),
-              a.data(), a.size(), a_sw,
-              b.data(), b.size(), b_sw,
-              workspace.data(), workspace.size());
-
-   const size_t r_size = std::max(r.sig_words(), c_sw);
-   bigint_add2(r.mutable_data(), r_size, c.data(), c_sw);
-   return r;
-   }
-
-/*
-* Subtract-Multiply Operation
-*/
-BigInt sub_mul(const BigInt& a, const BigInt& b, const BigInt& c)
-   {
-   if(a.is_negative() || b.is_negative())
-      throw Invalid_Argument("sub_mul: First two arguments must be >= 0");
-
-   BigInt r = a;
-   r -= b;
-   r *= c;
-   return r;
-   }
-
-/*
-* Multiply-Subtract Operation
-*/
-BigInt mul_sub(const BigInt& a, const BigInt& b, const BigInt& c)
-   {
-   if(c.is_negative() || c.is_zero())
-      throw Invalid_Argument("mul_sub: Third argument must be > 0");
-
-   BigInt r = a;
-   r *= b;
-   r -= c;
-   return r;
-   }
-
-}
diff --git a/src/lib/math/numbertheory/numthry.cpp b/src/lib/math/numbertheory/numthry.cpp
index f8e1ae8d9..18e372c58 100644
--- a/src/lib/math/numbertheory/numthry.cpp
+++ b/src/lib/math/numbertheory/numthry.cpp
@@ -1,6 +1,7 @@
 /*
 * Number Theory Functions
 * (C) 1999-2011,2016,2018,2019 Jack Lloyd
+* (C) 2007,2008 Falko Strenzke, FlexSecure GmbH
 *
 * Botan is released under the Simplified BSD License (see license.txt)
 */
@@ -34,6 +35,160 @@ void sub_abs(BigInt& z, const BigInt& x, const BigInt& y)
 }
 
 /*
+* Tonelli-Shanks algorithm
+*/
+BigInt ressol(const BigInt& a, const BigInt& p)
+   {
+   if(p <= 1 || p.is_even())
+      throw Invalid_Argument("ressol: invalid prime");
+
+   if(a == 0)
+      return 0;
+   else if(a < 0)
+      throw Invalid_Argument("ressol: value to solve for must be positive");
+   else if(a >= p)
+      throw Invalid_Argument("ressol: value to solve for must be less than p");
+
+   if(p == 2)
+      return a;
+
+   if(jacobi(a, p) != 1) // not a quadratic residue
+      return -BigInt(1);
+
+   if(p % 4 == 3) // The easy case
+      {
+      return power_mod(a, ((p+1) >> 2), p);
+      }
+
+   size_t s = low_zero_bits(p - 1);
+   BigInt q = p >> s;
+
+   q -= 1;
+   q >>= 1;
+
+   Modular_Reducer mod_p(p);
+
+   BigInt r = power_mod(a, q, p);
+   BigInt n = mod_p.multiply(a, mod_p.square(r));
+   r = mod_p.multiply(r, a);
+
+   if(n == 1)
+      return r;
+
+   // find random non quadratic residue z
+   BigInt z = 2;
+   while(jacobi(z, p) == 1) // while z quadratic residue
+      ++z;
+
+   BigInt c = power_mod(z, (q << 1) + 1, p);
+
+   while(n > 1)
+      {
+      q = n;
+
+      size_t i = 0;
+      while(q != 1)
+         {
+         q = mod_p.square(q);
+         ++i;
+
+         if(i >= s)
+            {
+            return -BigInt(1);
+            }
+         }
+
+      c = power_mod(c, BigInt::power_of_2(s-i-1), p);
+      r = mod_p.multiply(r, c);
+      c = mod_p.square(c);
+      n = mod_p.multiply(n, c);
+      s = i;
+      }
+
+   return r;
+   }
+
+/*
+* Calculate the Jacobi symbol
+*/
+int32_t jacobi(const BigInt& a, const BigInt& n)
+   {
+   if(n.is_even() || n < 2)
+      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
+
+   BigInt x = a % n;
+   BigInt y = n;
+   int32_t J = 1;
+
+   while(y > 1)
+      {
+      x %= y;
+      if(x > y / 2)
+         {
+         x = y - x;
+         if(y % 4 == 3)
+            J = -J;
+         }
+      if(x.is_zero())
+         return 0;
+
+      size_t shifts = low_zero_bits(x);
+      x >>= shifts;
+      if(shifts % 2)
+         {
+         word y_mod_8 = y % 8;
+         if(y_mod_8 == 3 || y_mod_8 == 5)
+            J = -J;
+         }
+
+      if(x % 4 == 3 && y % 4 == 3)
+         J = -J;
+      std::swap(x, y);
+      }
+   return J;
+   }
+
+/*
+* Multiply-Add Operation
+*/
+BigInt mul_add(const BigInt& a, const BigInt& b, const BigInt& c)
+   {
+   if(c.is_negative())
+      throw Invalid_Argument("mul_add: Third argument must be > 0");
+
+   BigInt::Sign sign = BigInt::Positive;
+   if(a.sign() != b.sign())
+      sign = BigInt::Negative;
+
+   const size_t a_sw = a.sig_words();
+   const size_t b_sw = b.sig_words();
+   const size_t c_sw = c.sig_words();
+
+   BigInt r(sign, std::max(a_sw + b_sw, c_sw) + 1);
+   secure_vector<word> workspace(r.size());
+
+   bigint_mul(r.mutable_data(), r.size(),
+              a.data(), a.size(), a_sw,
+              b.data(), b.size(), b_sw,
+              workspace.data(), workspace.size());
+
+   const size_t r_size = std::max(r.sig_words(), c_sw);
+   bigint_add2(r.mutable_data(), r_size, c.data(), c_sw);
+   return r;
+   }
+
+/*
+* Square a BigInt
+*/
+BigInt square(const BigInt& x)
+   {
+   BigInt z = x;
+   secure_vector<word> ws;
+   z.square(ws);
+   return z;
+   }
+
+/*
 * Return the number of 0 bits at the end of n
 */
 size_t low_zero_bits(const BigInt& n)
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index be9cd985e..4bf572f5b 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -27,30 +27,6 @@ BigInt BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Just use (a*b)+c")
            const BigInt& c);
 
 /**
-* Fused subtract-multiply
-* @param a an integer
-* @param b an integer
-* @param c an integer
-* @return (a-b)*c
-*/
-BigInt BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Just use (a-b)*c")
-   sub_mul(const BigInt& a,
-           const BigInt& b,
-           const BigInt& c);
-
-/**
-* Fused multiply-subtract
-* @param a an integer
-* @param b an integer
-* @param c an integer
-* @return (a*b)-c
-*/
-BigInt BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Just use (a*b)-c")
-   mul_sub(const BigInt& a,
-           const BigInt& b,
-           const BigInt& c);
-
-/**
 * Return the absolute value
 * @param n an integer
 * @return absolute value of n
@@ -95,36 +71,6 @@ BigInt BOTAN_PUBLIC_API(2,0) inverse_mod(const BigInt& x,
                                          const BigInt& modulus);
 
 /**
-* Deprecated modular inversion function. Use inverse_mod instead.
-* @param x a positive integer
-* @param modulus a positive integer
-* @return y st (x*y) % modulus == 1 or 0 if no such value
-*/
-BigInt BOTAN_DEPRECATED_API("Use inverse_mod") inverse_euclid(const BigInt& x, const BigInt& modulus);
-
-/**
-* Deprecated modular inversion function. Use inverse_mod instead.
-*/
-BigInt BOTAN_DEPRECATED_API("Use inverse_mod") ct_inverse_mod_odd_modulus(const BigInt& x, const BigInt& modulus);
-
-/**
-* Return a^-1 * 2^k mod b
-* Returns k, between n and 2n
-* Not const time
-*/
-size_t BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use inverse_mod")
-   almost_montgomery_inverse(BigInt& result,
-                             const BigInt& a,
-                             const BigInt& b);
-
-/**
-* Call almost_montgomery_inverse and correct the result to a^-1 mod b
-*/
-BigInt BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use inverse_mod")
-   normalized_montgomery_inverse(const BigInt& a, const BigInt& b);
-
-
-/**
 * Compute the Jacobi symbol. If n is prime, this is equivalent
 * to the Legendre symbol.
 * @see http://mathworld.wolfram.com/JacobiSymbol.html
@@ -156,14 +102,6 @@ BigInt BOTAN_PUBLIC_API(2,0) power_mod(const BigInt& b,
 */
 BigInt BOTAN_PUBLIC_API(2,0) ressol(const BigInt& x, const BigInt& p);
 
-/*
-* Compute -input^-1 mod 2^MP_WORD_BITS. Throws an exception if input
-* is even. If input is odd, then input and 2^n are relatively prime
-* and an inverse exists.
-*/
-word BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use inverse_mod")
-   monty_inverse(word input);
-
 /**
 * @param x an integer
 * @return count of the low zero bits in x, or, equivalently, the
@@ -194,18 +132,6 @@ bool BOTAN_PUBLIC_API(2,0) is_prime(const BigInt& n,
 */
 BigInt BOTAN_PUBLIC_API(2,8) is_perfect_square(const BigInt& x);
 
-inline bool BOTAN_DEPRECATED("Use is_prime")
-   quick_check_prime(const BigInt& n, RandomNumberGenerator& rng)
-   { return is_prime(n, rng, 32); }
-
-inline bool BOTAN_DEPRECATED("Use is_prime")
-   check_prime(const BigInt& n, RandomNumberGenerator& rng)
-   { return is_prime(n, rng, 56); }
-
-inline bool BOTAN_DEPRECATED("Use is_prime")
-   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
@@ -249,39 +175,6 @@ BigInt BOTAN_PUBLIC_API(2,0) random_safe_prime(RandomNumberGenerator& rng,
                                                size_t bits);
 
 /**
-* Generate DSA parameters using the FIPS 186 kosherizer
-* @param rng a random number generator
-* @param p_out where the prime p will be stored
-* @param q_out where the prime q will be stored
-* @param pbits how long p will be in bits
-* @param qbits how long q will be in bits
-* @return random seed used to generate this parameter set
-*/
-std::vector<uint8_t> BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use DL_Group")
-generate_dsa_primes(RandomNumberGenerator& rng,
-                    BigInt& p_out, BigInt& q_out,
-                    size_t pbits, size_t qbits);
-
-/**
-* Generate DSA parameters using the FIPS 186 kosherizer
-* @param rng a random number generator
-* @param p_out where the prime p will be stored
-* @param q_out where the prime q will be stored
-* @param pbits how long p will be in bits
-* @param qbits how long q will be in bits
-* @param seed the seed used to generate the parameters
-* @param offset optional offset from seed to start searching at
-* @return true if seed generated a valid DSA parameter set, otherwise
-          false. p_out and q_out are only valid if true was returned.
-*/
-bool BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use DL_Group")
-generate_dsa_primes(RandomNumberGenerator& rng,
-                    BigInt& p_out, BigInt& q_out,
-                    size_t pbits, size_t qbits,
-                    const std::vector<uint8_t>& seed,
-                    size_t offset = 0);
-
-/**
 * The size of the PRIMES[] array
 */
 const size_t PRIME_TABLE_SIZE = 6541;
diff --git a/src/lib/math/numbertheory/pow_mod.cpp b/src/lib/math/numbertheory/pow_mod.cpp
deleted file mode 100644
index 5e94fe586..000000000
--- a/src/lib/math/numbertheory/pow_mod.cpp
+++ /dev/null
@@ -1,328 +0,0 @@
-/*
-* Modular Exponentiation Proxy
-* (C) 1999-2007,2012,2018,2019 Jack Lloyd
-*     2016 Matthias Gierlings
-*
-* Botan is released under the Simplified BSD License (see license.txt)
-*/
-
-#include <botan/internal/pow_mod.h>
-#include <botan/numthry.h>
-#include <botan/reducer.h>
-#include <botan/internal/monty.h>
-#include <botan/internal/monty_exp.h>
-#include <botan/internal/rounding.h>
-#include <vector>
-
-namespace Botan {
-
-class Modular_Exponentiator
-   {
-   public:
-      virtual void set_base(const BigInt&) = 0;
-      virtual void set_exponent(const BigInt&) = 0;
-      virtual BigInt execute() const = 0;
-      virtual Modular_Exponentiator* copy() const = 0;
-
-      Modular_Exponentiator() = default;
-      Modular_Exponentiator(const Modular_Exponentiator&) = default;
-      Modular_Exponentiator & operator=(const Modular_Exponentiator&) = default;
-      virtual ~Modular_Exponentiator() = default;
-   };
-
-namespace {
-
-/**
-* Fixed Window Exponentiator
-*/
-class Fixed_Window_Exponentiator final : public Modular_Exponentiator
-   {
-   public:
-      void set_exponent(const BigInt& e) override { m_exp = e; }
-      void set_base(const BigInt&) override;
-      BigInt execute() const override;
-
-      Modular_Exponentiator* copy() const override
-         { return new Fixed_Window_Exponentiator(*this); }
-
-      Fixed_Window_Exponentiator(const BigInt&, Power_Mod::Usage_Hints);
-   private:
-      Modular_Reducer m_reducer;
-      BigInt m_exp;
-      size_t m_window_bits;
-      std::vector<BigInt> m_g;
-      Power_Mod::Usage_Hints m_hints;
-   };
-
-void Fixed_Window_Exponentiator::set_base(const BigInt& base)
-   {
-   m_window_bits = Power_Mod::window_bits(m_exp.bits(), base.bits(), m_hints);
-
-   m_g.resize(static_cast<size_t>(1) << m_window_bits);
-   m_g[0] = 1;
-   m_g[1] = m_reducer.reduce(base);
-
-   for(size_t i = 2; i != m_g.size(); ++i)
-      m_g[i] = m_reducer.multiply(m_g[i-1], m_g[1]);
-   }
-
-BigInt Fixed_Window_Exponentiator::execute() const
-   {
-   const size_t exp_nibbles = (m_exp.bits() + m_window_bits - 1) / m_window_bits;
-
-   BigInt x = 1;
-
-   for(size_t i = exp_nibbles; i > 0; --i)
-      {
-      for(size_t j = 0; j != m_window_bits; ++j)
-         x = m_reducer.square(x);
-
-      const uint32_t nibble = m_exp.get_substring(m_window_bits*(i-1), m_window_bits);
-
-      // not const time:
-      x = m_reducer.multiply(x, m_g[nibble]);
-      }
-   return x;
-   }
-
-/*
-* Fixed_Window_Exponentiator Constructor
-*/
-Fixed_Window_Exponentiator::Fixed_Window_Exponentiator(const BigInt& n,
-                                                       Power_Mod::Usage_Hints hints)
-   : m_reducer{Modular_Reducer(n)}, m_exp{}, m_window_bits{}, m_g{}, m_hints{hints}
-   {}
-
-class Montgomery_Exponentiator final : public Modular_Exponentiator
-   {
-   public:
-      void set_exponent(const BigInt& e) override { m_e = e; }
-      void set_base(const BigInt&) override;
-      BigInt execute() const override;
-
-      Modular_Exponentiator* copy() const override
-         { return new Montgomery_Exponentiator(*this); }
-
-      Montgomery_Exponentiator(const BigInt&, Power_Mod::Usage_Hints);
-   private:
-      BigInt m_p;
-      Modular_Reducer m_mod_p;
-      std::shared_ptr<const Montgomery_Params> m_monty_params;
-      std::shared_ptr<const Montgomery_Exponentation_State> m_monty;
-
-      BigInt m_e;
-      Power_Mod::Usage_Hints m_hints;
-   };
-
-void Montgomery_Exponentiator::set_base(const BigInt& base)
-   {
-   size_t window_bits = Power_Mod::window_bits(m_e.bits(), base.bits(), m_hints);
-   m_monty = monty_precompute(m_monty_params, m_mod_p.reduce(base), window_bits);
-   }
-
-BigInt Montgomery_Exponentiator::execute() const
-   {
-   /*
-   This leaks size of e via loop iterations, not possible to fix without
-   breaking this API. Round up to avoid leaking fine details.
-   */
-   return monty_execute(*m_monty, m_e, round_up(m_e.bits(), 8));
-   }
-
-Montgomery_Exponentiator::Montgomery_Exponentiator(const BigInt& mod,
-                                                   Power_Mod::Usage_Hints hints) :
-   m_p(mod),
-   m_mod_p(mod),
-   m_monty_params(std::make_shared<Montgomery_Params>(m_p, m_mod_p)),
-   m_hints(hints)
-   {
-   }
-
-}
-
-/*
-* Power_Mod Constructor
-*/
-Power_Mod::Power_Mod(const BigInt& n, Usage_Hints hints, bool disable_monty)
-   {
-   set_modulus(n, hints, disable_monty);
-   }
-
-Power_Mod::~Power_Mod() { /* for ~unique_ptr */ }
-
-/*
-* Power_Mod Copy Constructor
-*/
-Power_Mod::Power_Mod(const Power_Mod& other)
-   {
-   if(other.m_core.get())
-      m_core.reset(other.m_core->copy());
-   }
-
-/*
-* Power_Mod Assignment Operator
-*/
-Power_Mod& Power_Mod::operator=(const Power_Mod& other)
-   {
-   if(this != &other)
-      {
-      if(other.m_core)
-         m_core.reset(other.m_core->copy());
-      else
-         m_core.reset();
-      }
-   return (*this);
-   }
-
-/*
-* Set the modulus
-*/
-void Power_Mod::set_modulus(const BigInt& n, Usage_Hints hints, bool disable_monty) const
-   {
-   // Allow set_modulus(0) to mean "drop old state"
-
-   m_core.reset();
-
-   if(n != 0)
-      {
-      if(n.is_odd() && disable_monty == false)
-         m_core.reset(new Montgomery_Exponentiator(n, hints));
-      else
-         m_core.reset(new Fixed_Window_Exponentiator(n, hints));
-      }
-   }
-
-/*
-* Set the base
-*/
-void Power_Mod::set_base(const BigInt& b) const
-   {
-   if(b.is_negative())
-      throw Invalid_Argument("Power_Mod::set_base: arg must be non-negative");
-
-   if(!m_core)
-      throw Internal_Error("Power_Mod::set_base: m_core was NULL");
-   m_core->set_base(b);
-   }
-
-/*
-* Set the exponent
-*/
-void Power_Mod::set_exponent(const BigInt& e) const
-   {
-   if(e.is_negative())
-      throw Invalid_Argument("Power_Mod::set_exponent: arg must be > 0");
-
-   if(!m_core)
-      throw Internal_Error("Power_Mod::set_exponent: m_core was NULL");
-   m_core->set_exponent(e);
-   }
-
-/*
-* Compute the result
-*/
-BigInt Power_Mod::execute() const
-   {
-   if(!m_core)
-      throw Internal_Error("Power_Mod::execute: m_core was NULL");
-   return m_core->execute();
-   }
-
-/*
-* Try to choose a good window size
-*/
-size_t Power_Mod::window_bits(size_t exp_bits, size_t,
-                              Power_Mod::Usage_Hints hints)
-   {
-   static const size_t wsize[][2] = {
-      { 1434, 7 },
-      {  539, 6 },
-      {  197, 4 },
-      {   70, 3 },
-      {   17, 2 },
-      {    0, 0 }
-   };
-
-   size_t window_bits = 1;
-
-   if(exp_bits)
-      {
-      for(size_t j = 0; wsize[j][0]; ++j)
-         {
-         if(exp_bits >= wsize[j][0])
-            {
-            window_bits += wsize[j][1];
-            break;
-            }
-         }
-      }
-
-   if(hints & Power_Mod::BASE_IS_FIXED)
-      window_bits += 2;
-   if(hints & Power_Mod::EXP_IS_LARGE)
-      ++window_bits;
-
-   return window_bits;
-   }
-
-namespace {
-
-/*
-* Choose potentially useful hints
-*/
-Power_Mod::Usage_Hints choose_base_hints(const BigInt& b, const BigInt& n)
-   {
-   if(b == 2)
-      return Power_Mod::Usage_Hints(Power_Mod::BASE_IS_2 |
-                                    Power_Mod::BASE_IS_SMALL);
-
-   const size_t b_bits = b.bits();
-   const size_t n_bits = n.bits();
-
-   if(b_bits < n_bits / 32)
-      return Power_Mod::BASE_IS_SMALL;
-   if(b_bits > n_bits / 4)
-      return Power_Mod::BASE_IS_LARGE;
-
-   return Power_Mod::NO_HINTS;
-   }
-
-/*
-* Choose potentially useful hints
-*/
-Power_Mod::Usage_Hints choose_exp_hints(const BigInt& e, const BigInt& n)
-   {
-   const size_t e_bits = e.bits();
-   const size_t n_bits = n.bits();
-
-   if(e_bits < n_bits / 32)
-      return Power_Mod::BASE_IS_SMALL;
-   if(e_bits > n_bits / 4)
-      return Power_Mod::BASE_IS_LARGE;
-   return Power_Mod::NO_HINTS;
-   }
-
-}
-
-/*
-* Fixed_Exponent_Power_Mod Constructor
-*/
-Fixed_Exponent_Power_Mod::Fixed_Exponent_Power_Mod(const BigInt& e,
-                                                   const BigInt& n,
-                                                   Usage_Hints hints) :
-   Power_Mod(n, Usage_Hints(hints | EXP_IS_FIXED | choose_exp_hints(e, n)))
-   {
-   set_exponent(e);
-   }
-
-/*
-* Fixed_Base_Power_Mod Constructor
-*/
-Fixed_Base_Power_Mod::Fixed_Base_Power_Mod(const BigInt& b, const BigInt& n,
-                                           Usage_Hints hints) :
-   Power_Mod(n, Usage_Hints(hints | BASE_IS_FIXED | choose_base_hints(b, n)))
-   {
-   set_base(b);
-   }
-
-}
diff --git a/src/lib/math/numbertheory/pow_mod.h b/src/lib/math/numbertheory/pow_mod.h
deleted file mode 100644
index 5581752ea..000000000
--- a/src/lib/math/numbertheory/pow_mod.h
+++ /dev/null
@@ -1,120 +0,0 @@
-/*
-* Modular Exponentiator
-* (C) 1999-2007 Jack Lloyd
-*
-* Botan is released under the Simplified BSD License (see license.txt)
-*/
-
-#ifndef BOTAN_POWER_MOD_H_
-#define BOTAN_POWER_MOD_H_
-
-#include <botan/bigint.h>
-
-namespace Botan {
-
-class Modular_Exponentiator;
-
-/**
-* Modular Exponentiator Proxy
-*/
-class BOTAN_TEST_API Power_Mod
-   {
-   public:
-
-      enum Usage_Hints {
-         NO_HINTS        = 0x0000,
-
-         BASE_IS_FIXED   = 0x0001,
-         BASE_IS_SMALL   = 0x0002,
-         BASE_IS_LARGE   = 0x0004,
-         BASE_IS_2       = 0x0008,
-
-         EXP_IS_FIXED    = 0x0100,
-         EXP_IS_SMALL    = 0x0200,
-         EXP_IS_LARGE    = 0x0400
-      };
-
-      /*
-      * Try to choose a good window size
-      */
-      static size_t window_bits(size_t exp_bits, size_t base_bits,
-                                Power_Mod::Usage_Hints hints);
-
-      /**
-      * @param modulus the modulus
-      * @param hints Passed to set_modulus if modulus > 0
-      * @param disable_montgomery_arith Disables use of Montgomery
-      * representation. Likely only useful for testing.
-      */
-      void set_modulus(const BigInt& modulus,
-                       Usage_Hints hints = NO_HINTS,
-                       bool disable_montgomery_arith = false) const;
-
-      /**
-      * Set the base
-      */
-      void set_base(const BigInt& base) const;
-
-      /**
-      * Set the exponent
-      */
-      void set_exponent(const BigInt& exponent) const;
-
-      /**
-      * All three of the above functions must have already been called.
-      * @return result of g^x%p
-      */
-      BigInt execute() const;
-
-      Power_Mod& operator=(const Power_Mod&);
-
-      /**
-      * @param modulus Optionally call set_modulus
-      * @param hints Passed to set_modulus if modulus > 0
-      * @param disable_montgomery_arith Disables use of Montgomery
-      * representation. Likely only useful for testing.
-      */
-      Power_Mod(const BigInt& modulus = 0,
-                Usage_Hints hints = NO_HINTS,
-                bool disable_montgomery_arith = false);
-      Power_Mod(const Power_Mod&);
-      virtual ~Power_Mod();
-   private:
-      mutable std::unique_ptr<Modular_Exponentiator> m_core;
-   };
-
-/**
-* Fixed Exponent Modular Exponentiator Proxy
-*/
-class Fixed_Exponent_Power_Mod final : public Power_Mod
-   {
-   public:
-      BigInt operator()(const BigInt& b) const
-         { set_base(b); return execute(); }
-
-      Fixed_Exponent_Power_Mod() = default;
-
-      Fixed_Exponent_Power_Mod(const BigInt& exponent,
-                               const BigInt& modulus,
-                               Usage_Hints hints = NO_HINTS);
-   };
-
-/**
-* Fixed Base Modular Exponentiator Proxy
-*/
-class Fixed_Base_Power_Mod final : public Power_Mod
-   {
-   public:
-      BigInt operator()(const BigInt& e) const
-         { set_exponent(e); return execute(); }
-
-      Fixed_Base_Power_Mod() = default;
-
-      Fixed_Base_Power_Mod(const BigInt& base,
-                           const BigInt& modulus,
-                           Usage_Hints hints = NO_HINTS);
-   };
-
-}
-
-#endif
diff --git a/src/lib/math/numbertheory/primality.h b/src/lib/math/numbertheory/primality.h
index db7a76a74..8a602bbad 100644
--- a/src/lib/math/numbertheory/primality.h
+++ b/src/lib/math/numbertheory/primality.h
@@ -9,6 +9,7 @@
 
 #include <botan/types.h>
 #include <memory>
+#include <vector>
 
 namespace Botan {
 
@@ -95,6 +96,38 @@ bool BOTAN_TEST_API is_miller_rabin_probable_prime(const BigInt& n,
                                                    RandomNumberGenerator& rng,
                                                    size_t t);
 
+/**
+* Generate DSA parameters using the FIPS 186 kosherizer
+* @param rng a random number generator
+* @param p_out where the prime p will be stored
+* @param q_out where the prime q will be stored
+* @param pbits how long p will be in bits
+* @param qbits how long q will be in bits
+* @return random seed used to generate this parameter set
+*/
+std::vector<uint8_t>
+generate_dsa_primes(RandomNumberGenerator& rng,
+                    BigInt& p_out, BigInt& q_out,
+                    size_t pbits, size_t qbits);
+
+/**
+* Generate DSA parameters using the FIPS 186 kosherizer
+* @param rng a random number generator
+* @param p_out where the prime p will be stored
+* @param q_out where the prime q will be stored
+* @param pbits how long p will be in bits
+* @param qbits how long q will be in bits
+* @param seed the seed used to generate the parameters
+* @param offset optional offset from seed to start searching at
+* @return true if seed generated a valid DSA parameter set, otherwise
+          false. p_out and q_out are only valid if true was returned.
+*/
+bool BOTAN_TEST_API generate_dsa_primes(RandomNumberGenerator& rng,
+                                        BigInt& p_out, BigInt& q_out,
+                                        size_t pbits, size_t qbits,
+                                        const std::vector<uint8_t>& seed,
+                                        size_t offset = 0);
+
 }
 
 #endif
diff --git a/src/lib/math/numbertheory/ressol.cpp b/src/lib/math/numbertheory/ressol.cpp
deleted file mode 100644
index de89e0384..000000000
--- a/src/lib/math/numbertheory/ressol.cpp
+++ /dev/null
@@ -1,87 +0,0 @@
-/*
-* (C) 2007,2008 Falko Strenzke, FlexSecure GmbH
-* (C) 2008 Jack Lloyd
-*
-* Botan is released under the Simplified BSD License (see license.txt)
-*/
-
-#include <botan/numthry.h>
-#include <botan/reducer.h>
-
-namespace Botan {
-
-/*
-* Tonelli-Shanks algorithm
-*/
-BigInt ressol(const BigInt& a, const BigInt& p)
-   {
-   if(p <= 1 || p.is_even())
-      throw Invalid_Argument("ressol: invalid prime");
-
-   if(a == 0)
-      return 0;
-   else if(a < 0)
-      throw Invalid_Argument("ressol: value to solve for must be positive");
-   else if(a >= p)
-      throw Invalid_Argument("ressol: value to solve for must be less than p");
-
-   if(p == 2)
-      return a;
-
-   if(jacobi(a, p) != 1) // not a quadratic residue
-      return -BigInt(1);
-
-   if(p % 4 == 3) // The easy case
-      {
-      return power_mod(a, ((p+1) >> 2), p);
-      }
-
-   size_t s = low_zero_bits(p - 1);
-   BigInt q = p >> s;
-
-   q -= 1;
-   q >>= 1;
-
-   Modular_Reducer mod_p(p);
-
-   BigInt r = power_mod(a, q, p);
-   BigInt n = mod_p.multiply(a, mod_p.square(r));
-   r = mod_p.multiply(r, a);
-
-   if(n == 1)
-      return r;
-
-   // find random non quadratic residue z
-   BigInt z = 2;
-   while(jacobi(z, p) == 1) // while z quadratic residue
-      ++z;
-
-   BigInt c = power_mod(z, (q << 1) + 1, p);
-
-   while(n > 1)
-      {
-      q = n;
-
-      size_t i = 0;
-      while(q != 1)
-         {
-         q = mod_p.square(q);
-         ++i;
-
-         if(i >= s)
-            {
-            return -BigInt(1);
-            }
-         }
-
-      c = power_mod(c, BigInt::power_of_2(s-i-1), p);
-      r = mod_p.multiply(r, c);
-      c = mod_p.square(c);
-      n = mod_p.multiply(n, c);
-      s = i;
-      }
-
-   return r;
-   }
-
-}
diff --git a/src/lib/pubkey/dl_group/dl_group.cpp b/src/lib/pubkey/dl_group/dl_group.cpp
index f147c33fb..0de5a5de4 100644
--- a/src/lib/pubkey/dl_group/dl_group.cpp
+++ b/src/lib/pubkey/dl_group/dl_group.cpp
@@ -8,6 +8,7 @@
 #include <botan/dl_group.h>
 #include <botan/numthry.h>
 #include <botan/reducer.h>
+#include <botan/internal/primality.h>
 #include <botan/internal/monty.h>
 #include <botan/internal/divide.h>
 #include <botan/der_enc.h>
diff --git a/src/lib/pubkey/ec_group/curve_gfp.cpp b/src/lib/pubkey/ec_group/curve_gfp.cpp
index da049aebb..c4ac23e98 100644
--- a/src/lib/pubkey/ec_group/curve_gfp.cpp
+++ b/src/lib/pubkey/ec_group/curve_gfp.cpp
@@ -12,6 +12,7 @@
 #include <botan/reducer.h>
 #include <botan/internal/mp_core.h>
 #include <botan/internal/mp_asmi.h>
+#include <botan/internal/monty.h>
 
 namespace Botan {
 
diff --git a/src/tests/test_bigint.cpp b/src/tests/test_bigint.cpp
index cbca5667e..1552fad08 100644
--- a/src/tests/test_bigint.cpp
+++ b/src/tests/test_bigint.cpp
@@ -12,7 +12,6 @@
    #include <botan/internal/divide.h>
    #include <botan/internal/primality.h>
    #include <botan/reducer.h>
-   #include <botan/internal/pow_mod.h>
    #include <botan/internal/parsing.h>
    #include "test_rng.h"
 #endif
@@ -535,38 +534,6 @@ class BigInt_Powmod_Test final : public Text_Based_Test
          const BigInt expected = vars.get_req_bn("Output");
 
          result.test_eq("power_mod", Botan::power_mod(base, exponent, modulus), expected);
-
-         /*
-         * Only the basic power_mod interface supports negative base
-         */
-         if(base.is_negative())
-            return result;
-
-         Botan::Power_Mod pow_mod1(modulus);
-
-         pow_mod1.set_base(base);
-         pow_mod1.set_exponent(exponent);
-         result.test_eq("pow_mod1", pow_mod1.execute(), expected);
-
-         Botan::Power_Mod pow_mod2(modulus);
-
-         // Reverses ordering which affects window size
-         pow_mod2.set_exponent(exponent);
-         pow_mod2.set_base(base);
-         result.test_eq("pow_mod2", pow_mod2.execute(), expected);
-         result.test_eq("pow_mod2 #2", pow_mod2.execute(), expected);
-
-         if(modulus.is_odd())
-            {
-            // TODO: test different hints
-            // also TODO: remove bogus hinting arguments :)
-            Botan::Power_Mod pow_mod3(modulus, Botan::Power_Mod::NO_HINTS, /*disable_montgomery=*/true);
-
-            pow_mod3.set_exponent(exponent);
-            pow_mod3.set_base(base);
-            result.test_eq("pow_mod_fixed_window", pow_mod3.execute(), expected);
-            }
-
          return result;
          }
    };
@@ -673,13 +640,6 @@ class BigInt_InvMod_Test final : public Text_Based_Test
             }
          */
 
-         if(mod.is_odd() && a_inv != 0)
-            {
-            result.test_eq("normalized_montgomery_inverse",
-                           normalized_montgomery_inverse(a, mod),
-                           expected);
-            }
-
          return result;
          }
    };