Botan  2.15.0
Crypto and TLS for C++11
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 
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 */
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 
90 namespace {
91 
92 BigInt inverse_mod_odd_modulus(const BigInt& n, const BigInt& mod)
93  {
94  // Caller should assure these preconditions:
95  BOTAN_DEBUG_ASSERT(n.is_positive());
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 
205 BigInt 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 
250 BigInt 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:
318  {
319  return inverse_mod(x, modulus);
320  }
321 
323  {
324  return inverse_mod_odd_modulus(n, mod);
325  }
326 
327 word 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 }
void bigint_shr1(word x[], size_t x_size, size_t word_shift, size_t bit_shift)
Definition: mp_core.h:427
fe X
Definition: ge.cpp:27
bool is_negative() const
Definition: bigint.h:527
size_t low_zero_bits(const BigInt &n)
Definition: numthry.cpp:39
void carry(int64_t &h0, int64_t &h1)
size_t bits() const
Definition: bigint.cpp:297
void clear_mem(T *ptr, size_t n)
Definition: mem_ops.h:115
bool is_zero() const
Definition: bigint.h:421
BigInt ct_inverse_mod_odd_modulus(const BigInt &n, const BigInt &mod)
Definition: mod_inv.cpp:322
BigInt normalized_montgomery_inverse(const BigInt &a, const BigInt &p)
Definition: mod_inv.cpp:75
#define BOTAN_ASSERT_NOMSG(expr)
Definition: assert.h:68
BigInt ct_modulo(const BigInt &x, const BigInt &y)
Definition: divide.cpp:118
void poison(const T *p, size_t n)
Definition: ct_utils.h:48
bool is_even() const
Definition: bigint.h:403
bool is_nonzero() const
Definition: bigint.h:415
void mask_bits(size_t n)
Definition: bigint.h:455
const word * data() const
Definition: bigint.h:620
#define BOTAN_ASSERT(expr, assertion_made)
Definition: assert.h:55
BigInt const BigInt & mod
Definition: numthry.h:105
word bigint_cnd_sub(word cnd, word x[], size_t x_size, const word y[], size_t y_size)
Definition: mp_core.h:88
#define BOTAN_DEBUG_ASSERT(expr)
Definition: assert.h:123
void bigint_cnd_swap(word cnd, word x[], word y[], size_t size)
Definition: mp_core.h:29
BigInt const BigInt & modulus
Definition: numthry.h:100
word bigint_cnd_add(word cnd, word x[], word x_size, const word y[], size_t y_size)
Definition: mp_core.h:42
BigInt inverse_euclid(const BigInt &x, const BigInt &modulus)
Definition: mod_inv.cpp:317
bool is_odd() const
Definition: bigint.h:409
BigInt inverse_mod(const BigInt &n, const BigInt &mod)
Definition: mod_inv.cpp:250
size_t size() const
Definition: bigint.h:580
word bigint_add2_nc(word x[], size_t x_size, const word y[], size_t y_size)
Definition: mp_core.h:227
static BigInt power_of_2(size_t n)
Definition: bigint.h:758
void copy_mem(T *out, const T *in, size_t n)
Definition: mem_ops.h:133
Definition: alg_id.cpp:13
size_t sig_words() const
Definition: bigint.h:586
const word MP_WORD_MAX
Definition: mp_core.h:22
word monty_inverse(word a)
Definition: mod_inv.cpp:327
size_t almost_montgomery_inverse(BigInt &result, const BigInt &a, const BigInt &p)
Definition: mod_inv.cpp:27
void unpoison(const T *p, size_t n)
Definition: ct_utils.h:59
static Mask< T > set()
Definition: ct_utils.h:107
void bigint_cnd_abs(word cnd, word x[], size_t size)
Definition: mp_core.h:212
static Mask< T > is_zero(T x)
Definition: ct_utils.h:141
size_t round_up(size_t n, size_t align_to)
Definition: rounding.h:21
void set_sign(Sign sign)
Definition: bigint.h:563
static Mask< T > is_equal(T x, T y)
Definition: ct_utils.h:149
bool is_positive() const
Definition: bigint.h:533
void ct_cond_assign(bool predicate, const BigInt &other)
Definition: bigint.cpp:489