For odd moduli use a input-independent modular inverse algorithm.

Also adds a (not const time) implementation of almost Montgomery reduction.

7 files changed, 311 insertions, 62 deletions

diff --git a/src/cli/speed.cpp b/src/cli/speed.cpp
index 45a0b4717..568aa09dd 100644
--- a/src/cli/speed.cpp
+++ b/src/cli/speed.cpp
@@ -107,7 +107,21 @@ class Timer
             const uint64_t now = get_clock();
 
             if(now > m_timer_start)
-               m_time_used += (now - m_timer_start);
+               {
+               uint64_t dur = now - m_timer_start;
+
+               m_time_used += dur;
+
+               if(m_event_count == 0)
+                  {
+                  m_min_time = m_max_time = dur;
+                  }
+               else
+                  {
+                  m_max_time = std::max(m_max_time, dur);
+                  m_min_time = std::min(m_min_time, dur);
+                  }
+               }
 
             m_timer_start = 0;
             ++m_event_count;
@@ -156,12 +170,17 @@ class Timer
       const std::string& get_name() const { return m_name; }
       const std::string& doing() const { return m_doing; }
 
+      uint64_t min_time() const { return m_min_time; }
+      uint64_t max_time() const { return m_max_time; }
+
       static std::string result_string_bps(const Timer& t);
       static std::string result_string_ops(const Timer& t);
    private:
       std::string m_name, m_doing;
       uint64_t m_time_used = 0, m_timer_start = 0;
       uint64_t m_event_count = 0, m_event_mult = 0;
+
+      uint64_t m_max_time = 0, m_min_time = 0;
    };
 
 std::string Timer::result_string_bps(const Timer& timer)
@@ -616,20 +635,46 @@ class Speed final : public Command
 
       void bench_inverse_mod(const std::chrono::milliseconds runtime)
          {
-         const Botan::BigInt p = random_prime(rng(), 1024);
+         Botan::BigInt p;
+         p.set_bit(521);
+         p--;
 
          Timer invmod_timer("inverse_mod");
+         Timer monty_timer("montgomery_inverse");
+         Timer ct_invmod_timer("ct_inverse_mod");
+         Timer powm_timer("exponentiation");
+
+         Botan::Fixed_Exponent_Power_Mod powm_p(p - 2, p);
 
          while(invmod_timer.under(runtime))
             {
-            const Botan::BigInt x(rng(), 1023);
+            const Botan::BigInt x(rng(), p.bits() - 1);
 
             const Botan::BigInt x_inv1 = invmod_timer.run([&]{
-               return Botan::inverse_mod(x, p);
+               return Botan::inverse_mod(x + p, p);
                });
+
+            const Botan::BigInt x_inv2 = monty_timer.run([&]{
+               return Botan::normalized_montgomery_inverse(x, p);
+               });
+
+            const Botan::BigInt x_inv3 = ct_invmod_timer.run([&]{
+               return Botan::ct_inverse_mod_odd_modulus(x, p);
+               });
+
+            const Botan::BigInt x_inv4 = powm_timer.run([&]{
+               return powm_p(x);
+               });
+
+            BOTAN_ASSERT_EQUAL(x_inv1, x_inv2, "Same result");
+            BOTAN_ASSERT_EQUAL(x_inv1, x_inv3, "Same result");
+            BOTAN_ASSERT_EQUAL(x_inv1, x_inv4, "Same result");
             }
 
          output() << Timer::result_string_ops(invmod_timer);
+         output() << Timer::result_string_ops(monty_timer);
+         output() << Timer::result_string_ops(ct_invmod_timer);
+         output() << Timer::result_string_ops(powm_timer);
          }
 
       void bench_random_prime(const std::chrono::milliseconds runtime)
diff --git a/src/lib/math/mp/mp_asm.cpp b/src/lib/math/mp/mp_asm.cpp
index 6d60d7e77..cfbb027d7 100644
--- a/src/lib/math/mp/mp_asm.cpp
+++ b/src/lib/math/mp/mp_asm.cpp
@@ -24,9 +24,6 @@ void bigint_cnd_swap(word cnd, word x[], word y[], size_t size)
    {
    const word mask = CT::expand_mask(cnd);
 
-   CT::poison(x, size);
-   CT::poison(y, size);
-
    for(size_t i = 0; i != size; ++i)
       {
       word a = x[i];
@@ -34,9 +31,6 @@ void bigint_cnd_swap(word cnd, word x[], word y[], size_t size)
       x[i] = CT::select(mask, b, a);
       y[i] = CT::select(mask, a, b);
       }
-
-   CT::unpoison(x, size);
-   CT::unpoison(y, size);
    }
 
 /*
@@ -47,9 +41,6 @@ word bigint_cnd_add(word cnd, word x[], const word y[], size_t size)
    {
    const word mask = CT::expand_mask(cnd);
 
-   CT::poison(x, size);
-   CT::poison(y, size);
-
    word carry = 0;
    for(size_t i = 0; i != size; ++i)
       {
@@ -61,10 +52,6 @@ word bigint_cnd_add(word cnd, word x[], const word y[], size_t size)
       x[i] = CT::select(mask, z, x[i]);
       }
 
-   CT::unpoison(x, size);
-   CT::unpoison(y, size);
-   CT::unpoison(carry);
-
    return carry & mask;
    }
 
@@ -76,9 +63,6 @@ word bigint_cnd_sub(word cnd, word x[], const word y[], size_t size)
    {
    const word mask = CT::expand_mask(cnd);
 
-   CT::poison(x, size);
-   CT::poison(y, size);
-
    word carry = 0;
    for(size_t i = 0; i != size; ++i)
       {
@@ -86,13 +70,21 @@ word bigint_cnd_sub(word cnd, word x[], const word y[], size_t size)
       x[i] = CT::select(mask, z, x[i]);
       }
 
-   CT::unpoison(x, size);
-   CT::unpoison(y, size);
-   CT::unpoison(carry);
-
    return carry & mask;
    }
 
+void bigint_cnd_abs(word cnd, word x[], size_t size)
+   {
+   const word mask = CT::expand_mask(cnd);
+
+   word carry = mask & 1;
+   for(size_t i = 0; i != size; ++i)
+      {
+      const word z = word_add(~x[i], 0, &carry);
+      x[i] = CT::select(mask, z, x[i]);
+      }
+   }
+
 /*
 * Two Operand Addition, No Carry
 */
diff --git a/src/lib/math/mp/mp_core.h b/src/lib/math/mp/mp_core.h
index 86bc920cf..73f13742c 100644
--- a/src/lib/math/mp/mp_core.h
+++ b/src/lib/math/mp/mp_core.h
@@ -40,6 +40,14 @@ word bigint_cnd_add(word cnd, word x[], const word y[], size_t size);
 BOTAN_DLL
 word bigint_cnd_sub(word cnd, word x[], const word y[], size_t size);
 
+/*
+* 2s complement absolute value
+* If cond > 0 sets x to ~x + 1
+* Runs in constant time
+*/
+BOTAN_DLL
+void bigint_cnd_abs(word cnd, word x[], size_t size);
+
 /**
 * Two operand addition
 * @param x the first operand (and output)
diff --git a/src/lib/math/numbertheory/numthry.cpp b/src/lib/math/numbertheory/numthry.cpp
index 982e9264a..ae2d33524 100644
--- a/src/lib/math/numbertheory/numthry.cpp
+++ b/src/lib/math/numbertheory/numthry.cpp
@@ -1,6 +1,6 @@
 /*
 * Number Theory Functions
-* (C) 1999-2011 Jack Lloyd
+* (C) 1999-2011,2016 Jack Lloyd
 *
 * Botan is released under the Simplified BSD License (see license.txt)
 */
@@ -9,6 +9,7 @@
 #include <botan/reducer.h>
 #include <botan/internal/bit_ops.h>
 #include <botan/internal/mp_core.h>
+#include <botan/internal/ct_utils.h>
 #include <algorithm>
 
 namespace Botan {
@@ -75,52 +76,198 @@ BigInt lcm(const BigInt& a, const BigInt& b)
    }
 
 /*
-* If the modulus is odd, then we can avoid computing A and C. This is
-* a critical path algorithm in some instances and an odd modulus is
-* the common case for crypto, so worth special casing. See note 14.64
-* in Handbook of Applied Cryptography for more details.
+Sets result to a^-1 * 2^k mod a
+with n <= k <= 2n
+Returns k
+
+"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
+http://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
 */
-BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
+size_t almost_montgomery_inverse(BigInt& result,
+                                 const BigInt& a,
+                                 const BigInt& p)
    {
-   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");
+   size_t k = 0;
 
-   BigInt u = mod, v = n;
-   BigInt B = 0, D = 1;
+   BigInt u = p, v = a, r = 0, s = 1;
 
-   while(u.is_nonzero())
+   while(v > 0)
       {
-      const size_t u_zero_bits = low_zero_bits(u);
-      u >>= u_zero_bits;
-      for(size_t i = 0; i != u_zero_bits; ++i)
+      if(u.is_even())
          {
-         if(B.is_odd())
-            { B -= mod; }
-         B >>= 1;
+         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;
          }
 
-      const size_t v_zero_bits = low_zero_bits(v);
-      v >>= v_zero_bits;
-      for(size_t i = 0; i != v_zero_bits; ++i)
+      ++k;
+      }
+
+   if(r >= p)
+      {
+      r = 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;
+   }
+
+BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
+   {
+   if(n.is_negative() || mod.is_negative())
+      throw Invalid_Argument("ct_inverse_mod_odd_modulus: arguments must be non-negative");
+   if(mod < 3 || mod.is_even())
+      throw Invalid_Argument("Bad modulus to ct_inverse_mod_odd_modulus");
+
+   /*
+   This uses a modular inversion algorithm designed by Niels Möller
+   and implemented in Nettle. The same algorithm was later also
+   adapted to GMP in mpn_sec_invert.
+
+   It can be easily implemented in a way that does not depend on
+   secret branches or memory lookups, providing resistance against
+   some forms of side channel attack.
+
+   There is also a description of the algorithm in Appendix 5 of "Fast
+   Software Polynomial Multiplication on ARM Processors using the NEON Engine"
+   by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
+   Dahab in LNCS 8182
+      http://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
+
+   Thanks to Niels for creating the algorithm, explaining some things
+   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;
+
+   a.grow_to(mod_words);
+   u.grow_to(mod_words);
+   v.grow_to(mod_words);
+   mp1o2.grow_to(mod_words);
+
+   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();
+
+   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());
+
+   // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
+   size_t bits = 2 * mod.bits();
+
+   while(bits--)
+      {
+      /*
+      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)
          {
-         if(D.is_odd())
-            { D -= mod; }
-         D >>= 1;
+         std::swap(u, v);
          }
 
-      if(u >= v) { u -= v; B -= D; }
-      else       { v -= u; D -= B; }
+      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);
+
+      //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);
+
+      // a >>= 1
+      bigint_shr1(a_w.data(), mod_words, 0, 1);
+
+      //if(odd_a) u -= v;
+      word borrow = bigint_cnd_sub(odd_a, u_w.data(), v_w.data(), mod_words);
+
+      // if(borrow) u += p
+      bigint_cnd_add(borrow, u_w.data(), mod.data(), mod_words);
+
+      const word odd_u = u_w[0] & 1;
+
+      // u >>= 1
+      bigint_shr1(u_w.data(), mod_words, 0, 1);
+
+      //if(odd_u) u += mp1o2;
+      bigint_cnd_add(odd_u, u_w.data(), mp1o2.data(), mod_words);
       }
 
-   if(v != 1)
-      return 0; // no modular inverse
+   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());
 
-   while(D.is_negative()) D += mod;
-   while(D >= mod) D -= mod;
+   BOTAN_ASSERT(a.is_zero(), "A is zero");
 
-   return D;
+   if(b != 1)
+      return 0;
+
+   return v;
    }
 
 /*
@@ -137,7 +284,7 @@ BigInt inverse_mod(const BigInt& n, const BigInt& mod)
       return 0; // fast fail checks
 
    if(mod.is_odd())
-      return inverse_mod_odd_modulus(n, mod);
+      return ct_inverse_mod_odd_modulus(n, mod);
 
    BigInt u = mod, v = n;
    BigInt A = 1, B = 0, C = 0, D = 1;
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index b9a082c6c..e1e6c65f6 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -70,15 +70,30 @@ BigInt BOTAN_DLL square(const BigInt& x);
 * @param x a positive integer
 * @param modulus a positive integer
 * @return y st (x*y) % modulus == 1 or 0 if no such value
+* Not const time
 */
 BigInt BOTAN_DLL inverse_mod(const BigInt& x,
                              const BigInt& modulus);
 
 /**
-* As above but requires modulus be odd
+* Const time modular inversion
+* Requires the modulus be odd
 */
-BigInt BOTAN_DLL inverse_mod_odd_modulus(const BigInt& x,
-                                         const BigInt& modulus);
+BigInt BOTAN_DLL ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod);
+
+/**
+* Return a^-1 * 2^k mod b
+* Returns k, between n and 2n
+* Not const time
+*/
+size_t BOTAN_DLL 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_DLL normalized_montgomery_inverse(const BigInt& a, const BigInt& b);
 
 
 /**
diff --git a/src/tests/test_bigint.cpp b/src/tests/test_bigint.cpp
index 671c76bff..1a615c374 100644
--- a/src/tests/test_bigint.cpp
+++ b/src/tests/test_bigint.cpp
@@ -448,8 +448,18 @@ class BigInt_InvMod_Test : public Text_Based_Test
 
          if(a_inv > 1)
             {
-            const Botan::BigInt z = (a * a_inv) % mod;
-            result.test_eq("inverse ok", z, 1);
+            result.test_eq("inverse ok", (a * a_inv) % mod, 1);
+            }
+
+         if(mod.is_odd())
+            {
+            result.test_eq("normalized_montgomery_inverse",
+                           normalized_montgomery_inverse(a, mod),
+                           expected);
+
+            result.test_eq("ct_inverse_odd_modulus",
+                           ct_inverse_mod_odd_modulus(a, mod),
+                           expected);
             }
 
          return result;
diff --git a/src/tests/test_mp.cpp b/src/tests/test_mp.cpp
index b52d93406..cbaf465a4 100644
--- a/src/tests/test_mp.cpp
+++ b/src/tests/test_mp.cpp
@@ -26,6 +26,7 @@ class MP_Unit_Tests : public Test
          results.push_back(test_cnd_swap());
          results.push_back(test_cnd_add());
          results.push_back(test_cnd_sub());
+         results.push_back(test_cnd_abs());
 
          return results;
          }
@@ -75,6 +76,37 @@ class MP_Unit_Tests : public Test
          return result;
          }
 
+      Result test_cnd_abs()
+         {
+         Result result("bigint_cnd_abs");
+
+         using namespace Botan;
+
+         word x1 = MP_WORD_MAX;
+         bigint_cnd_abs(1, &x1, 1);
+         result.test_int_eq(x1, 1, "Abs");
+
+         x1 = 0;
+         bigint_cnd_abs(1, &x1, 1);
+         result.test_int_eq(x1, 0, "Abs");
+
+         x1 = 1;
+         bigint_cnd_abs(1, &x1, 1);
+         result.test_int_eq(x1, MP_WORD_MAX, "Abs");
+
+         x1 = 1;
+         bigint_cnd_abs(0, &x1, 1);
+         result.test_int_eq(x1, 1, "No change");
+
+         word x2[2] = { MP_WORD_MAX, MP_WORD_MAX };
+
+         bigint_cnd_abs(1, x2, 2);
+         result.test_int_eq(x2[0], 1, "Abs");
+         result.test_int_eq(x2[1], 0, "Abs");
+
+         return result;
+         }
+
       Result test_cnd_swap()
          {
          Result result("bigint_cnd_swap");