Coverage Report

Created: 2021-02-21 07:20

/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.7k
   {
19
   // Caller should assure these preconditions:
20
57.7k
   BOTAN_DEBUG_ASSERT(n.is_positive());
21
57.7k
   BOTAN_DEBUG_ASSERT(mod.is_positive());
22
57.7k
   BOTAN_DEBUG_ASSERT(n < mod);
23
57.7k
   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.7k
   const size_t mod_words = mod.sig_words();
45
57.7k
   BOTAN_ASSERT(mod_words > 0, "Not empty");
46
47
57.7k
   secure_vector<word> tmp_mem(5*mod_words);
48
49
57.7k
   word* v_w = &tmp_mem[0];
50
57.7k
   word* u_w = &tmp_mem[1*mod_words];
51
57.7k
   word* b_w = &tmp_mem[2*mod_words];
52
57.7k
   word* a_w = &tmp_mem[3*mod_words];
53
57.7k
   word* mp1o2 = &tmp_mem[4*mod_words];
54
55
57.7k
   CT::poison(tmp_mem.data(), tmp_mem.size());
56
57
57.7k
   copy_mem(a_w, n.data(), std::min(n.size(), mod_words));
58
57.7k
   copy_mem(b_w, mod.data(), std::min(mod.size(), mod_words));
59
57.7k
   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.7k
   copy_mem(mp1o2, mod.data(), std::min(mod.size(), mod_words));
65
57.7k
   bigint_shr1(mp1o2, mod_words, 0, 1);
66
57.7k
   word carry = bigint_add2_nc(mp1o2, mod_words, u_w, 1);
67
57.7k
   BOTAN_ASSERT_NOMSG(carry == 0);
68
69
   // Only n.bits() + mod.bits() iterations are required, but avoid leaking the size of n
70
57.7k
   const size_t execs = 2 * mod.bits();
71
72
40.2M
   for(size_t i = 0; i != execs; ++i)
73
40.1M
      {
74
40.1M
      const word odd_a = a_w[0] & 1;
75
76
      //if(odd_a) a -= b
77
40.1M
      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
40.1M
      bigint_cnd_add(underflow, b_w, a_w, mod_words);
81
40.1M
      bigint_cnd_abs(underflow, a_w, mod_words);
82
40.1M
      bigint_cnd_swap(underflow, u_w, v_w, mod_words);
83
84
      // a >>= 1
85
40.1M
      bigint_shr1(a_w, mod_words, 0, 1);
86
87
      //if(odd_a) u -= v;
88
40.1M
      word borrow = bigint_cnd_sub(odd_a, u_w, v_w, mod_words);
89
90
      // if(borrow) u += p
91
40.1M
      bigint_cnd_add(borrow, u_w, mod.data(), mod_words);
92
93
40.1M
      const word odd_u = u_w[0] & 1;
94
95
      // u >>= 1
96
40.1M
      bigint_shr1(u_w, mod_words, 0, 1);
97
98
      //if(odd_u) u += mp1o2;
99
40.1M
      bigint_cnd_add(odd_u, u_w, mp1o2, mod_words);
100
40.1M
      }
101
102
57.7k
   auto a_is_0 = CT::Mask<word>::set();
103
375k
   for(size_t i = 0; i != mod_words; ++i)
104
317k
      a_is_0 &= CT::Mask<word>::is_zero(a_w[i]);
105
106
57.7k
   auto b_is_1 = CT::Mask<word>::is_equal(b_w[0], 1);
107
317k
   for(size_t i = 1; i != mod_words; ++i)
108
259k
      b_is_1 &= CT::Mask<word>::is_zero(b_w[i]);
109
110
57.7k
   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.7k
   (~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.7k
   clear_mem(&tmp_mem[mod_words], 4*mod_words);
122
123
57.7k
   CT::unpoison(tmp_mem.data(), tmp_mem.size());
124
125
57.7k
   BigInt r;
126
57.7k
   r.swap_reg(tmp_mem);
127
57.7k
   return r;
128
57.7k
   }
129
130
BigInt inverse_mod_pow2(const BigInt& a1, size_t k)
131
1.28k
   {
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.28k
   if(a1.is_even())
138
0
      return 0;
139
140
1.28k
   BigInt a = a1;
141
1.28k
   a.mask_bits(k);
142
143
1.28k
   BigInt b = 1;
144
1.28k
   BigInt X = 0;
145
1.28k
   BigInt newb;
146
147
1.28k
   const size_t a_words = a.sig_words();
148
149
1.28k
   X.grow_to(round_up(k, BOTAN_MP_WORD_BITS) / BOTAN_MP_WORD_BITS);
150
1.28k
   b.grow_to(a_words);
151
152
   /*
153
   Hide the exact value of k. k is anyway known to word length
154
   granularity because of the length of a, so no point in doing more
155
   than this.
156
   */
157
1.28k
   const size_t iter = round_up(k, BOTAN_MP_WORD_BITS);
158
159
721k
   for(size_t i = 0; i != iter; ++i)
160
719k
      {
161
719k
      const bool b0 = b.get_bit(0);
162
719k
      X.conditionally_set_bit(i, b0);
163
719k
      newb = b - a;
164
719k
      b.ct_cond_assign(b0, newb);
165
719k
      b >>= 1;
166
719k
      }
167
168
1.28k
   X.mask_bits(k);
169
1.28k
   X.const_time_unpoison();
170
1.28k
   return X;
171
1.28k
   }
172
173
}
174
175
BigInt inverse_mod(const BigInt& n, const BigInt& mod)
176
57.8k
   {
177
57.8k
   if(mod.is_zero())
178
0
      throw Invalid_Argument("inverse_mod modulus cannot be zero");
179
57.8k
   if(mod.is_negative() || n.is_negative())
180
0
      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
181
57.8k
   if(n.is_zero() || (n.is_even() && mod.is_even()))
182
7
      return 0;
183
184
57.8k
   if(mod.is_odd())
185
57.1k
      {
186
      /*
187
      Fastpath for common case. This leaks if n is greater than mod or
188
      not, but we don't guarantee const time behavior in that case.
189
      */
190
57.1k
      if(n < mod)
191
57.0k
         return inverse_mod_odd_modulus(n, mod);
192
91
      else
193
91
         return inverse_mod_odd_modulus(ct_modulo(n, mod), mod);
194
716
      }
195
196
716
   const size_t mod_lz = low_zero_bits(mod);
197
716
   BOTAN_ASSERT_NOMSG(mod_lz > 0);
198
716
   const size_t mod_bits = mod.bits();
199
716
   BOTAN_ASSERT_NOMSG(mod_bits > mod_lz);
200
201
716
   if(mod_lz == mod_bits - 1)
202
50
      {
203
      // In this case we are performing an inversion modulo 2^k
204
50
      return inverse_mod_pow2(n, mod_lz);
205
50
      }
206
207
   /*
208
   * In this case we are performing an inversion modulo 2^k*o for
209
   * some k > 1 and some odd (not necessarily prime) integer.
210
   * Compute the inversions modulo 2^k and modulo o, then combine them
211
   * using CRT, which is possible because 2^k and o are relatively prime.
212
   */
213
214
666
   const BigInt o = mod >> mod_lz;
215
666
   const BigInt n_redc = ct_modulo(n, o);
216
666
   const BigInt inv_o = inverse_mod_odd_modulus(n_redc, o);
217
666
   const BigInt inv_2k = inverse_mod_pow2(n, mod_lz);
218
219
   // No modular inverse in this case:
220
666
   if(inv_o == 0 || inv_2k == 0)
221
100
      return 0;
222
223
566
   const BigInt m2k = BigInt::power_of_2(mod_lz);
224
   // Compute the CRT parameter
225
566
   const BigInt c = inverse_mod_pow2(o, mod_lz);
226
227
   // Compute h = c*(inv_2k-inv_o) mod 2^k
228
566
   BigInt h = c * (inv_2k - inv_o);
229
566
   const bool h_neg = h.is_negative();
230
566
   h.set_sign(BigInt::Positive);
231
566
   h.mask_bits(mod_lz);
232
566
   const bool h_nonzero = h.is_nonzero();
233
566
   h.ct_cond_assign(h_nonzero && h_neg, m2k - h);
234
235
   // Return result inv_o + h * o
236
566
   h *= o;
237
566
   h += inv_o;
238
566
   return h;
239
566
   }
240
241
}