/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 | 57.6k | { |
19 | | // Caller should assure these preconditions: |
20 | 57.6k | BOTAN_DEBUG_ASSERT(n.is_positive()); |
21 | 57.6k | BOTAN_DEBUG_ASSERT(mod.is_positive()); |
22 | 57.6k | BOTAN_DEBUG_ASSERT(n < mod); |
23 | 57.6k | 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 | 57.6k | const size_t mod_words = mod.sig_words(); |
45 | 57.6k | BOTAN_ASSERT(mod_words > 0, "Not empty"); |
46 | | |
47 | 57.6k | secure_vector<word> tmp_mem(5*mod_words); |
48 | | |
49 | 57.6k | word* v_w = &tmp_mem[0]; |
50 | 57.6k | word* u_w = &tmp_mem[1*mod_words]; |
51 | 57.6k | word* b_w = &tmp_mem[2*mod_words]; |
52 | 57.6k | word* a_w = &tmp_mem[3*mod_words]; |
53 | 57.6k | word* mp1o2 = &tmp_mem[4*mod_words]; |
54 | | |
55 | 57.6k | CT::poison(tmp_mem.data(), tmp_mem.size()); |
56 | | |
57 | 57.6k | copy_mem(a_w, n.data(), std::min(n.size(), mod_words)); |
58 | 57.6k | copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words)); |
59 | 57.6k | 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 | 57.6k | copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words)); |
65 | 57.6k | bigint_shr1(mp1o2, mod_words, 0, 1); |
66 | 57.6k | word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1); |
67 | 57.6k | BOTAN_ASSERT_NOMSG(carry == 0); |
68 | | |
69 | | // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n |
70 | 57.6k | const size_t execs = 2 * mod.bits(); |
71 | | |
72 | 41.3M | for(size_t i = 0; i != execs; ++i) |
73 | 41.2M | { |
74 | 41.2M | const word odd_a = a_w[0] & 1; |
75 | | |
76 | | //if(odd_a) a -= b |
77 | 41.2M | 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 | 41.2M | bigint_cnd_add(underflow, b_w, a_w, mod_words); |
81 | 41.2M | bigint_cnd_abs(underflow, a_w, mod_words); |
82 | 41.2M | bigint_cnd_swap(underflow, u_w, v_w, mod_words); |
83 | | |
84 | | // a >>= 1 |
85 | 41.2M | bigint_shr1(a_w, mod_words, 0, 1); |
86 | | |
87 | | //if(odd_a) u -= v; |
88 | 41.2M | word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words); |
89 | | |
90 | | // if(borrow) u += p |
91 | 41.2M | bigint_cnd_add(borrow, u_w, mod.data(), mod_words); |
92 | | |
93 | 41.2M | const word odd_u = u_w[0] & 1; |
94 | | |
95 | | // u >>= 1 |
96 | 41.2M | bigint_shr1(u_w, mod_words, 0, 1); |
97 | | |
98 | | //if(odd_u) u += mp1o2; |
99 | 41.2M | bigint_cnd_add(odd_u, u_w, mp1o2, mod_words); |
100 | 41.2M | } |
101 | | |
102 | 57.6k | auto a_is_0 = CT::Mask<word>::set(); |
103 | 384k | for(size_t i = 0; i != mod_words; ++i) |
104 | 326k | a_is_0 &= CT::Mask<word>::is_zero(a_w[i]); |
105 | | |
106 | 57.6k | auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1); |
107 | 326k | for(size_t i = 1; i != mod_words; ++i) |
108 | 269k | b_is_1 &= CT::Mask<word>::is_zero(b_w[i]); |
109 | | |
110 | 57.6k | 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 | 57.6k | (~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 | 57.6k | clear_mem(&tmp_mem[mod_words], 4*mod_words); |
122 | | |
123 | 57.6k | CT::unpoison(tmp_mem.data(), tmp_mem.size()); |
124 | | |
125 | 57.6k | BigInt r; |
126 | 57.6k | r.swap_reg(tmp_mem); |
127 | 57.6k | return r; |
128 | 57.6k | } |
129 | | |
130 | | BigInt inverse_mod_pow2(const BigInt& a1, size_t k) |
131 | 1.12k | { |
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 | 1.12k | if(a1.is_even() || k == 0) |
138 | 0 | return 0; |
139 | 1.12k | if(k == 1) |
140 | 8 | return 1; |
141 | | |
142 | 1.11k | BigInt a = a1; |
143 | 1.11k | a.mask_bits(k); |
144 | | |
145 | 1.11k | BigInt b = 1; |
146 | 1.11k | BigInt X = 0; |
147 | 1.11k | BigInt newb; |
148 | | |
149 | 1.11k | const size_t a_words = a.sig_words(); |
150 | | |
151 | 1.11k | X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS); |
152 | 1.11k | 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 | 1.11k | const size_t iter = round_up(k, BOTAN_MP_WORD_BITS); |
160 | | |
161 | 674k | for(size_t i = 0; i != iter; ++i) |
162 | 673k | { |
163 | 673k | const bool b0 = b.get_bit(0); |
164 | 673k | X.conditionally_set_bit(i, b0); |
165 | 673k | newb = b - a; |
166 | 673k | b.ct_cond_assign(b0, newb); |
167 | 673k | b >>= 1; |
168 | 673k | } |
169 | | |
170 | 1.11k | X.mask_bits(k); |
171 | 1.11k | X.const_time_unpoison(); |
172 | 1.11k | return X; |
173 | 1.11k | } |
174 | | |
175 | | } |
176 | | |
177 | | BigInt inverse_mod(const BigInt& n, const BigInt& mod) |
178 | 57.6k | { |
179 | 57.6k | if(mod.is_zero()) |
180 | 0 | throw Invalid_Argument("inverse_mod modulus cannot be zero"); |
181 | 57.6k | if(mod.is_negative() || n.is_negative()) |
182 | 0 | throw Invalid_Argument("inverse_mod: arguments must be non-negative"); |
183 | 57.6k | if(n.is_zero() || (n.is_even() && mod.is_even())) |
184 | 5 | return 0; |
185 | | |
186 | 57.6k | if(mod.is_odd()) |
187 | 56.9k | { |
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 | 56.9k | if(n < mod) |
193 | 56.8k | return inverse_mod_odd_modulus(n, mod); |
194 | 83 | else |
195 | 83 | return inverse_mod_odd_modulus(ct_modulo(n, mod), mod); |
196 | 754 | } |
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 | 754 | BOTAN_DEBUG_ASSERT(n.is_odd()); |
201 | | |
202 | 754 | const size_t mod_lz = low_zero_bits(mod); |
203 | 754 | BOTAN_ASSERT_NOMSG(mod_lz > 0); |
204 | 754 | const size_t mod_bits = mod.bits(); |
205 | 754 | BOTAN_ASSERT_NOMSG(mod_bits > mod_lz); |
206 | | |
207 | 754 | if(mod_lz == mod_bits - 1) |
208 | 51 | { |
209 | | // In this case we are performing an inversion modulo 2^k |
210 | 51 | return inverse_mod_pow2(n, mod_lz); |
211 | 51 | } |
212 | | |
213 | 703 | if(mod_lz == 1) |
214 | 120 | { |
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 | 120 | const BigInt o = mod >> 1; |
230 | 120 | const BigInt n_redc = ct_modulo(n, o); |
231 | 120 | const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o); |
232 | | |
233 | | // No modular inverse in this case: |
234 | 120 | if(inv_o == 0) |
235 | 15 | return 0; |
236 | | |
237 | 105 | BigInt h = inv_o; |
238 | 105 | h.ct_cond_add(!inv_o.get_bit(0), o); |
239 | 105 | return h; |
240 | 105 | } |
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 | 583 | const BigInt o = mod >> mod_lz; |
250 | 583 | const BigInt n_redc = ct_modulo(n, o); |
251 | 583 | const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o); |
252 | 583 | const BigInt inv_2k = inverse_mod_pow2(n, mod_lz); |
253 | | |
254 | | // No modular inverse in this case: |
255 | 583 | if(inv_o == 0 || inv_2k == 0) |
256 | 96 | return 0; |
257 | | |
258 | 487 | const BigInt m2k = BigInt::power_of_2(mod_lz); |
259 | | // Compute the CRT parameter |
260 | 487 | const BigInt c = inverse_mod_pow2(o, mod_lz); |
261 | | |
262 | | // Compute h = c*(inv_2k-inv_o) mod 2^k |
263 | 487 | BigInt h = c * (inv_2k - inv_o); |
264 | 487 | const bool h_neg = h.is_negative(); |
265 | 487 | h.set_sign(BigInt::Positive); |
266 | 487 | h.mask_bits(mod_lz); |
267 | 487 | const bool h_nonzero = h.is_nonzero(); |
268 | 487 | h.ct_cond_assign(h_nonzero && h_neg, m2k - h); |
269 | | |
270 | | // Return result inv_o + h * o |
271 | 487 | h *= o; |
272 | 487 | h += inv_o; |
273 | 487 | return h; |
274 | 487 | } |
275 | | |
276 | | } |