/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 | 51.1k | { |
94 | 51.1k | // Caller should assure these preconditions: |
95 | 51.1k | BOTAN_DEBUG_ASSERT(n.is_positive()); |
96 | 51.1k | BOTAN_DEBUG_ASSERT(mod.is_positive()); |
97 | 51.1k | BOTAN_DEBUG_ASSERT(n < mod); |
98 | 51.1k | BOTAN_DEBUG_ASSERT(mod >= 3 && mod.is_odd()); |
99 | 51.1k | |
100 | 51.1k | /* |
101 | 51.1k | This uses a modular inversion algorithm designed by Niels Möller |
102 | 51.1k | and implemented in Nettle. The same algorithm was later also |
103 | 51.1k | adapted to GMP in mpn_sec_invert. |
104 | 51.1k | |
105 | 51.1k | It can be easily implemented in a way that does not depend on |
106 | 51.1k | secret branches or memory lookups, providing resistance against |
107 | 51.1k | some forms of side channel attack. |
108 | 51.1k | |
109 | 51.1k | There is also a description of the algorithm in Appendix 5 of "Fast |
110 | 51.1k | Software Polynomial Multiplication on ARM Processors using the NEON Engine" |
111 | 51.1k | by Danilo Câmara, Conrado P. L. Gouvêa, Julio López, and Ricardo |
112 | 51.1k | Dahab in LNCS 8182 |
113 | 51.1k | https://conradoplg.cryptoland.net/files/2010/12/mocrysen13.pdf |
114 | 51.1k | |
115 | 51.1k | Thanks to Niels for creating the algorithm, explaining some things |
116 | 51.1k | about it, and the reference to the paper. |
117 | 51.1k | */ |
118 | 51.1k | |
119 | 51.1k | const size_t mod_words = mod.sig_words(); |
120 | 51.1k | BOTAN_ASSERT(mod_words > 0, "Not empty"); |
121 | 51.1k | |
122 | 51.1k | secure_vector<word> tmp_mem(5*mod_words); |
123 | 51.1k | |
124 | 51.1k | word* v_w = &tmp_mem[0]; |
125 | 51.1k | word* u_w = &tmp_mem[1*mod_words]; |
126 | 51.1k | word* b_w = &tmp_mem[2*mod_words]; |
127 | 51.1k | word* a_w = &tmp_mem[3*mod_words]; |
128 | 51.1k | word* mp1o2 = &tmp_mem[4*mod_words]; |
129 | 51.1k | |
130 | 51.1k | CT::poison(tmp_mem.data(), tmp_mem.size()); |
131 | 51.1k | |
132 | 51.1k | copy_mem(a_w, n.data(), std::min(n.size(), mod_words)); |
133 | 51.1k | copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words)); |
134 | 51.1k | u_w[0] = 1; |
135 | 51.1k | // v_w = 0 |
136 | 51.1k | |
137 | 51.1k | // compute (mod + 1) / 2 which [because mod is odd] is equal to |
138 | 51.1k | // (mod / 2) + 1 |
139 | 51.1k | copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words)); |
140 | 51.1k | bigint_shr1(mp1o2, mod_words, 0, 1); |
141 | 51.1k | word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1); |
142 | 51.1k | BOTAN_ASSERT_NOMSG(carry == 0); |
143 | 51.1k | |
144 | 51.1k | // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n |
145 | 51.1k | const size_t execs = 2 * mod.bits(); |
146 | 51.1k | |
147 | 43.3M | for(size_t i = 0; i != execs; ++i) |
148 | 43.3M | { |
149 | 43.3M | const word odd_a = a_w[0] & 1; |
150 | 43.3M | |
151 | 43.3M | //if(odd_a) a -= b |
152 | 43.3M | word underflow = bigint_cnd_sub(odd_a, a_w, b_w, mod_words); |
153 | 43.3M | |
154 | 43.3M | //if(underflow) { b -= a; a = abs(a); swap(u, v); } |
155 | 43.3M | bigint_cnd_add(underflow, b_w, a_w, mod_words); |
156 | 43.3M | bigint_cnd_abs(underflow, a_w, mod_words); |
157 | 43.3M | bigint_cnd_swap(underflow, u_w, v_w, mod_words); |
158 | 43.3M | |
159 | 43.3M | // a >>= 1 |
160 | 43.3M | bigint_shr1(a_w, mod_words, 0, 1); |
161 | 43.3M | |
162 | 43.3M | //if(odd_a) u -= v; |
163 | 43.3M | word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words); |
164 | 43.3M | |
165 | 43.3M | // if(borrow) u += p |
166 | 43.3M | bigint_cnd_add(borrow, u_w, mod.data(), mod_words); |
167 | 43.3M | |
168 | 43.3M | const word odd_u = u_w[0] & 1; |
169 | 43.3M | |
170 | 43.3M | // u >>= 1 |
171 | 43.3M | bigint_shr1(u_w, mod_words, 0, 1); |
172 | 43.3M | |
173 | 43.3M | //if(odd_u) u += mp1o2; |
174 | 43.3M | bigint_cnd_add(odd_u, u_w, mp1o2, mod_words); |
175 | 43.3M | } |
176 | 51.1k | |
177 | 51.1k | auto a_is_0 = CT::Mask<word>::set(); |
178 | 395k | for(size_t i = 0; i != mod_words; ++i) |
179 | 344k | a_is_0 &= CT::Mask<word>::is_zero(a_w[i]); |
180 | 51.1k | |
181 | 51.1k | auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1); |
182 | 344k | for(size_t i = 1; i != mod_words; ++i) |
183 | 293k | b_is_1 &= CT::Mask<word>::is_zero(b_w[i]); |
184 | 51.1k | |
185 | 51.1k | BOTAN_ASSERT(a_is_0.is_set(), "A is zero"); |
186 | 51.1k | |
187 | 51.1k | // if b != 1 then gcd(n,mod) > 1 and inverse does not exist |
188 | 51.1k | // in which case zero out the result to indicate this |
189 | 51.1k | (~b_is_1).if_set_zero_out(v_w, mod_words); |
190 | 51.1k | |
191 | 51.1k | /* |
192 | 51.1k | * We've placed the result in the lowest words of the temp buffer. |
193 | 51.1k | * So just clear out the other values and then give that buffer to a |
194 | 51.1k | * BigInt. |
195 | 51.1k | */ |
196 | 51.1k | clear_mem(&tmp_mem[mod_words], 4*mod_words); |
197 | 51.1k | |
198 | 51.1k | CT::unpoison(tmp_mem.data(), tmp_mem.size()); |
199 | 51.1k | |
200 | 51.1k | BigInt r; |
201 | 51.1k | r.swap_reg(tmp_mem); |
202 | 51.1k | return r; |
203 | 51.1k | } |
204 | | |
205 | | BigInt inverse_mod_pow2(const BigInt& a1, size_t k) |
206 | 1.23k | { |
207 | 1.23k | /* |
208 | 1.23k | * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç |
209 | 1.23k | * https://eprint.iacr.org/2017/411.pdf sections 5 and 7. |
210 | 1.23k | */ |
211 | 1.23k | |
212 | 1.23k | if(a1.is_even()) |
213 | 0 | return 0; |
214 | 1.23k | |
215 | 1.23k | BigInt a = a1; |
216 | 1.23k | a.mask_bits(k); |
217 | 1.23k | |
218 | 1.23k | BigInt b = 1; |
219 | 1.23k | BigInt X = 0; |
220 | 1.23k | BigInt newb; |
221 | 1.23k | |
222 | 1.23k | const size_t a_words = a.sig_words(); |
223 | 1.23k | |
224 | 1.23k | X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS); |
225 | 1.23k | b.grow_to(a_words); |
226 | 1.23k | |
227 | 1.23k | /* |
228 | 1.23k | Hide the exact value of k. k is anyway known to word length |
229 | 1.23k | granularity because of the length of a, so no point in doing more |
230 | 1.23k | than this. |
231 | 1.23k | */ |
232 | 1.23k | const size_t iter = round_up(k, BOTAN_MP_WORD_BITS); |
233 | 1.23k | |
234 | 748k | for(size_t i = 0; i != iter; ++i) |
235 | 747k | { |
236 | 747k | const bool b0 = b.get_bit(0); |
237 | 747k | X.conditionally_set_bit(i, b0); |
238 | 747k | newb = b - a; |
239 | 747k | b.ct_cond_assign(b0, newb); |
240 | 747k | b >>= 1; |
241 | 747k | } |
242 | 1.23k | |
243 | 1.23k | X.mask_bits(k); |
244 | 1.23k | X.const_time_unpoison(); |
245 | 1.23k | return X; |
246 | 1.23k | } |
247 | | |
248 | | } |
249 | | |
250 | | BigInt inverse_mod(const BigInt& n, const BigInt& mod) |
251 | 51.3k | { |
252 | 51.3k | if(mod.is_zero()) |
253 | 0 | throw BigInt::DivideByZero(); |
254 | 51.3k | if(mod.is_negative() || n.is_negative()) |
255 | 0 | throw Invalid_Argument("inverse_mod: arguments must be non-negative"); |
256 | 51.3k | if(n.is_zero() || (n.is_even() && mod.is_even())) |
257 | 163 | return 0; |
258 | 51.2k | |
259 | 51.2k | if(mod.is_odd()) |
260 | 50.5k | { |
261 | 50.5k | /* |
262 | 50.5k | Fastpath for common case. This leaks information if n > mod |
263 | 50.5k | but we don't guarantee const time behavior in that case. |
264 | 50.5k | */ |
265 | 50.5k | if(n < mod) |
266 | 50.4k | return inverse_mod_odd_modulus(n, mod); |
267 | 66 | else |
268 | 66 | return inverse_mod_odd_modulus(ct_modulo(n, mod), mod); |
269 | 701 | } |
270 | 701 | |
271 | 701 | const size_t mod_lz = low_zero_bits(mod); |
272 | 701 | BOTAN_ASSERT_NOMSG(mod_lz > 0); |
273 | 701 | const size_t mod_bits = mod.bits(); |
274 | 701 | BOTAN_ASSERT_NOMSG(mod_bits > mod_lz); |
275 | 701 | |
276 | 701 | if(mod_lz == mod_bits - 1) |
277 | 53 | { |
278 | 53 | // In this case we are performing an inversion modulo 2^k |
279 | 53 | return inverse_mod_pow2(n, mod_lz); |
280 | 53 | } |
281 | 648 | |
282 | 648 | /* |
283 | 648 | * In this case we are performing an inversion modulo 2^k*o for |
284 | 648 | * some k > 1 and some odd (not necessarily prime) integer. |
285 | 648 | * Compute the inversions modulo 2^k and modulo o, then combine them |
286 | 648 | * using CRT, which is possible because 2^k and o are relatively prime. |
287 | 648 | */ |
288 | 648 | |
289 | 648 | const BigInt o = mod >> mod_lz; |
290 | 648 | const BigInt n_redc = ct_modulo(n, o); |
291 | 648 | const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o); |
292 | 648 | const BigInt inv_2k = inverse_mod_pow2(n, mod_lz); |
293 | 648 | |
294 | 648 | // No modular inverse in this case: |
295 | 648 | if(inv_o == 0 || inv_2k == 0) |
296 | 111 | return 0; |
297 | 537 | |
298 | 537 | const BigInt m2k = BigInt::power_of_2(mod_lz); |
299 | 537 | // Compute the CRT parameter |
300 | 537 | const BigInt c = inverse_mod_pow2(o, mod_lz); |
301 | 537 | |
302 | 537 | // Compute h = c*(inv_2k-inv_o) mod 2^k |
303 | 537 | BigInt h = c * (inv_2k - inv_o); |
304 | 537 | const bool h_neg = h.is_negative(); |
305 | 537 | h.set_sign(BigInt::Positive); |
306 | 537 | h.mask_bits(mod_lz); |
307 | 537 | const bool h_nonzero = h.is_nonzero(); |
308 | 537 | h.ct_cond_assign(h_nonzero && h_neg, m2k - h); |
309 | 537 | |
310 | 537 | // Return result inv_o + h * o |
311 | 537 | h *= o; |
312 | 537 | h += inv_o; |
313 | 537 | return h; |
314 | 537 | } |
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 | 85.6k | { |
329 | 85.6k | if(a % 2 == 0) |
330 | 0 | throw Invalid_Argument("monty_inverse only valid for odd integers"); |
331 | 85.6k | |
332 | 85.6k | /* |
333 | 85.6k | * From "A New Algorithm for Inversion mod p^k" by Çetin Kaya Koç |
334 | 85.6k | * https://eprint.iacr.org/2017/411.pdf sections 5 and 7. |
335 | 85.6k | */ |
336 | 85.6k | |
337 | 85.6k | word b = 1; |
338 | 85.6k | word r = 0; |
339 | 85.6k | |
340 | 5.57M | for(size_t i = 0; i != BOTAN_MP_WORD_BITS; ++i) |
341 | 5.48M | { |
342 | 5.48M | const word bi = b % 2; |
343 | 5.48M | r >>= 1; |
344 | 5.48M | r += bi << (BOTAN_MP_WORD_BITS - 1); |
345 | 5.48M | |
346 | 5.48M | b -= a * bi; |
347 | 5.48M | b >>= 1; |
348 | 5.48M | } |
349 | 85.6k | |
350 | 85.6k | // Now invert in addition space |
351 | 85.6k | r = (MP_WORD_MAX - r) + 1; |
352 | 85.6k | |
353 | 85.6k | return r; |
354 | 85.6k | } |
355 | | |
356 | | } |