Coverage Report

Created: 2022-06-23 06:44

/src/botan/src/lib/math/numbertheory/mod_inv.cpp
Line
Count
Source (jump to first uncovered line)
1
/*
2
* (C) 1999-2011,2016,2018,2019,2020 Jack Lloyd
3
*
4
* Botan is released under the Simplified BSD License (see license.txt)
5
*/
6
7
#include <botan/numthry.h>
8
#include <botan/internal/divide.h>
9
#include <botan/internal/ct_utils.h>
10
#include <botan/internal/mp_core.h>
11
#include <botan/internal/rounding.h>
12
13
namespace Botan {
14
15
namespace {
16
17
BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
18
87.8k
   {
19
   // Caller should assure these preconditions:
20
87.8k
   BOTAN_DEBUG_ASSERT(n.is_positive());
21
87.8k
   BOTAN_DEBUG_ASSERT(mod.is_positive());
22
87.8k
   BOTAN_DEBUG_ASSERT(n < mod);
23
87.8k
   BOTAN_DEBUG_ASSERT(mod >= 3 && mod.is_odd());
24
25
   /*
26
   This uses a modular inversion algorithm designed by Niels Möller
27
   and implemented in Nettle. The same algorithm was later also
28
   adapted to GMP in mpn_sec_invert.
29
30
   It can be easily implemented in a way that does not depend on
31
   secret branches or memory lookups, providing resistance against
32
   some forms of side channel attack.
33
34
   There is also a description of the algorithm in Appendix 5 of "Fast
35
   Software Polynomial Multiplication on ARM Processors using the NEON Engine"
36
   by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
37
   Dahab in LNCS 8182
38
      https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
39
40
   Thanks to Niels for creating the algorithm, explaining some things
41
   about it, and the reference to the paper.
42
   */
43
44
87.8k
   const size_t mod_words = mod.sig_words();
45
87.8k
   BOTAN_ASSERT(mod_words > 0, "Not empty");
46
47
87.8k
   secure_vector<word> tmp_mem(5*mod_words);
48
49
87.8k
   word* v_w = &tmp_mem[0];
50
87.8k
   word* u_w = &tmp_mem[1*mod_words];
51
87.8k
   word* b_w = &tmp_mem[2*mod_words];
52
87.8k
   word* a_w = &tmp_mem[3*mod_words];
53
87.8k
   word* mp1o2 = &tmp_mem[4*mod_words];
54
55
87.8k
   CT::poison(tmp_mem.data(), tmp_mem.size());
56
57
87.8k
   copy_mem(a_w, n.data(), std::min(n.size(), mod_words));
58
87.8k
   copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words));
59
87.8k
   u_w[0] = 1;
60
   // v_w = 0
61
62
   // compute (mod + 1) / 2 which [because mod is odd] is equal to
63
   // (mod / 2) + 1
64
87.8k
   copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words));
65
87.8k
   bigint_shr1(mp1o2, mod_words, 0, 1);
66
87.8k
   word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1);
67
87.8k
   BOTAN_ASSERT_NOMSG(carry == 0);
68
69
   // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
70
87.8k
   const size_t execs = 2 * mod.bits();
71
72
59.8M
   for(size_t i = 0; i != execs; ++i)
73
59.7M
      {
74
59.7M
      const word odd_a = a_w[0] & 1;
75
76
      //if(odd_a) a -= b
77
59.7M
      word underflow = bigint_cnd_sub(odd_a, a_w, b_w, mod_words);
78
79
      //if(underflow) { b -= a; a = abs(a); swap(u, v); }
80
59.7M
      bigint_cnd_add(underflow, b_w, a_w, mod_words);
81
59.7M
      bigint_cnd_abs(underflow, a_w, mod_words);
82
59.7M
      bigint_cnd_swap(underflow, u_w, v_w, mod_words);
83
84
      // a >>= 1
85
59.7M
      bigint_shr1(a_w, mod_words, 0, 1);
86
87
      //if(odd_a) u -= v;
88
59.7M
      word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words);
89
90
      // if(borrow) u += p
91
59.7M
      bigint_cnd_add(borrow, u_w, mod.data(), mod_words);
92
93
59.7M
      const word odd_u = u_w[0] & 1;
94
95
      // u >>= 1
96
59.7M
      bigint_shr1(u_w, mod_words, 0, 1);
97
98
      //if(odd_u) u += mp1o2;
99
59.7M
      bigint_cnd_add(odd_u, u_w, mp1o2, mod_words);
100
59.7M
      }
101
102
87.8k
   auto a_is_0 = CT::Mask<word>::set();
103
560k
   for(size_t i = 0; i != mod_words; ++i)
104
473k
      a_is_0 &= CT::Mask<word>::is_zero(a_w[i]);
105
106
87.8k
   auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1);
107
473k
   for(size_t i = 1; i != mod_words; ++i)
108
385k
      b_is_1 &= CT::Mask<word>::is_zero(b_w[i]);
109
110
87.8k
   BOTAN_ASSERT(a_is_0.is_set(), "A is zero");
111
112
   // if b != 1 then gcd(n,mod) > 1 and inverse does not exist
113
   // in which case zero out the result to indicate this
114
87.8k
   (~b_is_1).if_set_zero_out(v_w, mod_words);
115
116
   /*
117
   * We've placed the result in the lowest words of the temp buffer.
118
   * So just clear out the other values and then give that buffer to a
119
   * BigInt.
120
   */
121
87.8k
   clear_mem(&tmp_mem[mod_words], 4*mod_words);
122
123
87.8k
   CT::unpoison(tmp_mem.data(), tmp_mem.size());
124
125
87.8k
   BigInt r;
126
87.8k
   r.swap_reg(tmp_mem);
127
87.8k
   return r;
128
87.8k
   }
129
130
BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
131
2.34k
   {
132
   /*
133
   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
134
   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
135
   */
136
137
2.34k
   if(a1.is_even() || k == 0)
138
0
      return BigInt::zero();
139
2.34k
   if(k == 1)
140
8
      return BigInt::one();
141
142
2.33k
   BigInt a = a1;
143
2.33k
   a.mask_bits(k);
144
145
2.33k
   BigInt b = BigInt::one();
146
2.33k
   BigInt X = BigInt::zero();
147
2.33k
   BigInt newb;
148
149
2.33k
   const size_t a_words = a.sig_words();
150
151
2.33k
   X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS);
152
2.33k
   b.grow_to(a_words);
153
154
   /*
155
   Hide the exact value of k. k is anyway known to word length
156
   granularity because of the length of a, so no point in doing more
157
   than this.
158
   */
159
2.33k
   const size_t iter = round_up(k, BOTAN_MP_WORD_BITS);
160
161
858k
   for(size_t i = 0; i != iter; ++i)
162
856k
      {
163
856k
      const bool b0 = b.get_bit(0);
164
856k
      X.conditionally_set_bit(i, b0);
165
856k
      newb = b - a;
166
856k
      b.ct_cond_assign(b0, newb);
167
856k
      b >>= 1;
168
856k
      }
169
170
2.33k
   X.mask_bits(k);
171
2.33k
   X.const_time_unpoison();
172
2.33k
   return X;
173
2.34k
   }
174
175
}
176
177
BigInt inverse_mod(const BigInt& n, const BigInt& mod)
178
87.9k
   {
179
87.9k
   if(mod.is_zero())
180
0
      throw Invalid_Argument("inverse_mod modulus cannot be zero");
181
87.9k
   if(mod.is_negative() || n.is_negative())
182
0
      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
183
87.9k
   if(n.is_zero() || (n.is_even() && mod.is_even()))
184
22
      return BigInt::zero();
185
186
87.9k
   if(mod.is_odd())
187
86.5k
      {
188
      /*
189
      Fastpath for common case. This leaks if n is greater than mod or
190
      not, but we don't guarantee const time behavior in that case.
191
      */
192
86.5k
      if(n < mod)
193
86.2k
         return inverse_mod_odd_modulus(n, mod);
194
273
      else
195
273
         return inverse_mod_odd_modulus(ct_modulo(n, mod), mod);
196
86.5k
      }
197
198
   // If n is even and mod is even we already returned 0
199
   // If n is even and mod is odd we jumped directly to odd-modulus algo
200
1.40k
   BOTAN_DEBUG_ASSERT(n.is_odd());
201
202
1.40k
   const size_t mod_lz = low_zero_bits(mod);
203
1.40k
   BOTAN_ASSERT_NOMSG(mod_lz > 0);
204
1.40k
   const size_t mod_bits = mod.bits();
205
1.40k
   BOTAN_ASSERT_NOMSG(mod_bits > mod_lz);
206
207
1.40k
   if(mod_lz == mod_bits - 1)
208
70
      {
209
      // In this case we are performing an inversion modulo 2^k
210
70
      return inverse_mod_pow2(n, mod_lz);
211
70
      }
212
213
1.33k
   if(mod_lz == 1)
214
147
      {
215
      /*
216
      Inversion modulo 2*o is an easier special case of CRT
217
218
      This is exactly the main CRT flow below but taking advantage of
219
      the fact that any odd number ^-1 modulo 2 is 1. As a result both
220
      inv_2k and c can be taken to be 1, m2k is 2, and h is always
221
      either 0 or 1, and its value depends only on the low bit of inv_o.
222
223
      This is worth special casing because we generate RSA primes such
224
      that phi(n) is of this form. However this only works for keys
225
      that we generated in this way; pre-existing keys will typically
226
      fall back to the general algorithm below.
227
      */
228
229
147
      const BigInt o = mod >> 1;
230
147
      const BigInt n_redc = ct_modulo(n, o);
231
147
      const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o);
232
233
      // No modular inverse in this case:
234
147
      if(inv_o == 0)
235
23
         return BigInt::zero();
236
237
124
      BigInt h = inv_o;
238
124
      h.ct_cond_add(!inv_o.get_bit(0), o);
239
124
      return h;
240
147
      }
241
242
   /*
243
   * In this case we are performing an inversion modulo 2^k*o for
244
   * some k >= 2 and some odd (not necessarily prime) integer.
245
   * Compute the inversions modulo 2^k and modulo o, then combine them
246
   * using CRT, which is possible because 2^k and o are relatively prime.
247
   */
248
249
1.18k
   const BigInt o = mod >> mod_lz;
250
1.18k
   const BigInt n_redc = ct_modulo(n, o);
251
1.18k
   const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o);
252
1.18k
   const BigInt inv_2k = inverse_mod_pow2(n, mod_lz);
253
254
   // No modular inverse in this case:
255
1.18k
   if(inv_o == 0 || inv_2k == 0)
256
102
      return BigInt::zero();
257
258
1.08k
   const BigInt m2k = BigInt::power_of_2(mod_lz);
259
   // Compute the CRT parameter
260
1.08k
   const BigInt c = inverse_mod_pow2(o, mod_lz);
261
262
   // Compute h = c*(inv_2k-inv_o) mod 2^k
263
1.08k
   BigInt h = c * (inv_2k - inv_o);
264
1.08k
   const bool h_neg = h.is_negative();
265
1.08k
   h.set_sign(BigInt::Positive);
266
1.08k
   h.mask_bits(mod_lz);
267
1.08k
   const bool h_nonzero = h.is_nonzero();
268
1.08k
   h.ct_cond_assign(h_nonzero && h_neg, m2k - h);
269
270
   // Return result inv_o + h * o
271
1.08k
   h *= o;
272
1.08k
   h += inv_o;
273
1.08k
   return h;
274
1.18k
   }
275
276
}