Optimize inverse_mod

About 25% faster

3 files changed, 95 insertions, 95 deletions

diff --git a/src/cli/speed.cpp b/src/cli/speed.cpp
index e876b7205..06c1e6ab9 100644
--- a/src/cli/speed.cpp
+++ b/src/cli/speed.cpp
@@ -1435,8 +1435,8 @@ class Speed final : public Command
 
             while(f_timer->under(runtime))
                {
-               e_timer->run([&]() { Botan::power_mod(group.get_g(), random_e, group.get_p()); });
-               f_timer->run([&]() { Botan::power_mod(group.get_g(), random_f, group.get_p()); });
+               e_timer->run([&]() { group.power_g_p(random_e); });
+               f_timer->run([&]() { group.power_g_p(random_f); });
                }
 
             record_result(e_timer);
@@ -1532,31 +1532,33 @@ class Speed final : public Command
 
       void bench_inverse_mod(const std::chrono::milliseconds runtime)
          {
-
-         for(size_t bits : { 256, 384, 512 })
+         for(size_t bits : { 256, 384, 512, 1024, 2048 })
             {
-            const Botan::BigInt p = Botan::random_prime(rng(), bits);
-
             const std::string bit_str = std::to_string(bits);
 
-            std::unique_ptr<Timer> monty_timer = make_timer("monty-" + bit_str);
-            std::unique_ptr<Timer> ct_invmod_timer = make_timer("ct-" + bit_str);
+            std::unique_ptr<Timer> timer = make_timer("inverse_mod-" + bit_str);
 
-            while(monty_timer->under(runtime))
+            while(timer->under(runtime))
                {
-               const Botan::BigInt x(rng(), p.bits() - 1);
+               const Botan::BigInt x(rng(), bits - 1);
+               Botan::BigInt mod(rng(), bits);
 
-               const Botan::BigInt x_inv1 = monty_timer->run(
-                  [&] { return Botan::normalized_montgomery_inverse(x, p); });
+               const Botan::BigInt x_inv = timer->run(
+                  [&] { return Botan::inverse_mod(x, mod); });
 
-               const Botan::BigInt x_inv2 = ct_invmod_timer->run(
-                  [&] { return Botan::inverse_mod(x, p); });
-
-               BOTAN_ASSERT_EQUAL(x_inv1, x_inv2, "Same result");
+               if(x_inv == 0)
+                  {
+                  const Botan::BigInt g = gcd(x, mod);
+                  BOTAN_ASSERT(g != 1, "Inversion only fails if gcd(x, mod) > 1");
+                  }
+               else
+                  {
+                  const Botan::BigInt check = (x_inv*x) % mod;
+                  BOTAN_ASSERT_EQUAL(check, 1, "Const time inversion correct");
+                  }
                }
 
-            record_result(monty_timer);
-            record_result(ct_invmod_timer);
+            record_result(timer);
             }
          }
 
@@ -1856,6 +1858,8 @@ class Speed final : public Command
                std::unique_ptr<Botan::Private_Key> key(keygen_timer->run([&] {
                   return Botan::create_private_key("RSA", rng(), std::to_string(keylen));
                   }));
+
+               BOTAN_ASSERT(key->check_key(rng(), true), "Key is ok");
                }
 
             record_result(keygen_timer);
diff --git a/src/lib/math/numbertheory/mod_inv.cpp b/src/lib/math/numbertheory/mod_inv.cpp
index 8a23707e0..8bc073d9c 100644
--- a/src/lib/math/numbertheory/mod_inv.cpp
+++ b/src/lib/math/numbertheory/mod_inv.cpp
@@ -8,6 +8,7 @@
 #include <botan/divide.h>
 #include <botan/internal/ct_utils.h>
 #include <botan/internal/mp_core.h>
+#include <botan/internal/rounding.h>
 
 namespace Botan {
 
@@ -101,12 +102,11 @@ namespace {
 
 BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
    {
-   if(n.is_negative() || mod.is_negative())
-      throw Invalid_Argument("inverse_mod_odd_modulus: arguments must be non-negative");
-   if(mod < 3 || mod.is_even())
-      throw Invalid_Argument("Bad modulus to inverse_mod_odd_modulus");
-   if(n >= mod)
-      throw Invalid_Argument("inverse_mod_odd_modulus n >= mod not supported");
+   // Caller should assure these preconditions:
+   BOTAN_DEBUG_ASSERT(n.is_positive());
+   BOTAN_DEBUG_ASSERT(mod.is_positive());
+   BOTAN_DEBUG_ASSERT(n < mod);
+   BOTAN_DEBUG_ASSERT(mod >= 3 && mod.is_odd());
 
    /*
    This uses a modular inversion algorithm designed by Niels Möller
@@ -127,98 +127,85 @@ BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
    about it, and the reference to the paper.
    */
 
-   // todo allow this to be pre-calculated and passed in as arg
-   BigInt mp1o2 = (mod + 1) >> 1;
-
    const size_t mod_words = mod.sig_words();
    BOTAN_ASSERT(mod_words > 0, "Not empty");
 
-   BigInt a = n;
-   BigInt b = mod;
-   BigInt u = 1, v = 0;
+   secure_vector<word> tmp_mem(5*mod_words);
+
+   word* v_w = &tmp_mem[0];
+   word* u_w = &tmp_mem[1*mod_words];
+   word* b_w = &tmp_mem[2*mod_words];
+   word* a_w = &tmp_mem[3*mod_words];
+   word* mp1o2 = &tmp_mem[4*mod_words];
 
-   a.grow_to(mod_words);
-   u.grow_to(mod_words);
-   v.grow_to(mod_words);
-   mp1o2.grow_to(mod_words);
+   copy_mem(a_w, n.data(), std::min(n.size(), mod_words));
+   copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words));
+   u_w[0] = 1;
+   // v_w = 0
 
-   secure_vector<word>& a_w = a.get_word_vector();
-   secure_vector<word>& b_w = b.get_word_vector();
-   secure_vector<word>& u_w = u.get_word_vector();
-   secure_vector<word>& v_w = v.get_word_vector();
+   // compute (mod + 1) / 2 which [because mod is odd] is equal to
+   // (mod / 2) + 1
+   copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words));
+   bigint_shr1(mp1o2, mod_words, 0, 1);
+   word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1);
+   BOTAN_ASSERT_NOMSG(carry == 0);
 
-   CT::poison(a_w.data(), a_w.size());
-   CT::poison(b_w.data(), b_w.size());
-   CT::poison(u_w.data(), u_w.size());
-   CT::poison(v_w.data(), v_w.size());
+   CT::poison(tmp_mem.data(), tmp_mem.size());
 
    // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
-   size_t bits = 2 * mod.bits();
+   size_t execs = 2 * mod.bits();
 
-   while(bits--)
+   while(execs--)
       {
-      /*
-      const word odd = a.is_odd();
-      a -= odd * b;
-      const word underflow = a.is_negative();
-      b += a * underflow;
-      a.set_sign(BigInt::Positive);
-
-      a >>= 1;
-
-      if(underflow)
-         {
-         std::swap(u, v);
-         }
-
-      u -= odd * v;
-      u += u.is_negative() * mod;
-
-      const word odd_u = u.is_odd();
-
-      u >>= 1;
-      u += mp1o2 * odd_u;
-      */
-
       const word odd_a = a_w[0] & 1;
 
       //if(odd_a) a -= b
-      word underflow = bigint_cnd_sub(odd_a, a_w.data(), b_w.data(), mod_words);
+      word underflow = bigint_cnd_sub(odd_a, a_w, b_w, mod_words);
 
       //if(underflow) { b -= a; a = abs(a); swap(u, v); }
-      bigint_cnd_add(underflow, b_w.data(), a_w.data(), mod_words);
-      bigint_cnd_abs(underflow, a_w.data(), mod_words);
-      bigint_cnd_swap(underflow, u_w.data(), v_w.data(), mod_words);
+      bigint_cnd_add(underflow, b_w, a_w, mod_words);
+      bigint_cnd_abs(underflow, a_w, mod_words);
+      bigint_cnd_swap(underflow, u_w, v_w, mod_words);
 
       // a >>= 1
-      bigint_shr1(a_w.data(), mod_words, 0, 1);
+      bigint_shr1(a_w, mod_words, 0, 1);
 
       //if(odd_a) u -= v;
-      word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
+      word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words);
 
       // if(borrow) u += p
-      bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
+      bigint_cnd_add(borrow, u_w, mod.data(), mod_words);
 
       const word odd_u = u_w[0] & 1;
 
       // u >>= 1
-      bigint_shr1(u_w.data(), mod_words, 0, 1);
+      bigint_shr1(u_w, mod_words, 0, 1);
 
       //if(odd_u) u += mp1o2;
-      bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
+      bigint_cnd_add(odd_u, u_w, mp1o2, mod_words);
       }
 
-   CT::unpoison(a_w.data(), a_w.size());
-   CT::unpoison(b_w.data(), b_w.size());
-   CT::unpoison(u_w.data(), u_w.size());
-   CT::unpoison(v_w.data(), v_w.size());
+   auto a_is_0 = CT::Mask<word>::set();
+   for(size_t i = 0; i != mod_words; ++i)
+      a_is_0 &= CT::Mask<word>::is_zero(a_w[i]);
 
-   BOTAN_ASSERT(a.is_zero(), "A is zero");
+   auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1);
+   for(size_t i = 1; i != mod_words; ++i)
+      b_is_1 &= CT::Mask<word>::is_zero(b_w[i]);
 
-   if(b != 1)
-      return 0;
+   BOTAN_ASSERT(a_is_0.is_set(), "A is zero");
+
+   // if b != 1 then gcd(n,mod) > 1 and inverse does not exist
+   // in which case zero out the result to indicate this
+   (~b_is_1).if_set_zero_out(v_w, mod_words);
 
-   return v;
+   clear_mem(&tmp_mem[mod_words], 4*mod_words);
+
+   CT::unpoison(tmp_mem.data(), tmp_mem.size());
+
+   BigInt r;
+   r.swap_reg(tmp_mem);
+   return r;
    }
 
 BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
@@ -237,15 +224,25 @@ BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
    BigInt b = 1;
    BigInt X = 0;
 
-   for(size_t i = 0; i != k; ++i)
+   const size_t a_words = 1 + a.sig_words();
+
+   /*
+   Hide the exact value of k. k is anyway known to word length
+   granularity because of the length of a, so no point in doing more
+   than this.
+   */
+   const size_t iter = round_up(k, BOTAN_MP_WORD_BITS);
+
+   for(size_t i = 0; i != iter; ++i)
       {
       const bool b0 = b.get_bit(0);
-
       X.conditionally_set_bit(i, b0);
-      b -= a * b0;
-      b /= 2;
+      b.grow_to(a_words);
+      bigint_cnd_sub(b0, b.mutable_data(), b.size(), a.data(), a.sig_words());
+      b >>= 1;
       }
 
+   X.mask_bits(k);
    return X;
    }
 
@@ -257,9 +254,7 @@ BigInt inverse_mod(const BigInt& n, const BigInt& mod)
       throw BigInt::DivideByZero();
    if(mod.is_negative() || n.is_negative())
       throw Invalid_Argument("inverse_mod: arguments must be non-negative");
-   if(n.is_zero())
-      return 0;
-   if(n.is_even() && mod.is_even())
+   if(n.is_zero() || (n.is_even() && mod.is_even()))
       return 0;
 
    if(mod.is_odd())
@@ -306,10 +301,10 @@ BigInt inverse_mod(const BigInt& n, const BigInt& mod)
    // Compute h = c*(inv_2k-inv_o) mod 2^k
    BigInt h = c * (inv_2k - inv_o);
    const bool h_neg = h.is_negative();
-   h.mask_bits(mod_lz);
    h.set_sign(BigInt::Positive);
-   if(h_neg && h > 0)
-      h = BigInt::power_of_2(mod_lz) - h;
+   h.mask_bits(mod_lz);
+   const bool h_nonzero = h.is_nonzero();
+   h.ct_cond_assign(h_nonzero && h_neg, BigInt::power_of_2(mod_lz) - h);
 
    // Return result inv_o + h * o
    h *= o;
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index b476db600..e818aa0bb 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -81,7 +81,8 @@ BigInt BOTAN_PUBLIC_API(2,0) square(const BigInt& x);
 * as long as x is less than modulus. It also avoids leaking
 * information about the modulus, except that it does leak which of 3
 * categories the modulus is in: an odd integer, a power of 2, or some
-* other even number.
+* other even number, and if the modulus is even, leaks the power of 2
+* which divides the modulus.
 *
 * @param x a positive integer
 * @param modulus a positive integer