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.

3 files changed, 79 insertions, 4 deletions

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