Coverage Report

Created: 2020-05-23 13:54

/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/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
/*
16
Sets result to a^-1 * 2^k mod a
17
with n <= k <= 2n
18
Returns k
19
20
"The Montgomery Modular Inverse - Revisited" Çetin Koç, E. Savas
21
https://citeseerx.ist.psu.edu/viewdoc/citations?doi=10.1.1.75.8377
22
23
A const time implementation of this algorithm is described in
24
"Constant Time Modular Inversion" Joppe W. Bos
25
http://www.joppebos.com/files/CTInversion.pdf
26
*/
27
size_t almost_montgomery_inverse(BigInt& result,
28
                                 const BigInt& a,
29
                                 const BigInt& p)
30
0
   {
31
0
   size_t k = 0;
32
0
33
0
   BigInt u = p, v = a, r = 0, s = 1;
34
0
35
0
   while(v > 0)
36
0
      {
37
0
      if(u.is_even())
38
0
         {
39
0
         u >>= 1;
40
0
         s <<= 1;
41
0
         }
42
0
      else if(v.is_even())
43
0
         {
44
0
         v >>= 1;
45
0
         r <<= 1;
46
0
         }
47
0
      else if(u > v)
48
0
         {
49
0
         u -= v;
50
0
         u >>= 1;
51
0
         r += s;
52
0
         s <<= 1;
53
0
         }
54
0
      else
55
0
         {
56
0
         v -= u;
57
0
         v >>= 1;
58
0
         s += r;
59
0
         r <<= 1;
60
0
         }
61
0
62
0
      ++k;
63
0
      }
64
0
65
0
   if(r >= p)
66
0
      {
67
0
      r -= p;
68
0
      }
69
0
70
0
   result = p - r;
71
0
72
0
   return k;
73
0
   }
74
75
BigInt normalized_montgomery_inverse(const BigInt& a, const BigInt& p)
76
0
   {
77
0
   BigInt r;
78
0
   size_t k = almost_montgomery_inverse(r, a, p);
79
0
80
0
   for(size_t i = 0; i != k; ++i)
81
0
      {
82
0
      if(r.is_odd())
83
0
         r += p;
84
0
      r >>= 1;
85
0
      }
86
0
87
0
   return r;
88
0
   }
89
90
namespace {
91
92
BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
93
42.8k
   {
94
42.8k
   // Caller should assure these preconditions:
95
42.8k
   BOTAN_DEBUG_ASSERT(n.is_positive());
96
42.8k
   BOTAN_DEBUG_ASSERT(mod.is_positive());
97
42.8k
   BOTAN_DEBUG_ASSERT(n < mod);
98
42.8k
   BOTAN_DEBUG_ASSERT(mod >= 3 && mod.is_odd());
99
42.8k
100
42.8k
   /*
101
42.8k
   This uses a modular inversion algorithm designed by Niels Möller
102
42.8k
   and implemented in Nettle. The same algorithm was later also
103
42.8k
   adapted to GMP in mpn_sec_invert.
104
42.8k
105
42.8k
   It can be easily implemented in a way that does not depend on
106
42.8k
   secret branches or memory lookups, providing resistance against
107
42.8k
   some forms of side channel attack.
108
42.8k
109
42.8k
   There is also a description of the algorithm in Appendix 5 of "Fast
110
42.8k
   Software Polynomial Multiplication on ARM Processors using the NEON Engine"
111
42.8k
   by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo
112
42.8k
   Dahab in LNCS 8182
113
42.8k
      https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf
114
42.8k
115
42.8k
   Thanks to Niels for creating the algorithm, explaining some things
116
42.8k
   about it, and the reference to the paper.
117
42.8k
   */
118
42.8k
119
42.8k
   const size_t mod_words = mod.sig_words();
120
42.8k
   BOTAN_ASSERT(mod_words > 0, "Not empty");
121
42.8k
122
42.8k
   secure_vector<word> tmp_mem(5*mod_words);
123
42.8k
124
42.8k
   word* v_w = &tmp_mem[0];
125
42.8k
   word* u_w = &tmp_mem[1*mod_words];
126
42.8k
   word* b_w = &tmp_mem[2*mod_words];
127
42.8k
   word* a_w = &tmp_mem[3*mod_words];
128
42.8k
   word* mp1o2 = &tmp_mem[4*mod_words];
129
42.8k
130
42.8k
   CT::poison(tmp_mem.data(), tmp_mem.size());
131
42.8k
132
42.8k
   copy_mem(a_w, n.data(), std::min(n.size(), mod_words));
133
42.8k
   copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words));
134
42.8k
   u_w[0] = 1;
135
42.8k
   // v_w = 0
136
42.8k
137
42.8k
   // compute (mod + 1) / 2 which [because mod is odd] is equal to
138
42.8k
   // (mod / 2) + 1
139
42.8k
   copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words));
140
42.8k
   bigint_shr1(mp1o2, mod_words, 0, 1);
141
42.8k
   word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1);
142
42.8k
   BOTAN_ASSERT_NOMSG(carry == 0);
143
42.8k
144
42.8k
   // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
145
42.8k
   const size_t execs = 2 * mod.bits();
146
42.8k
147
38.2M
   for(size_t i = 0; i != execs; ++i)
148
38.1M
      {
149
38.1M
      const word odd_a = a_w[0] & 1;
150
38.1M
151
38.1M
      //if(odd_a) a -= b
152
38.1M
      word underflow = bigint_cnd_sub(odd_a, a_w, b_w, mod_words);
153
38.1M
154
38.1M
      //if(underflow) { b -= a; a = abs(a); swap(u, v); }
155
38.1M
      bigint_cnd_add(underflow, b_w, a_w, mod_words);
156
38.1M
      bigint_cnd_abs(underflow, a_w, mod_words);
157
38.1M
      bigint_cnd_swap(underflow, u_w, v_w, mod_words);
158
38.1M
159
38.1M
      // a >>= 1
160
38.1M
      bigint_shr1(a_w, mod_words, 0, 1);
161
38.1M
162
38.1M
      //if(odd_a) u -= v;
163
38.1M
      word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words);
164
38.1M
165
38.1M
      // if(borrow) u += p
166
38.1M
      bigint_cnd_add(borrow, u_w, mod.data(), mod_words);
167
38.1M
168
38.1M
      const word odd_u = u_w[0] & 1;
169
38.1M
170
38.1M
      // u >>= 1
171
38.1M
      bigint_shr1(u_w, mod_words, 0, 1);
172
38.1M
173
38.1M
      //if(odd_u) u += mp1o2;
174
38.1M
      bigint_cnd_add(odd_u, u_w, mp1o2, mod_words);
175
38.1M
      }
176
42.8k
177
42.8k
   auto a_is_0 = CT::Mask<word>::set();
178
346k
   for(size_t i = 0; i != mod_words; ++i)
179
303k
      a_is_0 &= CT::Mask<word>::is_zero(a_w[i]);
180
42.8k
181
42.8k
   auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1);
182
303k
   for(size_t i = 1; i != mod_words; ++i)
183
260k
      b_is_1 &= CT::Mask<word>::is_zero(b_w[i]);
184
42.8k
185
42.8k
   BOTAN_ASSERT(a_is_0.is_set(), "A is zero");
186
42.8k
187
42.8k
   // if b != 1 then gcd(n,mod) > 1 and inverse does not exist
188
42.8k
   // in which case zero out the result to indicate this
189
42.8k
   (~b_is_1).if_set_zero_out(v_w, mod_words);
190
42.8k
191
42.8k
   /*
192
42.8k
   * We've placed the result in the lowest words of the temp buffer.
193
42.8k
   * So just clear out the other values and then give that buffer to a
194
42.8k
   * BigInt.
195
42.8k
   */
196
42.8k
   clear_mem(&tmp_mem[mod_words], 4*mod_words);
197
42.8k
198
42.8k
   CT::unpoison(tmp_mem.data(), tmp_mem.size());
199
42.8k
200
42.8k
   BigInt r;
201
42.8k
   r.swap_reg(tmp_mem);
202
42.8k
   return r;
203
42.8k
   }
204
205
BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
206
1.09k
   {
207
1.09k
   /*
208
1.09k
   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
209
1.09k
   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
210
1.09k
   */
211
1.09k
212
1.09k
   if(a1.is_even())
213
0
      return 0;
214
1.09k
215
1.09k
   BigInt a = a1;
216
1.09k
   a.mask_bits(k);
217
1.09k
218
1.09k
   BigInt b = 1;
219
1.09k
   BigInt X = 0;
220
1.09k
   BigInt newb;
221
1.09k
222
1.09k
   const size_t a_words = a.sig_words();
223
1.09k
224
1.09k
   X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS);
225
1.09k
   b.grow_to(a_words);
226
1.09k
227
1.09k
   /*
228
1.09k
   Hide the exact value of k. k is anyway known to word length
229
1.09k
   granularity because of the length of a, so no point in doing more
230
1.09k
   than this.
231
1.09k
   */
232
1.09k
   const size_t iter = round_up(k, BOTAN_MP_WORD_BITS);
233
1.09k
234
602k
   for(size_t i = 0; i != iter; ++i)
235
601k
      {
236
601k
      const bool b0 = b.get_bit(0);
237
601k
      X.conditionally_set_bit(i, b0);
238
601k
      newb = b - a;
239
601k
      b.ct_cond_assign(b0, newb);
240
601k
      b >>= 1;
241
601k
      }
242
1.09k
243
1.09k
   X.mask_bits(k);
244
1.09k
   X.const_time_unpoison();
245
1.09k
   return X;
246
1.09k
   }
247
248
}
249
250
BigInt inverse_mod(const BigInt& n, const BigInt& mod)
251
43.0k
   {
252
43.0k
   if(mod.is_zero())
253
0
      throw BigInt::DivideByZero();
254
43.0k
   if(mod.is_negative() || n.is_negative())
255
0
      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
256
43.0k
   if(n.is_zero() || (n.is_even() && mod.is_even()))
257
179
      return 0;
258
42.9k
259
42.9k
   if(mod.is_odd())
260
42.2k
      {
261
42.2k
      /*
262
42.2k
      Fastpath for common case. This leaks information if n > mod
263
42.2k
      but we don't guarantee const time behavior in that case.
264
42.2k
      */
265
42.2k
      if(n < mod)
266
42.2k
         return inverse_mod_odd_modulus(n, mod);
267
65
      else
268
65
         return inverse_mod_odd_modulus(ct_modulo(n, mod), mod);
269
618
      }
270
618
271
618
   const size_t mod_lz = low_zero_bits(mod);
272
618
   BOTAN_ASSERT_NOMSG(mod_lz > 0);
273
618
   const size_t mod_bits = mod.bits();
274
618
   BOTAN_ASSERT_NOMSG(mod_bits > mod_lz);
275
618
276
618
   if(mod_lz == mod_bits - 1)
277
54
      {
278
54
      // In this case we are performing an inversion modulo 2^k
279
54
      return inverse_mod_pow2(n, mod_lz);
280
54
      }
281
564
282
564
   /*
283
564
   * In this case we are performing an inversion modulo 2^k*o for
284
564
   * some k > 1 and some odd (not necessarily prime) integer.
285
564
   * Compute the inversions modulo 2^k and modulo o, then combine them
286
564
   * using CRT, which is possible because 2^k and o are relatively prime.
287
564
   */
288
564
289
564
   const BigInt o = mod >> mod_lz;
290
564
   const BigInt n_redc = ct_modulo(n, o);
291
564
   const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o);
292
564
   const BigInt inv_2k = inverse_mod_pow2(n, mod_lz);
293
564
294
564
   // No modular inverse in this case:
295
564
   if(inv_o == 0 || inv_2k == 0)
296
88
      return 0;
297
476
298
476
   const BigInt m2k = BigInt::power_of_2(mod_lz);
299
476
   // Compute the CRT parameter
300
476
   const BigInt c = inverse_mod_pow2(o, mod_lz);
301
476
302
476
   // Compute h = c*(inv_2k-inv_o) mod 2^k
303
476
   BigInt h = c * (inv_2k - inv_o);
304
476
   const bool h_neg = h.is_negative();
305
476
   h.set_sign(BigInt::Positive);
306
476
   h.mask_bits(mod_lz);
307
476
   const bool h_nonzero = h.is_nonzero();
308
476
   h.ct_cond_assign(h_nonzero && h_neg, m2k - h);
309
476
310
476
   // Return result inv_o + h * o
311
476
   h *= o;
312
476
   h += inv_o;
313
476
   return h;
314
476
   }
315
316
// Deprecated forwarding functions:
317
BigInt inverse_euclid(const BigInt& x, const BigInt& modulus)
318
0
   {
319
0
   return inverse_mod(x, modulus);
320
0
   }
321
322
BigInt ct_inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
323
0
   {
324
0
   return inverse_mod_odd_modulus(n, mod);
325
0
   }
326
327
word monty_inverse(word a)
328
76.8k
   {
329
76.8k
   if(a % 2 == 0)
330
0
      throw Invalid_Argument("monty_inverse only valid for odd integers");
331
76.8k
332
76.8k
   /*
333
76.8k
   * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç
334
76.8k
   * https://eprint.iacr.org/2017/411.pdf sections 5 and 7.
335
76.8k
   */
336
76.8k
337
76.8k
   word b = 1;
338
76.8k
   word r = 0;
339
76.8k
340
4.99M
   for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i)
341
4.91M
      {
342
4.91M
      const word bi = b % 2;
343
4.91M
      r >>= 1;
344
4.91M
      r += bi << (BOTAN_MP_WORD_BITS - 1);
345
4.91M
346
4.91M
      b -= a * bi;
347
4.91M
      b >>= 1;
348
4.91M
      }
349
76.8k
350
76.8k
   // Now invert in addition space
351
76.8k
   r = (MP_WORD_MAX - r) + 1;
352
76.8k
353
76.8k
   return r;
354
76.8k
   }
355
356
}