Repository: randombit/botan
Pull request #2287 (2020-03-01 12:35:46)

Remove use of Binary Extended Euclidean Algorithm for inversion
Instead use two specialized algorithms, one for odd modulus and the other for power of 2 modulus, then combine the results using CRT.

codecov-io: # [Codecov](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`master@f34cfff`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `87.66%`.

[![Impacted file tree graph](https://codecov.io/gh/randombit/botan/pull/2287/graphs/tree.svg?width=650&token=VR6iYmSPH8&height=150&src=pr)](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##             master    #2287   +/-   ##
=========================================
  Coverage          ?   92.59%           
=========================================
  Files             ?      557           
  Lines             ?    60616           
  Branches          ?     6427           
=========================================
  Hits              ?    56126           
  Misses            ?     4490           
  Partials          ?        0
```


| [Impacted Files](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [src/lib/math/numbertheory/numthry.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2xpYi9tYXRoL251bWJlcnRoZW9yeS9udW10aHJ5LmNwcA==) | `92% <ø> (ø)` | |
| [src/lib/math/numbertheory/monty.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2xpYi9tYXRoL251bWJlcnRoZW9yeS9tb250eS5jcHA=) | `65.51% <0%> (ø)` | |
| [src/lib/pubkey/rsa/rsa.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2xpYi9wdWJrZXkvcnNhL3JzYS5jcHA=) | `98.02% <100%> (ø)` | |
| [src/tests/test\_bigint.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL3Rlc3RzL3Rlc3RfYmlnaW50LmNwcA==) | `98.43% <100%> (ø)` | |
| [src/lib/math/numbertheory/make\_prm.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2xpYi9tYXRoL251bWJlcnRoZW9yeS9tYWtlX3BybS5jcHA=) | `91.39% <100%> (ø)` | |
| [src/cli/speed.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2NsaS9zcGVlZC5jcHA=) | `94.6% <21.42%> (ø)` | |
| [src/lib/math/numbertheory/mod\_inv.cpp](https://codecov.io/gh/randombit/botan/pull/2287/diff?src=pr&el=tree#diff-c3JjL2xpYi9tYXRoL251bWJlcnRoZW9yeS9tb2RfaW52LmNwcA==) | `94.81% <94.81%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=footer). Last update [f34cfff...6748ce2](https://codecov.io/gh/randombit/botan/pull/2287?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/src/cli/speed.cpp b/src/cli/speed.cpp
index b33dfc157..e876b7205 100644
--- a/src/cli/speed.cpp
+++ b/src/cli/speed.cpp
@@ -559,6 +559,10 @@ class Speed final : public Command
                {
                bench_rsa(provider, msec);
                }
+            else if(algo == "RSA_keygen")
+               {
+               bench_rsa_keygen(provider, msec);
+               }
 #endif
 #if defined(BOTAN_HAS_ECDSA)
             else if(algo == "ECDSA")
@@ -1535,28 +1539,22 @@ class Speed final : public Command
 
             const std::string bit_str = std::to_string(bits);
 
-            std::unique_ptr<Timer> invmod_timer = make_timer("binext-" + bit_str);
             std::unique_ptr<Timer> monty_timer = make_timer("monty-" + bit_str);
             std::unique_ptr<Timer> ct_invmod_timer = make_timer("ct-" + bit_str);
 
-            while(invmod_timer->under(runtime))
+            while(monty_timer->under(runtime))
                {
                const Botan::BigInt x(rng(), p.bits() - 1);
 
-               const Botan::BigInt x_inv1 = invmod_timer->run(
-                  [&] { return Botan::inverse_euclid(x, p); });
-
-               const Botan::BigInt x_inv2 = monty_timer->run(
+               const Botan::BigInt x_inv1 = 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_inv2 = ct_invmod_timer->run(
+                  [&] { return Botan::inverse_mod(x, p); });
 
                BOTAN_ASSERT_EQUAL(x_inv1, x_inv2, "Same result");
-               BOTAN_ASSERT_EQUAL(x_inv1, x_inv3, "Same result");
                }
 
-            record_result(invmod_timer);
             record_result(monty_timer);
             record_result(ct_invmod_timer);
             }
@@ -1845,6 +1843,25 @@ class Speed final : public Command
 #endif
 
 #if defined(BOTAN_HAS_RSA)
+      void bench_rsa_keygen(const std::string& provider,
+                            std::chrono::milliseconds msec)
+         {
+         for(size_t keylen : { 1024, 2048, 3072, 4096 })
+            {
+            const std::string nm = "RSA-" + std::to_string(keylen);
+            std::unique_ptr<Timer> keygen_timer = make_timer(nm, provider, "keygen");
+
+            while(keygen_timer->under(msec))
+               {
+               std::unique_ptr<Botan::Private_Key> key(keygen_timer->run([&] {
+                  return Botan::create_private_key("RSA", rng(), std::to_string(keylen));
+                  }));
+               }
+
+            record_result(keygen_timer);
+            }
+         }
+
       void bench_rsa(const std::string& provider,
                      std::chrono::milliseconds msec)
          {
diff --git a/src/fuzzer/invert.cpp b/src/fuzzer/invert.cpp
index 81923b07e..bd426a17e 100644
--- a/src/fuzzer/invert.cpp
+++ b/src/fuzzer/invert.cpp
@@ -6,6 +6,52 @@
 #include "fuzzers.h"
 #include <botan/numthry.h>
 
+namespace {
+
+Botan::BigInt ref_inverse_mod(const Botan::BigInt& n, const Botan::BigInt& mod)
+   {
+   if(n == 0 || mod < 2)
+      return 0;
+   if(n.is_even() && mod.is_even())
+      return 0;
+   Botan::BigInt u = mod, v = n;
+   Botan::BigInt A = 1, B = 0, C = 0, D = 1;
+
+   while(u.is_nonzero())
+      {
+      const size_t u_zero_bits = Botan::low_zero_bits(u);
+      u >>= u_zero_bits;
+      for(size_t i = 0; i != u_zero_bits; ++i)
+         {
+         if(A.is_odd() || B.is_odd())
+            { A += n; B -= mod; }
+         A >>= 1; B >>= 1;
+         }
+
+      const size_t v_zero_bits = Botan::low_zero_bits(v);
+      v >>= v_zero_bits;
+      for(size_t i = 0; i != v_zero_bits; ++i)
+         {
+         if(C.is_odd() || D.is_odd())
+            { C += n; D -= mod; }
+         C >>= 1; D >>= 1;
+         }
+
+      if(u >= v) { u -= v; A -= C; B -= D; }
+      else       { v -= u; C -= A; D -= B; }
+      }
+
+   if(v != 1)
+      return 0; // no modular inverse
+
+   while(D.is_negative()) D += mod;
+   while(D >= mod) D -= mod;
+
+   return D;
+   }
+
+}
+
 void fuzz(const uint8_t in[], size_t len)
    {
    if(len % 2 == 1 || len > 2*4096/8)
@@ -16,24 +62,21 @@ void fuzz(const uint8_t in[], size_t len)
    const Botan::BigInt x = Botan::BigInt::decode(in, part_len);
    Botan::BigInt mod = Botan::BigInt::decode(in + part_len, part_len);
 
-   mod.set_bit(0);
-
-   if(mod < 3 || x >= mod)
+   if(mod < 2)
       return;
 
-   const Botan::BigInt ref = Botan::inverse_euclid(x, mod);
-   const Botan::BigInt ct = Botan::ct_inverse_mod_odd_modulus(x, mod);
-   //Botan::BigInt mon = Botan::normalized_montgomery_inverse(x, mod);
+   const Botan::BigInt ref = ref_inverse_mod(x, mod);
+   const Botan::BigInt lib = Botan::inverse_mod(x, mod);
 
-   if(ref != ct)
+   if(ref != lib)
       {
       FUZZER_WRITE_AND_CRASH("X = " << x << "\n"
                              << "P = " << mod << "\n"
                              << "GCD = " << gcd(x, mod) << "\n"
                              << "Ref = " << ref << "\n"
-                             << "CT  = " << ct << "\n"
-                             << "RefCheck = " << (ref*ref)%mod << "\n"
-                             << "CTCheck  = " << (ct*ct)%mod << "\n");
+                             << "Lib  = " << lib << "\n"
+                             << "RefCheck = " << (x*ref)%mod << "\n"
+                             << "LibCheck  = " << (x*lib)%mod << "\n");
       }
    }
 
diff --git a/src/lib/math/numbertheory/make_prm.cpp b/src/lib/math/numbertheory/make_prm.cpp
index 00055d594..6c2c20718 100644
--- a/src/lib/math/numbertheory/make_prm.cpp
+++ b/src/lib/math/numbertheory/make_prm.cpp
@@ -171,21 +171,11 @@ BigInt random_prime(RandomNumberGenerator& rng,
          if(coprime > 1)
             {
             /*
-            * Check if gcd(p - 1, coprime) != 1 by attempting to compute a
-            * modular inverse. We are assured coprime is odd (since if it was
-            * even, it would always have a common factor with p - 1), so
-            * take off the factors of 2 from (p-1) then compute the inverse
-            * of coprime modulo that integer.
+            * Check if p - 1 and coprime are relatively prime, using gcd.
+            * The gcd computation is const-time
             */
-            BigInt p1 = p - 1;
-            p1 >>= low_zero_bits(p1);
-            if(ct_inverse_mod_odd_modulus(coprime, p1).is_zero())
-               {
-               BOTAN_DEBUG_ASSERT(gcd(p - 1, coprime) > 1);
+            if(gcd(p - 1, coprime) > 1)
                continue;
-               }
-
-            BOTAN_DEBUG_ASSERT(gcd(p - 1, coprime) == 1);
             }
 
          if(p.bits() > bits)
@@ -256,26 +246,10 @@ BigInt generate_rsa_prime(RandomNumberGenerator& keygen_rng,
             continue;
 
          /*
-         * Check if p - 1 and coprime are relatively prime by computing the inverse.
-         *
-         * We avoid gcd here because that algorithm is not constant time.
-         * Modular inverse is (for odd modulus anyway).
-         *
-         * We earlier verified that coprime argument is odd. Thus the factors of 2
-         * in (p - 1) cannot possibly be factors in coprime, so remove them from p - 1.
-         * Using an odd modulus allows the const time algorithm to be used.
-         *
-         * This assumes coprime < p - 1 which is always true for RSA.
+         * Check if p - 1 and coprime are relatively prime.
          */
-         BigInt p1 = p - 1;
-         p1 >>= low_zero_bits(p1);
-         if(ct_inverse_mod_odd_modulus(coprime, p1).is_zero())
-            {
-            BOTAN_DEBUG_ASSERT(gcd(p - 1, coprime) > 1);
+         if(gcd(p - 1, coprime) > 1)
             continue;
-            }
-
-         BOTAN_DEBUG_ASSERT(gcd(p - 1, coprime) == 1);
 
          if(p.bits() > bits)
             break;
diff --git a/src/lib/math/numbertheory/mod_inv.cpp b/src/lib/math/numbertheory/mod_inv.cpp
new file mode 100644
index 000000000..1827e8d53
--- /dev/null
+++ b/src/lib/math/numbertheory/mod_inv.cpp
@@ -0,0 +1,338 @@
+/*
+* (C) 1999-2011,2016,2018,2019,2020 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include <botan/numthry.h>
+#include <botan/divide.h>
+#include <botan/internal/ct_utils.h>
+#include <botan/internal/mp_core.h>
+
+namespace Botan {
+
+/*
+Sets result to a^-1 * 2^k mod a
+with n <= k <= 2n
+Returns k
+
+"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
+https://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
+*/
+size_t almost_montgomery_inverse(BigInt& result,
+                                 const BigInt& a,
+                                 const BigInt& p)
+   {
+   size_t k = 0;
+
+   BigInt u = p, v = a, r = 0, s = 1;
+
+   while(v > 0)
+      {
+      if(u.is_even())
+         {
+         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;
+         }
+
+      ++k;
+      }
+
+   if(r >= p)
+      {
+      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;
+   }
+
+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");
+
+   /*
+   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
+      https://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)
+         {
+         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);
+
+      //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);
+      }
+
+   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());
+
+   BOTAN_ASSERT(a.is_zero(), "A is zero");
+
+   if(b != 1)
+      return 0;
+
+   return v;
+   }
+
+BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
+   {
+   /*
+   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
+   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
+   */
+
+   if(a1.is_even())
+      return 0;
+
+   BigInt a = a1;
+   a.mask_bits(k);
+
+   BigInt b = 1;
+   BigInt X = 0;
+
+   for(size_t i = 0; i != k; ++i)
+      {
+      const bool b0 = b.get_bit(0);
+
+      X.conditionally_set_bit(i, b0);
+      b -= a * b0;
+      b /= 2;
+      }
+
+   return X;
+   }
+
+}
+
+BigInt inverse_mod(const BigInt& n, const BigInt& mod)
+   {
+   if(mod.is_zero())
+      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())
+      return 0;
+
+   if(mod.is_odd())
+      {
+      /*
+      Fastpath for common case. This leaks information if n > mod
+      but we don't guarantee const time behavior in that case.
+      */
+      if(n < mod)
+         return inverse_mod_odd_modulus(n, mod);
+      else
+         return inverse_mod_odd_modulus(ct_modulo(n, mod), mod);
+      }
+
+   const size_t mod_lz = low_zero_bits(mod);
+   BOTAN_ASSERT_NOMSG(mod_lz > 0);
+   const size_t mod_bits = mod.bits();
+   BOTAN_ASSERT_NOMSG(mod_bits > mod_lz);
+
+   if(mod_lz == mod_bits - 1)
+      {
+      // In this case we are performing an inversion modulo 2^k
+      return inverse_mod_pow2(n, mod_lz);
+      }
+
+   /*
+   * In this case we are performing an inversion modulo 2^k*o for
+   * some k > 1 and some odd (not necessarily prime) integer.
+   * Compute the inversions modulo 2^k and modulo o, then combine them
+   * using CRT, which is possible because 2^k and o are relatively prime.
+   */
+
+   const BigInt o = mod >> mod_lz;
+   const BigInt inv_o = inverse_mod_odd_modulus(ct_modulo(n, o), o);
+   const BigInt inv_2k = inverse_mod_pow2(n, mod_lz);
+
+   // No modular inverse in this case:
+   if(inv_o == 0 || inv_2k == 0)
+      return 0;
+
+   // Compute the CRT parameter
+   const BigInt c = inverse_mod_pow2(o, mod_lz);
+
+   // 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;
+
+   // Return result inv_o + h * o
+   h *= o;
+   h += inv_o;
+   return h;
+   }
+
+word monty_inverse(word a)
+   {
+   if(a % 2 == 0)
+      throw Invalid_Argument("monty_inverse only valid for odd integers");
+
+   /*
+   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
+   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
+   */
+
+   word b = 1;
+   word r = 0;
+
+   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
+      {
+      const word bi = b % 2;
+      r >>= 1;
+      r += bi << (BOTAN_MP_WORD_BITS - 1);
+
+      b -= a * bi;
+      b >>= 1;
+      }
+
+   // Now invert in addition space
+   r = (MP_WORD_MAX - r) + 1;
+
+   return r;
+   }
+
+}
diff --git a/src/lib/math/numbertheory/monty.cpp b/src/lib/math/numbertheory/monty.cpp
index b05bdd30b..ef75a587a 100644
--- a/src/lib/math/numbertheory/monty.cpp
+++ b/src/lib/math/numbertheory/monty.cpp
@@ -49,7 +49,8 @@ Montgomery_Params::Montgomery_Params(const BigInt& p)
 
 BigInt Montgomery_Params::inv_mod_p(const BigInt& x) const
    {
-   return ct_inverse_mod_odd_modulus(x, p());
+   // TODO use Montgomery inverse here?
+   return inverse_mod(x, p());
    }
 
 BigInt Montgomery_Params::redc(const BigInt& x, secure_vector<word>& ws) const
diff --git a/src/lib/math/numbertheory/numthry.cpp b/src/lib/math/numbertheory/numthry.cpp
index 90f1279b6..9fbaf8429 100644
--- a/src/lib/math/numbertheory/numthry.cpp
+++ b/src/lib/math/numbertheory/numthry.cpp
@@ -10,9 +10,7 @@
 #include <botan/monty.h>
 #include <botan/divide.h>
 #include <botan/rng.h>
-#include <botan/internal/bit_ops.h>
 #include <botan/internal/mp_core.h>
-#include <botan/internal/ct_utils.h>
 #include <botan/internal/monty_exp.h>
 #include <botan/internal/primality.h>
 #include <algorithm>
@@ -122,329 +120,6 @@ BigInt lcm(const BigInt& a, const BigInt& b)
    return ct_divide(a * b, gcd(a, b));
    }
 
-/*
-Sets result to a^-1 * 2^k mod a
-with n <= k <= 2n
-Returns k
-
-"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
-https://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
-*/
-size_t almost_montgomery_inverse(BigInt& result,
-                                 const BigInt& a,
-                                 const BigInt& p)
-   {
-   size_t k = 0;
-
-   BigInt u = p, v = a, r = 0, s = 1;
-
-   while(v > 0)
-      {
-      if(u.is_even())
-         {
-         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;
-         }
-
-      ++k;
-      }
-
-   if(r >= p)
-      {
-      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");
-   if(n >= mod)
-      throw Invalid_Argument("ct_inverse_mod_odd_modulus n >= mod not supported");
-
-   /*
-   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
-      https://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)
-         {
-         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);
-
-      //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);
-      }
-
-   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());
-
-   BOTAN_ASSERT(a.is_zero(), "A is zero");
-
-   if(b != 1)
-      return 0;
-
-   return v;
-   }
-
-/*
-* Find the Modular Inverse
-*/
-BigInt inverse_mod(const BigInt& n, const BigInt& mod)
-   {
-   if(mod.is_zero())
-      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(mod.is_odd() && n < mod)
-      return ct_inverse_mod_odd_modulus(n, mod);
-
-   return inverse_euclid(n, mod);
-   }
-
-BigInt inverse_euclid(const BigInt& n, const BigInt& mod)
-   {
-   if(mod.is_zero())
-      throw BigInt::DivideByZero();
-   if(mod.is_negative() || n.is_negative())
-      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
-
-   if(n.is_zero() || (n.is_even() && mod.is_even()))
-      return 0; // fast fail checks
-
-   BigInt u = mod, v = n;
-   BigInt A = 1, B = 0, C = 0, D = 1;
-   BigInt T0, T1, T2;
-
-   while(u.is_nonzero())
-      {
-      const size_t u_zero_bits = low_zero_bits(u);
-      u >>= u_zero_bits;
-
-      const size_t v_zero_bits = low_zero_bits(v);
-      v >>= v_zero_bits;
-
-      const bool u_gte_v = (u >= v);
-
-      for(size_t i = 0; i != u_zero_bits; ++i)
-         {
-         const bool needs_adjust = A.is_odd() || B.is_odd();
-
-         T0 = A + n;
-         T1 = B - mod;
-
-         A.ct_cond_assign(needs_adjust, T0);
-         B.ct_cond_assign(needs_adjust, T1);
-
-         A >>= 1;
-         B >>= 1;
-         }
-
-      for(size_t i = 0; i != v_zero_bits; ++i)
-         {
-         const bool needs_adjust = C.is_odd() || D.is_odd();
-         T0 = C + n;
-         T1 = D - mod;
-
-         C.ct_cond_assign(needs_adjust, T0);
-         D.ct_cond_assign(needs_adjust, T1);
-
-         C >>= 1;
-         D >>= 1;
-         }
-
-      T0 = u - v;
-      T1 = A - C;
-      T2 = B - D;
-
-      T0.cond_flip_sign(!u_gte_v);
-      T1.cond_flip_sign(!u_gte_v);
-      T2.cond_flip_sign(!u_gte_v);
-
-      u.ct_cond_assign(u_gte_v, T0);
-      A.ct_cond_assign(u_gte_v, T1);
-      B.ct_cond_assign(u_gte_v, T2);
-
-      v.ct_cond_assign(!u_gte_v, T0);
-      C.ct_cond_assign(!u_gte_v, T1);
-      D.ct_cond_assign(!u_gte_v, T2);
-      }
-
-   if(v != 1)
-      return 0; // no modular inverse
-
-   while(D.is_negative())
-      D += mod;
-   while(D >= mod)
-      D -= mod;
-
-   return D;
-   }
-
-word monty_inverse(word a)
-   {
-   if(a % 2 == 0)
-      throw Invalid_Argument("monty_inverse only valid for odd integers");
-
-   /*
-   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
-   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
-   */
-
-   word b = 1;
-   word r = 0;
-
-   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
-      {
-      const word bi = b % 2;
-      r >>= 1;
-      r += bi << (BOTAN_MP_WORD_BITS - 1);
-
-      b -= a * bi;
-      b >>= 1;
-      }
-
-   // Now invert in addition space
-   r = (MP_WORD_MAX - r) + 1;
-
-   return r;
-   }
-
 /*
 * Modular Exponentiation
 */
diff --git a/src/lib/math/numbertheory/numthry.h b/src/lib/math/numbertheory/numthry.h
index aab62a838..4266543cf 100644
--- a/src/lib/math/numbertheory/numthry.h
+++ b/src/lib/math/numbertheory/numthry.h
@@ -77,30 +77,40 @@ BigInt BOTAN_PUBLIC_API(2,0) lcm(const BigInt& x, const BigInt& y);
 BigInt BOTAN_PUBLIC_API(2,0) square(const BigInt& x);
 
 /**
-* Modular inversion
+* Modular inversion. This algorithm is const time with respect to 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.
+*
 * @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_PUBLIC_API(2,0) inverse_mod(const BigInt& x,
                                          const BigInt& modulus);
 
 /**
-* Modular inversion using extended binary Euclidian algorithm
+* Modular inversion
 * @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_PUBLIC_API(2,5) inverse_euclid(const BigInt& x,
-                                            const BigInt& modulus);
+inline BigInt BOTAN_DEPRECATED("Use inverse_mod")
+   inverse_euclid(const BigInt& x, const BigInt& modulus)
+   {
+   return inverse_mod(x, modulus);
+   }
 
 /**
 * Const time modular inversion
 * Requires the modulus be odd
 */
-BigInt BOTAN_PUBLIC_API(2,0) ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod);
+inline BigInt BOTAN_DEPRECATED("Use inverse_mod")
+   ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
+   {
+   return inverse_mod(n, mod);
+   }
 
 /**
 * Return a^-1 * 2^k mod b
diff --git a/src/lib/pubkey/rsa/rsa.cpp b/src/lib/pubkey/rsa/rsa.cpp
index bff1a1c15..bce6fae0f 100644
--- a/src/lib/pubkey/rsa/rsa.cpp
+++ b/src/lib/pubkey/rsa/rsa.cpp
@@ -298,11 +298,10 @@ RSA_PrivateKey::RSA_PrivateKey(RandomNumberGenerator& rng,
    const BigInt q_minus_1 = q - 1;
 
    const BigInt phi_n = lcm(p_minus_1, q_minus_1);
-   // FIXME: this uses binary ext gcd because phi_n is even
    d = inverse_mod(e, phi_n);
    d1 = ct_modulo(d, p_minus_1);
    d2 = ct_modulo(d, q_minus_1);
-   c = inverse_mod(q, p); // p odd, so uses const time algorithm
+   c = inverse_mod(q, p);
 
    RSA_PublicKey::init(std::move(n), std::move(e));
 


Test patch:
diff --git a/src/scripts/test_cli.py b/src/scripts/test_cli.py
index de6d92255..cafab6a09 100755
--- a/src/scripts/test_cli.py
+++ b/src/scripts/test_cli.py
@@ -1133,7 +1133,7 @@ def cli_speed_pk_tests(_tmp_dir):
 
     pk_algos = ["ECDSA", "ECDH", "SM2", "ECKCDSA", "ECGDSA", "GOST-34.10",
                 "DH", "DSA", "ElGamal", "Ed25519", "Curve25519", "NEWHOPE", "McEliece",
-                "RSA", "XMSS"]
+                "RSA", "RSA_keygen", "XMSS"]
 
     output = test_cli("speed", ["--msec=%d" % (msec)] + pk_algos, None).split('\n')
 
diff --git a/src/tests/data/bn/invmod.vec b/src/tests/data/bn/invmod.vec
index 99900b30b..9bcca8ee3 100644
--- a/src/tests/data/bn/invmod.vec
+++ b/src/tests/data/bn/invmod.vec
@@ -487,3 +487,44 @@ Input = 657900513264442578215729259525804042708991731028255426638688946278845827
 Modulus = 1460598742727828605836860681156746006166272494387112693708892744343548055514972863039476201659720247220048664250204314648520085411164461712526657028588699682983099362771
 Output = 370134468809235355809931021496418560257927761651782631431723120822241707286126622619231228275210678415132318220494260963802381448709723310690465171935975287188943190781
 
+# Moduli power of 2
+Input = 3
+Modulus = 0x10
+Output  = 11
+
+Input = 11
+Modulus = 0x100
+Output  = 163
+
+Input = 163
+Modulus = 0x1000
+Output  = 779
+
+Input = 2731
+Modulus = 0x10000
+Output  = 40963
+
+Input = 40963
+Modulus = 0x100000
+Output  = 461483
+
+Input = 461483
+Modulus = 0x1000000
+Output  = 12623875
+
+Input = 12623875
+Modulus = 0x10000000
+Output  = 50793131
+
+Input = 50793131
+Modulus = 0x100000000
+Output  = 2965413891
+
+Input = 2965413891
+Modulus = 0x1000000000
+Output  = 43000466091
+
+Input = 43000466091
+Modulus = 0x10000000000
+Output  = 140404367363
+
diff --git a/src/tests/test_bigint.cpp b/src/tests/test_bigint.cpp
index dbe308079..e612747ca 100644
--- a/src/tests/test_bigint.cpp
+++ b/src/tests/test_bigint.cpp
@@ -658,7 +658,7 @@ class BigInt_InvMod_Test final : public Text_Based_Test
          const Botan::BigInt mod = vars.get_req_bn("Modulus");
          const Botan::BigInt expected = vars.get_req_bn("Output");
 
-         const Botan::BigInt a_inv = Botan::inverse_euclid(a, mod);
+         const Botan::BigInt a_inv = Botan::inverse_mod(a, mod);
 
          result.test_eq("inverse_mod", a_inv, expected);
 
@@ -666,13 +666,12 @@ class BigInt_InvMod_Test final : public Text_Based_Test
             {
             result.test_eq("inverse ok", (a * a_inv) % mod, 1);
             }
-
-         if(mod.is_odd() && a < mod)
+         /*
+         else if((a % mod) > 0)
             {
-            result.test_eq("ct_inverse_odd_modulus",
-                           ct_inverse_mod_odd_modulus(a, mod),
-                           expected);
+            result.confirm("no inverse with gcd > 1", gcd(a, mod) > 1);
             }
+         */
 
          if(mod.is_odd() && a_inv != 0)
             {