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