Optimize Barrett reduction

OSS-Fuzz 6570 flagged an issue with slow modular exponentation.
It turned out the problem was not in the library version but the
simple square-and-multiply algorithm. Computing g^x % p with all
three integers being dense (high Hamming weight) numbers took about
1.5 seconds on a fast machine with almost all of the time taken
by the Barrett reductions. With these changes, same testcase
now takes only a tiny fraction of a second.

5 files changed, 134 insertions, 21 deletions

diff --git a/src/fuzzer/barrett.cpp b/src/fuzzer/barrett.cpp
new file mode 100644
index 000000000..1c5d88f87
--- /dev/null
+++ b/src/fuzzer/barrett.cpp
@@ -0,0 +1,41 @@
+/*
+* (C) 2018 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include "fuzzers.h"
+#include <botan/numthry.h>
+#include <botan/reducer.h>
+
+void fuzz(const uint8_t in[], size_t len)
+   {
+   static const size_t max_bits = 2048;
+
+   if(len % 2 != 0)
+      return;
+
+   const size_t part_size = len / 2;
+
+   if(part_size * 8 > max_bits)
+      return;
+
+   const Botan::BigInt x = Botan::BigInt::decode(in, part_size);
+   const Botan::BigInt p = Botan::BigInt::decode(in + part_size, part_size);
+
+   if(p.is_zero())
+      return;
+
+   const Botan::BigInt ref = x % p;
+
+   const Botan::Modular_Reducer mod_p(p);
+   const Botan::BigInt z = mod_p.reduce(x);
+
+   if(ref != z)
+      {
+      FUZZER_WRITE_AND_CRASH("X = " << x << "\n"
+                             << "P = " << p << "\n"
+                             << "Z = " << z << "\n"
+                             << "R = " << ref << "\n");
+      }
+   }
diff --git a/src/lib/math/bigint/big_ops2.cpp b/src/lib/math/bigint/big_ops2.cpp
index 7653884c9..fc6135c22 100644
--- a/src/lib/math/bigint/big_ops2.cpp
+++ b/src/lib/math/bigint/big_ops2.cpp
@@ -92,12 +92,54 @@ BigInt& BigInt::operator-=(const BigInt& y)
    return (*this);
    }
 
+BigInt& BigInt::rev_sub(const word y[], size_t y_sw, secure_vector<word>& ws)
+   {
+   /*
+   *this = BigInt(y, y_sw) - *this;
+   return *this;
+   */
+   if(this->sign() != BigInt::Positive)
+      throw Invalid_State("BigInt::sub_rev requires this is positive");
+
+   const size_t x_sw = this->sig_words();
+
+   const int32_t relative_size = bigint_cmp(y, y_sw, this->data(), x_sw);
+
+   ws.resize(std::max(y_sw, x_sw) + 1);
+   clear_mem(ws.data(), ws.size());
+
+   if(relative_size < 0)
+      {
+      bigint_sub3(ws.data(), this->data(), x_sw, y, y_sw);
+      this->flip_sign();
+      }
+   else if(relative_size == 0)
+      {
+      ws.clear();
+      }
+   else if(relative_size > 0)
+      {
+      bigint_sub3(ws.data(), y, y_sw, this->data(), x_sw);
+      }
+
+   m_reg.swap(ws);
+
+   return (*this);
+   }
+
 /*
 * Multiplication Operator
 */
 BigInt& BigInt::operator*=(const BigInt& y)
    {
-   const size_t x_sw = sig_words(), y_sw = y.sig_words();
+   secure_vector<word> ws;
+   return this->mul(y, ws);
+   }
+
+BigInt& BigInt::mul(const BigInt& y, secure_vector<word>& ws)
+   {
+   const size_t x_sw = sig_words();
+   const size_t y_sw = y.sig_words();
    set_sign((sign() == y.sign()) ? Positive : Negative);
 
    if(x_sw == 0 || y_sw == 0)
@@ -117,9 +159,16 @@ BigInt& BigInt::operator*=(const BigInt& y)
       }
    else
       {
-      grow_to(size() + y.size());
-      secure_vector<word> workspace(size());
-      bigint_mul(*this, BigInt(*this), y, workspace.data(), workspace.size());
+      const size_t new_size = x_sw + y_sw + 1;
+      ws.resize(new_size);
+      secure_vector<word> z_reg(new_size);
+
+      bigint_mul(z_reg.data(), z_reg.size(),
+                 data(), size(), x_sw,
+                 y.data(), y.size(), y_sw,
+                 ws.data(), ws.size());
+
+      z_reg.swap(m_reg);
       }
 
    return (*this);
diff --git a/src/lib/math/bigint/bigint.cpp b/src/lib/math/bigint/bigint.cpp
index c822a94e1..e99ddb50a 100644
--- a/src/lib/math/bigint/bigint.cpp
+++ b/src/lib/math/bigint/bigint.cpp
@@ -13,6 +13,11 @@
 
 namespace Botan {
 
+BigInt::BigInt(const word words[], size_t length)
+   {
+   m_reg.assign(words, words + length);
+   }
+
 /*
 * Construct a BigInt from a regular number
 */
diff --git a/src/lib/math/bigint/bigint.h b/src/lib/math/bigint/bigint.h
index 2c1d5e905..ca35bd07d 100644
--- a/src/lib/math/bigint/bigint.h
+++ b/src/lib/math/bigint/bigint.h
@@ -79,6 +79,13 @@ class BOTAN_PUBLIC_API(2,0) BigInt final
      BigInt(const uint8_t buf[], size_t length, Base base = Binary);
 
      /**
+     * Create a BigInt from an array of words
+     * @param words the words
+     * @param length number of words
+     */
+     BigInt(const word words[], size_t length);
+
+     /**
      * \brief Create a random BigInt of the specified size
      *
      * @param rng random number generator
@@ -222,6 +229,20 @@ class BOTAN_PUBLIC_API(2,0) BigInt final
      bool operator !() const { return (!is_nonzero()); }
 
      /**
+     * Multiply this with y
+     * @param y the BigInt to multiply with this
+     * @param ws a temp workspace
+     */
+     BigInt& mul(const BigInt& y, secure_vector<word>& ws);
+
+     /**
+     * Set *this to y - *this
+     * @param y the BigInt to subtract from
+     * @param ws a temp workspace
+     */
+     BigInt& rev_sub(const word y[], size_t y_size, secure_vector<word>& ws);
+
+     /**
      * Return *this below mod
      *
      * Assumes that *this is (if anything) only slightly larger than
diff --git a/src/lib/math/numbertheory/reducer.cpp b/src/lib/math/numbertheory/reducer.cpp
index d5f1666e1..1d7c2259a 100644
--- a/src/lib/math/numbertheory/reducer.cpp
+++ b/src/lib/math/numbertheory/reducer.cpp
@@ -34,7 +34,9 @@ BigInt Modular_Reducer::reduce(const BigInt& x) const
    if(m_mod_words == 0)
       throw Invalid_State("Modular_Reducer: Never initalized");
 
-   if(x.cmp(m_modulus, false) < 0)
+   const size_t x_sw = x.sig_words();
+
+   if(x_sw < m_mod_words || x.cmp(m_modulus, false) < 0)
       {
       if(x.is_negative())
          return x + m_modulus; // make positive
@@ -42,34 +44,29 @@ BigInt Modular_Reducer::reduce(const BigInt& x) const
       }
    else if(x.cmp(m_modulus_2, false) < 0)
       {
-      BigInt t1 = x;
-      t1.set_sign(BigInt::Positive);
-      t1 >>= (MP_WORD_BITS * (m_mod_words - 1));
-      t1 *= m_mu;
+      secure_vector<word> ws;
+
+      BigInt t1(x.data() + m_mod_words - 1, x_sw - (m_mod_words - 1));
 
+      t1.mul(m_mu, ws);
       t1 >>= (MP_WORD_BITS * (m_mod_words + 1));
-      t1 *= m_modulus;
 
+      t1.mul(m_modulus, ws);
       t1.mask_bits(MP_WORD_BITS * (m_mod_words + 1));
 
-      BigInt t2 = x;
-      t2.set_sign(BigInt::Positive);
-      t2.mask_bits(MP_WORD_BITS * (m_mod_words + 1));
-
-      t2 -= t1;
+      t1.rev_sub(x.data(), std::min(x_sw, m_mod_words + 1), ws);
 
-      if(t2.is_negative())
+      if(t1.is_negative())
          {
-         t2 += BigInt::power_of_2(MP_WORD_BITS * (m_mod_words + 1));
+         t1 += BigInt::power_of_2(MP_WORD_BITS * (m_mod_words + 1));
          }
 
-      while(t2 >= m_modulus)
-         t2 -= m_modulus;
+      t1.reduce_below(m_modulus, ws);
 
       if(x.is_positive())
-         return t2;
+         return t1;
       else
-         return (m_modulus - t2);
+         return (m_modulus - t1);
       }
    else
       {