Fix exponentiation bug, related fixes

GH #754 exposed a bug in the non-Montgomery exponentiation case.
It turned out then when the fixed window was picked to any value
> 1, the result would be incorrect due to an off by one. This is
the one line fix in powm_fw.cpp

Also fix a bug in bigint_mul which caused incorrect results,
because the output BigInt was not being zeroed out before use. This
is only exposed in rare cases, found (somewhat indirectly) in
OSS-Fuzz #287.

Add more modular exponentiation tests, which would have caught
these issues earlier.

10 files changed, 134 insertions, 37 deletions

diff --git a/src/lib/math/mp/mp_karat.cpp b/src/lib/math/mp/mp_karat.cpp
index 7a763e2a9..62a52b88c 100644
--- a/src/lib/math/mp/mp_karat.cpp
+++ b/src/lib/math/mp/mp_karat.cpp
@@ -255,53 +255,58 @@ size_t karatsuba_size(size_t z_size, size_t x_size, size_t x_sw)
 */
 void bigint_mul(BigInt& z, const BigInt& x, const BigInt& y, word workspace[])
    {
-   if(x.sig_words() == 1)
+   const size_t x_sig_words = x.sig_words();
+   const size_t y_sig_words = y.sig_words();
+
+   clear_mem(z.mutable_data(), z.size());
+
+   if(x_sig_words == 1)
       {
-      bigint_linmul3(z.mutable_data(), y.data(), y.sig_words(), x.data()[0]);
+      bigint_linmul3(z.mutable_data(), y.data(), y_sig_words, x.data()[0]);
       }
-   else if(y.sig_words() == 1)
+   else if(y_sig_words == 1)
       {
-      bigint_linmul3(z.mutable_data(), x.data(), x.sig_words(), y.data()[0]);
+      bigint_linmul3(z.mutable_data(), x.data(), x_sig_words, y.data()[0]);
       }
-   else if(x.sig_words() <= 4 && x.size() >= 4 &&
-           y.sig_words() <= 4 && y.size() >= 4 && z.size() >= 8)
+   else if(x_sig_words <= 4 && x.size() >= 4 &&
+           y_sig_words <= 4 && y.size() >= 4 && z.size() >= 8)
       {
       bigint_comba_mul4(z.mutable_data(), x.data(), y.data());
       }
-   else if(x.sig_words() <= 6 && x.size() >= 6 &&
-           y.sig_words() <= 6 && y.size() >= 6 && z.size() >= 12)
+   else if(x_sig_words <= 6 && x.size() >= 6 &&
+           y_sig_words <= 6 && y.size() >= 6 && z.size() >= 12)
       {
       bigint_comba_mul6(z.mutable_data(), x.data(), y.data());
       }
-   else if(x.sig_words() <= 8 && x.size() >= 8 &&
-           y.sig_words() <= 8 && y.size() >= 8 && z.size() >= 16)
+   else if(x_sig_words <= 8 && x.size() >= 8 &&
+           y_sig_words <= 8 && y.size() >= 8 && z.size() >= 16)
       {
       bigint_comba_mul8(z.mutable_data(), x.data(), y.data());
       }
-   else if(x.sig_words() <= 9 && x.size() >= 9 &&
-           y.sig_words() <= 9 && y.size() >= 9 && z.size() >= 18)
+   else if(x_sig_words <= 9 && x.size() >= 9 &&
+           y_sig_words <= 9 && y.size() >= 9 && z.size() >= 18)
       {
       bigint_comba_mul9(z.mutable_data(), x.data(), y.data());
       }
-   else if(x.sig_words() <= 16 && x.size() >= 16 &&
-           y.sig_words() <= 16 && y.size() >= 16 && z.size() >= 32)
+   else if(x_sig_words <= 16 && x.size() >= 16 &&
+           y_sig_words <= 16 && y.size() >= 16 && z.size() >= 32)
       {
       bigint_comba_mul16(z.mutable_data(), x.data(), y.data());
       }
-   else if(x.sig_words() < KARATSUBA_MULTIPLY_THRESHOLD ||
-           y.sig_words() < KARATSUBA_MULTIPLY_THRESHOLD ||
+   else if(x_sig_words < KARATSUBA_MULTIPLY_THRESHOLD ||
+           y_sig_words < KARATSUBA_MULTIPLY_THRESHOLD ||
            !workspace)
       {
-      basecase_mul(z.mutable_data(), x.data(), x.sig_words(), y.data(), y.sig_words());
+      basecase_mul(z.mutable_data(), x.data(), x_sig_words, y.data(), y_sig_words);
       }
    else
       {
-      const size_t N = karatsuba_size(z.size(), x.size(), x.sig_words(), y.size(), y.sig_words());
+      const size_t N = karatsuba_size(z.size(), x.size(), x_sig_words, y.size(), y_sig_words);
 
       if(N)
          karatsuba_mul(z.mutable_data(), x.data(), y.data(), N, workspace);
       else
-         basecase_mul(z.mutable_data(), x.data(), x.sig_words(), y.data(), y.sig_words());
+         basecase_mul(z.mutable_data(), x.data(), x_sig_words, y.data(), y_sig_words);
       }
    }
 
diff --git a/src/lib/math/numbertheory/def_powm.h b/src/lib/math/numbertheory/def_powm.h
index d60ca8173..10ae8aa5b 100644
--- a/src/lib/math/numbertheory/def_powm.h
+++ b/src/lib/math/numbertheory/def_powm.h
@@ -52,6 +52,7 @@ class Montgomery_Exponentiator : public Modular_Exponentiator
       Montgomery_Exponentiator(const BigInt&, Power_Mod::Usage_Hints);
    private:
       BigInt m_exp, m_modulus, m_R_mod, m_R2_mod;
+      Modular_Reducer m_reducer;
       word m_mod_prime;
       size_t m_mod_words, m_exp_bits, m_window_bits;
       Power_Mod::Usage_Hints m_hints;
diff --git a/src/lib/math/numbertheory/numthry.cpp b/src/lib/math/numbertheory/numthry.cpp
index 6c3d2c931..71dbf6aba 100644
--- a/src/lib/math/numbertheory/numthry.cpp
+++ b/src/lib/math/numbertheory/numthry.cpp
@@ -358,6 +358,9 @@ word monty_inverse(word input)
       y1 = y;
       }
 
+   const word check = y2 * input;
+   BOTAN_ASSERT_EQUAL(check, 1, "monty_inverse result is inverse of input");
+
    // Now invert in addition space
    y2 = (MP_WORD_MAX - y2) + 1;
 
diff --git a/src/lib/math/numbertheory/pow_mod.cpp b/src/lib/math/numbertheory/pow_mod.cpp
index 5503f313c..9f7459ac5 100644
--- a/src/lib/math/numbertheory/pow_mod.cpp
+++ b/src/lib/math/numbertheory/pow_mod.cpp
@@ -13,10 +13,10 @@ namespace Botan {
 /*
 * Power_Mod Constructor
 */
-Power_Mod::Power_Mod(const BigInt& n, Usage_Hints hints)
+Power_Mod::Power_Mod(const BigInt& n, Usage_Hints hints, bool disable_monty)
    {
    m_core = nullptr;
-   set_modulus(n, hints);
+   set_modulus(n, hints, disable_monty);
    }
 
 /*
@@ -58,14 +58,21 @@ Power_Mod::~Power_Mod()
 /*
 * Set the modulus
 */
-void Power_Mod::set_modulus(const BigInt& n, Usage_Hints hints) const
+void Power_Mod::set_modulus(const BigInt& n, Usage_Hints hints, bool disable_monty) const
    {
+   // Allow set_modulus(0) to mean "drop old state"
+
    delete m_core;
+   m_core = nullptr;
+   disable_monty =true;
+   if(n != 0)
+      {
+      if(n.is_odd() && disable_monty == false)
+         m_core = new Montgomery_Exponentiator(n, hints);
 
-   if(n.is_odd())
-      m_core = new Montgomery_Exponentiator(n, hints);
-   else if(n != 0)
-      m_core = new Fixed_Window_Exponentiator(n, hints);
+      if(!m_core)
+         m_core = new Fixed_Window_Exponentiator(n, hints);
+      }
    }
 
 /*
diff --git a/src/lib/math/numbertheory/pow_mod.h b/src/lib/math/numbertheory/pow_mod.h
index 4f94fd62d..194d00ec8 100644
--- a/src/lib/math/numbertheory/pow_mod.h
+++ b/src/lib/math/numbertheory/pow_mod.h
@@ -51,15 +51,43 @@ class BOTAN_DLL Power_Mod
       static size_t window_bits(size_t exp_bits, size_t base_bits,
                                 Power_Mod::Usage_Hints hints);
 
-      void set_modulus(const BigInt&, Usage_Hints = NO_HINTS) const;
-      void set_base(const BigInt&) const;
-      void set_exponent(const BigInt&) const;
+      /**
+      * @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 = 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&);
 
-      Power_Mod(const BigInt& = 0, Usage_Hints = NO_HINTS);
+      /**
+      * @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:
diff --git a/src/lib/math/numbertheory/powm_fw.cpp b/src/lib/math/numbertheory/powm_fw.cpp
index 7369959a9..7d69a2602 100644
--- a/src/lib/math/numbertheory/powm_fw.cpp
+++ b/src/lib/math/numbertheory/powm_fw.cpp
@@ -31,7 +31,7 @@ void Fixed_Window_Exponentiator::set_base(const BigInt& base)
    m_g[1] = base;
 
    for(size_t i = 2; i != m_g.size(); ++i)
-      m_g[i] = m_reducer.multiply(m_g[i-1], m_g[0]);
+      m_g[i] = m_reducer.multiply(m_g[i-1], m_g[1]);
    }
 
 /*
diff --git a/src/lib/math/numbertheory/powm_mnt.cpp b/src/lib/math/numbertheory/powm_mnt.cpp
index 572f0de98..546a2739a 100644
--- a/src/lib/math/numbertheory/powm_mnt.cpp
+++ b/src/lib/math/numbertheory/powm_mnt.cpp
@@ -41,7 +41,7 @@ void Montgomery_Exponentiator::set_base(const BigInt& base)
                     workspace.data());
    m_g[0] = z;
 
-   m_g[1] = (base >= m_modulus) ? (base % m_modulus) : base;
+   m_g[1] = m_reducer.reduce(base);
 
    bigint_monty_mul(z, m_g[1], m_R2_mod,
                     m_modulus.data(), m_mod_words, m_mod_prime,
@@ -112,6 +112,7 @@ BigInt Montgomery_Exponentiator::execute() const
 Montgomery_Exponentiator::Montgomery_Exponentiator(const BigInt& mod,
                                                    Power_Mod::Usage_Hints hints) :
    m_modulus(mod),
+   m_reducer(m_modulus),
    m_mod_words(m_modulus.sig_words()),
    m_window_bits(1),
    m_hints(hints)
@@ -123,8 +124,8 @@ Montgomery_Exponentiator::Montgomery_Exponentiator(const BigInt& mod,
    m_mod_prime = monty_inverse(mod.word_at(0));
 
    const BigInt r = BigInt::power_of_2(m_mod_words * BOTAN_MP_WORD_BITS);
-   m_R_mod = r % m_modulus;
-   m_R2_mod = (m_R_mod * m_R_mod) % m_modulus;
+   m_R_mod = m_reducer.reduce(r);
+   m_R2_mod = m_reducer.square(m_R_mod);
    m_exp_bits = 0;
    }
 
diff --git a/src/tests/data/bn/mul.vec b/src/tests/data/bn/mul.vec
index d70428112..a7e5601c7 100644
--- a/src/tests/data/bn/mul.vec
+++ b/src/tests/data/bn/mul.vec
@@ -410,3 +410,7 @@ Output = 0xED478AFD970F5D8E22B1E85E2C186CD98172870A148C78475D57F7B524BD7752FD6F7
 In1 = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
 In2 = 0x4000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001
 Output = 0x4000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000BFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
+
+In1 = 1024
+In2 = 0x918100000000000181fe7e0000000000fffe0001000
+Output = 0x2460400000000000607f9f80000000003fff8000400000
diff --git a/src/tests/data/bn/powmod.vec b/src/tests/data/bn/powmod.vec
index a84659768..024a2d36b 100644
--- a/src/tests/data/bn/powmod.vec
+++ b/src/tests/data/bn/powmod.vec
@@ -1,4 +1,5 @@
 [ModExp]
+
 Base = 0x1
 Exponent = 0x0
 Modulus = 0x2
@@ -193,3 +194,25 @@ Base = 0xF9FFFFF000000FFFFFFFFFFFFFFFC0000000
 Exponent = 0x83FFFF000000000000000003FFFFFFFFFFFF
 Modulus = 0xFFE007FFF9F83FFFFF8F000FFFFFFFFFFFFF
 Output = 0xA917797602DADCC854BD67D27E86BB1D6575
+
+# OSS-Fuzz #287 followed by some variations
+Base = 1024
+Exponent = 0x1000000000000000000000000000000000000000000000000000000000030400000000000004000
+Modulus = 0x40000000000000000000000000000000000000000000000000000000000c100000000000000ffff
+Output = 32
+
+Base = 1024
+Exponent = 0x1000000000000000000000000000000000000000000000000000000000030400000000000004001
+Modulus = 0x40000000000000000000000000000000000000000000000000000000000c100000000000000ffff
+Output = 32768
+
+Base = 1024
+Exponent = 0x1000000000000000000000000000000000000000000000000000000000030400000000000004001
+Modulus = 0x40000000000000000000000000000000000000000000000000000000000c100000000000000fffd
+Output = 0x3a130cf2e4c28710c6661c071e38d642150f28479017f313d2a855564458f9a958753fb286d6b6e
+
+Base = 1024
+Exponent = 0x20000000000000000000000000000000000000000000000000000000000608000000000000080
+Modulus  = 0x40000000000000000000000000000000000000000000000000000000000c100000000000000fffd
+Output = 0x9943fa648c48cb4cd01756bed11e3382aca74d84fb0bf8cf8d56cd4524f80538d4888cbd23b8e2
+
diff --git a/src/tests/test_bigint.cpp b/src/tests/test_bigint.cpp
index 27de5cfcb..ad837cc91 100644
--- a/src/tests/test_bigint.cpp
+++ b/src/tests/test_bigint.cpp
@@ -404,12 +404,37 @@ class BigInt_Powmod_Test : public Text_Based_Test
          {
          Test::Result result("BigInt Powmod");
 
-         const BigInt value = get_req_bn(vars, "Base");
+         const BigInt base = get_req_bn(vars, "Base");
          const BigInt exponent = get_req_bn(vars, "Exponent");
          const BigInt modulus = get_req_bn(vars, "Modulus");
-         const BigInt output = get_req_bn(vars, "Output");
+         const BigInt expected = get_req_bn(vars, "Output");
+
+         result.test_eq("power_mod", Botan::power_mod(base, exponent, modulus), expected);
+
+         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);
 
-         result.test_eq("power_mod", Botan::power_mod(value, exponent, modulus), output);
+         // 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;
          }