diff options
Diffstat (limited to 'comm/third_party/botan/src/lib/math/numbertheory/mod_inv.cpp')
-rw-r--r-- | comm/third_party/botan/src/lib/math/numbertheory/mod_inv.cpp | 356 |
1 files changed, 356 insertions, 0 deletions
diff --git a/comm/third_party/botan/src/lib/math/numbertheory/mod_inv.cpp b/comm/third_party/botan/src/lib/math/numbertheory/mod_inv.cpp new file mode 100644 index 0000000000..ec3bb33f00 --- /dev/null +++ b/comm/third_party/botan/src/lib/math/numbertheory/mod_inv.cpp @@ -0,0 +1,356 @@ +/* +* (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> +#include <botan/internal/rounding.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) + { + // 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 + 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. + */ + + const size_t mod_words = mod.sig_words(); + BOTAN_ASSERT(mod_words > 0, "Not empty"); + + 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]; + + CT::poison(tmp_mem.data(), tmp_mem.size()); + + 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 + + // 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); + + // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n + const size_t execs = 2 * mod.bits(); + + for(size_t i = 0; i != execs; ++i) + { + const word odd_a = a_w[0] & 1; + + //if(odd_a) a -= b + 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, 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, mod_words, 0, 1); + + //if(odd_a) u -= v; + word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words); + + // if(borrow) u += p + bigint_cnd_add(borrow, u_w, mod.data(), mod_words); + + const word odd_u = u_w[0] & 1; + + // u >>= 1 + bigint_shr1(u_w, mod_words, 0, 1); + + //if(odd_u) u += mp1o2; + bigint_cnd_add(odd_u, u_w, mp1o2, mod_words); + } + + 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]); + + 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]); + + 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); + + /* + * We've placed the result in the lowest words of the temp buffer. + * So just clear out the other values and then give that buffer to a + * BigInt. + */ + 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) + { + /* + * 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; + BigInt newb; + + const size_t a_words = a.sig_words(); + + X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS); + b.grow_to(a_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); + newb = b - a; + b.ct_cond_assign(b0, newb); + b >>= 1; + } + + X.mask_bits(k); + X.const_time_unpoison(); + 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() || (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 n_redc = ct_modulo(n, o); + const BigInt inv_o = inverse_mod_odd_modulus(n_redc, 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; + + const BigInt m2k = BigInt::power_of_2(mod_lz); + // 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.set_sign(BigInt::Positive); + h.mask_bits(mod_lz); + const bool h_nonzero = h.is_nonzero(); + h.ct_cond_assign(h_nonzero && h_neg, m2k - h); + + // Return result inv_o + h * o + h *= o; + h += inv_o; + return h; + } + +// Deprecated forwarding functions: +BigInt inverse_euclid(const BigInt& x, const BigInt& modulus) + { + return inverse_mod(x, modulus); + } + +BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod) + { + return inverse_mod_odd_modulus(n, mod); + } + +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; + } + +} |