Branch data Line data Source code
1 : : /*
2 : : * Core bignum functions
3 : : *
4 : : * Copyright The Mbed TLS Contributors
5 : : * SPDX-License-Identifier: Apache-2.0 OR GPL-2.0-or-later
6 : : */
7 : :
8 : : #include "common.h"
9 : :
10 : : #if defined(MBEDTLS_BIGNUM_C)
11 : :
12 : : #include <string.h>
13 : :
14 : : #include "mbedtls/error.h"
15 : : #include "mbedtls/platform_util.h"
16 : : #include "constant_time_internal.h"
17 : :
18 : : #include "mbedtls/platform.h"
19 : :
20 : : #include "bignum_core.h"
21 : : #include "bignum_core_invasive.h"
22 : : #include "bn_mul.h"
23 : : #include "constant_time_internal.h"
24 : :
25 : 480830 : size_t mbedtls_mpi_core_clz(mbedtls_mpi_uint a)
26 : : {
27 : : #if defined(__has_builtin)
28 : : #if (MBEDTLS_MPI_UINT_MAX == UINT_MAX) && __has_builtin(__builtin_clz)
29 : : #define core_clz __builtin_clz
30 : : #elif (MBEDTLS_MPI_UINT_MAX == ULONG_MAX) && __has_builtin(__builtin_clzl)
31 : : #define core_clz __builtin_clzl
32 : : #elif (MBEDTLS_MPI_UINT_MAX == ULLONG_MAX) && __has_builtin(__builtin_clzll)
33 : : #define core_clz __builtin_clzll
34 : : #endif
35 : : #endif
36 : : #if defined(core_clz)
37 : 480830 : return (size_t) core_clz(a);
38 : : #else
39 : : size_t j;
40 : : mbedtls_mpi_uint mask = (mbedtls_mpi_uint) 1 << (biL - 1);
41 : :
42 : : for (j = 0; j < biL; j++) {
43 : : if (a & mask) {
44 : : break;
45 : : }
46 : :
47 : : mask >>= 1;
48 : : }
49 : :
50 : : return j;
51 : : #endif
52 : : }
53 : :
54 : 480830 : size_t mbedtls_mpi_core_bitlen(const mbedtls_mpi_uint *A, size_t A_limbs)
55 : : {
56 : 480830 : int i;
57 : 480830 : size_t j;
58 : :
59 [ + - ]: 3081099 : for (i = ((int) A_limbs) - 1; i >= 0; i--) {
60 [ + + ]: 3081099 : if (A[i] != 0) {
61 : 480830 : j = biL - mbedtls_mpi_core_clz(A[i]);
62 : 480830 : return (i * biL) + j;
63 : : }
64 : : }
65 : :
66 : : return 0;
67 : : }
68 : :
69 : 720 : static mbedtls_mpi_uint mpi_bigendian_to_host(mbedtls_mpi_uint a)
70 : : {
71 : 720 : if (MBEDTLS_IS_BIG_ENDIAN) {
72 : : /* Nothing to do on bigendian systems. */
73 : : return a;
74 : : } else {
75 : : #if defined(MBEDTLS_HAVE_INT32)
76 : 720 : return (mbedtls_mpi_uint) MBEDTLS_BSWAP32(a);
77 : : #elif defined(MBEDTLS_HAVE_INT64)
78 : : return (mbedtls_mpi_uint) MBEDTLS_BSWAP64(a);
79 : : #endif
80 : : }
81 : : }
82 : :
83 : 90 : void mbedtls_mpi_core_bigendian_to_host(mbedtls_mpi_uint *A,
84 : : size_t A_limbs)
85 : : {
86 : 90 : mbedtls_mpi_uint *cur_limb_left;
87 : 90 : mbedtls_mpi_uint *cur_limb_right;
88 [ + - ]: 90 : if (A_limbs == 0) {
89 : : return;
90 : : }
91 : :
92 : : /*
93 : : * Traverse limbs and
94 : : * - adapt byte-order in each limb
95 : : * - swap the limbs themselves.
96 : : * For that, simultaneously traverse the limbs from left to right
97 : : * and from right to left, as long as the left index is not bigger
98 : : * than the right index (it's not a problem if limbs is odd and the
99 : : * indices coincide in the last iteration).
100 : : */
101 : 90 : for (cur_limb_left = A, cur_limb_right = A + (A_limbs - 1);
102 [ + + ]: 450 : cur_limb_left <= cur_limb_right;
103 : 360 : cur_limb_left++, cur_limb_right--) {
104 : 360 : mbedtls_mpi_uint tmp;
105 : : /* Note that if cur_limb_left == cur_limb_right,
106 : : * this code effectively swaps the bytes only once. */
107 : 360 : tmp = mpi_bigendian_to_host(*cur_limb_left);
108 : 360 : *cur_limb_left = mpi_bigendian_to_host(*cur_limb_right);
109 : 360 : *cur_limb_right = tmp;
110 : : }
111 : : }
112 : :
113 : : /* Whether min <= A, in constant time.
114 : : * A_limbs must be at least 1. */
115 : 12 : mbedtls_ct_condition_t mbedtls_mpi_core_uint_le_mpi(mbedtls_mpi_uint min,
116 : : const mbedtls_mpi_uint *A,
117 : : size_t A_limbs)
118 : : {
119 : : /* min <= least significant limb? */
120 : 12 : mbedtls_ct_condition_t min_le_lsl = mbedtls_ct_uint_ge(A[0], min);
121 : :
122 : : /* limbs other than the least significant one are all zero? */
123 : 12 : mbedtls_ct_condition_t msll_mask = MBEDTLS_CT_FALSE;
124 [ + + ]: 96 : for (size_t i = 1; i < A_limbs; i++) {
125 : 84 : msll_mask = mbedtls_ct_bool_or(msll_mask, mbedtls_ct_bool(A[i]));
126 : : }
127 : :
128 : : /* min <= A iff the lowest limb of A is >= min or the other limbs
129 : : * are not all zero. */
130 : 12 : return mbedtls_ct_bool_or(msll_mask, min_le_lsl);
131 : : }
132 : :
133 : 28684 : mbedtls_ct_condition_t mbedtls_mpi_core_lt_ct(const mbedtls_mpi_uint *A,
134 : : const mbedtls_mpi_uint *B,
135 : : size_t limbs)
136 : : {
137 : 28684 : mbedtls_ct_condition_t ret = MBEDTLS_CT_FALSE, cond = MBEDTLS_CT_FALSE, done = MBEDTLS_CT_FALSE;
138 : :
139 [ + + ]: 258156 : for (size_t i = limbs; i > 0; i--) {
140 : : /*
141 : : * If B[i - 1] < A[i - 1] then A < B is false and the result must
142 : : * remain 0.
143 : : *
144 : : * Again even if we can make a decision, we just mark the result and
145 : : * the fact that we are done and continue looping.
146 : : */
147 : 229472 : cond = mbedtls_ct_uint_lt(B[i - 1], A[i - 1]);
148 : 229472 : done = mbedtls_ct_bool_or(done, cond);
149 : :
150 : : /*
151 : : * If A[i - 1] < B[i - 1] then A < B is true.
152 : : *
153 : : * Again even if we can make a decision, we just mark the result and
154 : : * the fact that we are done and continue looping.
155 : : */
156 : 229472 : cond = mbedtls_ct_uint_lt(A[i - 1], B[i - 1]);
157 : 229472 : ret = mbedtls_ct_bool_or(ret, mbedtls_ct_bool_and(cond, mbedtls_ct_bool_not(done)));
158 : 229472 : done = mbedtls_ct_bool_or(done, cond);
159 : : }
160 : :
161 : : /*
162 : : * If all the limbs were equal, then the numbers are equal, A < B is false
163 : : * and leaving the result 0 is correct.
164 : : */
165 : :
166 : 28684 : return ret;
167 : : }
168 : :
169 : 200612 : void mbedtls_mpi_core_cond_assign(mbedtls_mpi_uint *X,
170 : : const mbedtls_mpi_uint *A,
171 : : size_t limbs,
172 : : mbedtls_ct_condition_t assign)
173 : : {
174 [ + - ]: 200612 : if (X == A) {
175 : : return;
176 : : }
177 : :
178 : : /* This function is very performance-sensitive for RSA. For this reason
179 : : * we have the loop below, instead of calling mbedtls_ct_memcpy_if
180 : : * (this is more optimal since here we don't have to handle the case where
181 : : * we copy awkwardly sized data).
182 : : */
183 [ + + ]: 1805508 : for (size_t i = 0; i < limbs; i++) {
184 : 1604896 : X[i] = mbedtls_ct_mpi_uint_if(assign, A[i], X[i]);
185 : : }
186 : : }
187 : :
188 : 57344 : void mbedtls_mpi_core_cond_swap(mbedtls_mpi_uint *X,
189 : : mbedtls_mpi_uint *Y,
190 : : size_t limbs,
191 : : mbedtls_ct_condition_t swap)
192 : : {
193 [ + - ]: 57344 : if (X == Y) {
194 : : return;
195 : : }
196 : :
197 [ + + ]: 516096 : for (size_t i = 0; i < limbs; i++) {
198 : 458752 : mbedtls_mpi_uint tmp = X[i];
199 : 458752 : X[i] = mbedtls_ct_mpi_uint_if(swap, Y[i], X[i]);
200 : 458752 : Y[i] = mbedtls_ct_mpi_uint_if(swap, tmp, Y[i]);
201 : : }
202 : : }
203 : :
204 : 0 : int mbedtls_mpi_core_read_le(mbedtls_mpi_uint *X,
205 : : size_t X_limbs,
206 : : const unsigned char *input,
207 : : size_t input_length)
208 : : {
209 : 0 : const size_t limbs = CHARS_TO_LIMBS(input_length);
210 : :
211 [ # # ]: 0 : if (X_limbs < limbs) {
212 : : return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
213 : : }
214 : :
215 [ # # ]: 0 : if (X != NULL) {
216 : 0 : memset(X, 0, X_limbs * ciL);
217 : :
218 [ # # ]: 0 : for (size_t i = 0; i < input_length; i++) {
219 : 0 : size_t offset = ((i % ciL) << 3);
220 : 0 : X[i / ciL] |= ((mbedtls_mpi_uint) input[i]) << offset;
221 : : }
222 : : }
223 : :
224 : : return 0;
225 : : }
226 : :
227 : 78 : int mbedtls_mpi_core_read_be(mbedtls_mpi_uint *X,
228 : : size_t X_limbs,
229 : : const unsigned char *input,
230 : : size_t input_length)
231 : : {
232 : 78 : const size_t limbs = CHARS_TO_LIMBS(input_length);
233 : :
234 [ + - ]: 78 : if (X_limbs < limbs) {
235 : : return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
236 : : }
237 : :
238 : : /* If X_limbs is 0, input_length must also be 0 (from previous test).
239 : : * Nothing to do. */
240 [ - + ]: 78 : if (X_limbs == 0) {
241 : : return 0;
242 : : }
243 : :
244 [ + - ]: 78 : memset(X, 0, X_limbs * ciL);
245 : :
246 : : /* memcpy() with (NULL, 0) is undefined behaviour */
247 [ + - ]: 78 : if (input_length != 0) {
248 : 78 : size_t overhead = (X_limbs * ciL) - input_length;
249 : 78 : unsigned char *Xp = (unsigned char *) X;
250 : 78 : memcpy(Xp + overhead, input, input_length);
251 : : }
252 : :
253 : 78 : mbedtls_mpi_core_bigendian_to_host(X, X_limbs);
254 : :
255 : 78 : return 0;
256 : : }
257 : :
258 : 0 : int mbedtls_mpi_core_write_le(const mbedtls_mpi_uint *A,
259 : : size_t A_limbs,
260 : : unsigned char *output,
261 : : size_t output_length)
262 : : {
263 : 0 : size_t stored_bytes = A_limbs * ciL;
264 : 0 : size_t bytes_to_copy;
265 : :
266 [ # # ]: 0 : if (stored_bytes < output_length) {
267 : : bytes_to_copy = stored_bytes;
268 : : } else {
269 : 0 : bytes_to_copy = output_length;
270 : :
271 : : /* The output buffer is smaller than the allocated size of A.
272 : : * However A may fit if its leading bytes are zero. */
273 [ # # ]: 0 : for (size_t i = bytes_to_copy; i < stored_bytes; i++) {
274 [ # # ]: 0 : if (GET_BYTE(A, i) != 0) {
275 : : return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
276 : : }
277 : : }
278 : : }
279 : :
280 [ # # ]: 0 : for (size_t i = 0; i < bytes_to_copy; i++) {
281 : 0 : output[i] = GET_BYTE(A, i);
282 : : }
283 : :
284 [ # # ]: 0 : if (stored_bytes < output_length) {
285 : : /* Write trailing 0 bytes */
286 : 0 : memset(output + stored_bytes, 0, output_length - stored_bytes);
287 : : }
288 : :
289 : : return 0;
290 : : }
291 : :
292 : 32 : int mbedtls_mpi_core_write_be(const mbedtls_mpi_uint *X,
293 : : size_t X_limbs,
294 : : unsigned char *output,
295 : : size_t output_length)
296 : : {
297 : 32 : size_t stored_bytes;
298 : 32 : size_t bytes_to_copy;
299 : 32 : unsigned char *p;
300 : :
301 : 32 : stored_bytes = X_limbs * ciL;
302 : :
303 [ - + ]: 32 : if (stored_bytes < output_length) {
304 : : /* There is enough space in the output buffer. Write initial
305 : : * null bytes and record the position at which to start
306 : : * writing the significant bytes. In this case, the execution
307 : : * trace of this function does not depend on the value of the
308 : : * number. */
309 : 0 : bytes_to_copy = stored_bytes;
310 : 0 : p = output + output_length - stored_bytes;
311 : 0 : memset(output, 0, output_length - stored_bytes);
312 : : } else {
313 : : /* The output buffer is smaller than the allocated size of X.
314 : : * However X may fit if its leading bytes are zero. */
315 : 288 : bytes_to_copy = output_length;
316 : 288 : p = output;
317 [ + + ]: 288 : for (size_t i = bytes_to_copy; i < stored_bytes; i++) {
318 [ + - ]: 256 : if (GET_BYTE(X, i) != 0) {
319 : : return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
320 : : }
321 : : }
322 : : }
323 : :
324 [ + + ]: 1056 : for (size_t i = 0; i < bytes_to_copy; i++) {
325 : 1024 : p[bytes_to_copy - i - 1] = GET_BYTE(X, i);
326 : : }
327 : :
328 : : return 0;
329 : : }
330 : :
331 : 135204 : void mbedtls_mpi_core_shift_r(mbedtls_mpi_uint *X, size_t limbs,
332 : : size_t count)
333 : : {
334 : 135204 : size_t i, v0, v1;
335 : 135204 : mbedtls_mpi_uint r0 = 0, r1;
336 : :
337 : 135204 : v0 = count / biL;
338 : 135204 : v1 = count & (biL - 1);
339 : :
340 [ + - - + ]: 135204 : if (v0 > limbs || (v0 == limbs && v1 > 0)) {
341 : 0 : memset(X, 0, limbs * ciL);
342 : 0 : return;
343 : : }
344 : :
345 : : /*
346 : : * shift by count / limb_size
347 : : */
348 [ + + ]: 135204 : if (v0 > 0) {
349 [ + + ]: 389240 : for (i = 0; i < limbs - v0; i++) {
350 : 350316 : X[i] = X[i + v0];
351 : : }
352 : :
353 [ + + ]: 350348 : for (; i < limbs; i++) {
354 : 311424 : X[i] = 0;
355 : : }
356 : : }
357 : :
358 : : /*
359 : : * shift by count % limb_size
360 : : */
361 [ + + ]: 135204 : if (v1 > 0) {
362 [ + + ]: 1216760 : for (i = limbs; i > 0; i--) {
363 : 1120492 : r1 = X[i - 1] << (biL - v1);
364 : 1120492 : X[i - 1] >>= v1;
365 : 1120492 : X[i - 1] |= r0;
366 : 1120492 : r0 = r1;
367 : : }
368 : : }
369 : : }
370 : :
371 : 441752 : void mbedtls_mpi_core_shift_l(mbedtls_mpi_uint *X, size_t limbs,
372 : : size_t count)
373 : : {
374 : 441752 : size_t i, v0, v1;
375 : 441752 : mbedtls_mpi_uint r0 = 0, r1;
376 : :
377 : 441752 : v0 = count / (biL);
378 : 441752 : v1 = count & (biL - 1);
379 : :
380 : : /*
381 : : * shift by count / limb_size
382 : : */
383 [ + + ]: 441752 : if (v0 > 0) {
384 [ + + ]: 4476980 : for (i = limbs; i > v0; i--) {
385 : 4165556 : X[i - 1] = X[i - v0 - 1];
386 : : }
387 : :
388 [ + + ]: 1713088 : for (; i > 0; i--) {
389 : 1401664 : X[i - 1] = 0;
390 : : }
391 : : }
392 : :
393 : : /*
394 : : * shift by count % limb_size
395 : : */
396 [ + + ]: 441752 : if (v1 > 0) {
397 [ + + ]: 1320272 : for (i = v0; i < limbs; i++) {
398 : 1228868 : r1 = X[i] >> (biL - v1);
399 : 1228868 : X[i] <<= v1;
400 : 1228868 : X[i] |= r0;
401 : 1228868 : r0 = r1;
402 : : }
403 : : }
404 : 441752 : }
405 : :
406 : 3124 : mbedtls_mpi_uint mbedtls_mpi_core_add(mbedtls_mpi_uint *X,
407 : : const mbedtls_mpi_uint *A,
408 : : const mbedtls_mpi_uint *B,
409 : : size_t limbs)
410 : : {
411 : 3124 : mbedtls_mpi_uint c = 0;
412 : :
413 [ + + ]: 27994 : for (size_t i = 0; i < limbs; i++) {
414 : 24870 : mbedtls_mpi_uint t = c + A[i];
415 : 24870 : c = (t < A[i]);
416 : 24870 : t += B[i];
417 : 24870 : c += (t < B[i]);
418 : 24870 : X[i] = t;
419 : : }
420 : :
421 : 3124 : return c;
422 : : }
423 : :
424 : 57344 : mbedtls_mpi_uint mbedtls_mpi_core_add_if(mbedtls_mpi_uint *X,
425 : : const mbedtls_mpi_uint *A,
426 : : size_t limbs,
427 : : unsigned cond)
428 : : {
429 : 57344 : mbedtls_mpi_uint c = 0;
430 : :
431 : 57344 : mbedtls_ct_condition_t do_add = mbedtls_ct_bool(cond);
432 : :
433 [ + + ]: 516096 : for (size_t i = 0; i < limbs; i++) {
434 : 458752 : mbedtls_mpi_uint add = mbedtls_ct_mpi_uint_if_else_0(do_add, A[i]);
435 : 458752 : mbedtls_mpi_uint t = c + X[i];
436 : 458752 : c = (t < X[i]);
437 : 458752 : t += add;
438 : 458752 : c += (t < add);
439 : 458752 : X[i] = t;
440 : : }
441 : :
442 : 57344 : return c;
443 : : }
444 : :
445 : 415617 : mbedtls_mpi_uint mbedtls_mpi_core_sub(mbedtls_mpi_uint *X,
446 : : const mbedtls_mpi_uint *A,
447 : : const mbedtls_mpi_uint *B,
448 : : size_t limbs)
449 : : {
450 : 415617 : mbedtls_mpi_uint c = 0;
451 : :
452 [ + + ]: 5453206 : for (size_t i = 0; i < limbs; i++) {
453 : 5037589 : mbedtls_mpi_uint z = (A[i] < c);
454 : 5037589 : mbedtls_mpi_uint t = A[i] - c;
455 : 5037589 : c = (t < B[i]) + z;
456 : 5037589 : X[i] = t - B[i];
457 : : }
458 : :
459 : 415617 : return c;
460 : : }
461 : :
462 : 1085130 : mbedtls_mpi_uint mbedtls_mpi_core_mla(mbedtls_mpi_uint *d, size_t d_len,
463 : : const mbedtls_mpi_uint *s, size_t s_len,
464 : : mbedtls_mpi_uint b)
465 : : {
466 : 1085130 : mbedtls_mpi_uint c = 0; /* carry */
467 : : /*
468 : : * It is a documented precondition of this function that d_len >= s_len.
469 : : * If that's not the case, we swap these round: this turns what would be
470 : : * a buffer overflow into an incorrect result.
471 : : */
472 : 1085130 : if (d_len < s_len) {
473 : : s_len = d_len;
474 : : }
475 : 1085130 : size_t excess_len = d_len - s_len;
476 : 1085130 : size_t steps_x8 = s_len / 8;
477 : 1085130 : size_t steps_x1 = s_len & 7;
478 : :
479 [ + + ]: 1711112 : while (steps_x8--) {
480 : 625982 : MULADDC_X8_INIT
481 : 625982 : MULADDC_X8_CORE
482 : : MULADDC_X8_STOP
483 : : }
484 : :
485 [ + + ]: 2313828 : while (steps_x1--) {
486 : 1228698 : MULADDC_X1_INIT
487 : 1228698 : MULADDC_X1_CORE
488 : : MULADDC_X1_STOP
489 : : }
490 : :
491 [ + + ]: 11542746 : while (excess_len--) {
492 : 10457616 : *d += c;
493 : 10457616 : c = (*d < c);
494 : 10457616 : d++;
495 : : }
496 : :
497 : 1085130 : return c;
498 : : }
499 : :
500 : 39160 : void mbedtls_mpi_core_mul(mbedtls_mpi_uint *X,
501 : : const mbedtls_mpi_uint *A, size_t A_limbs,
502 : : const mbedtls_mpi_uint *B, size_t B_limbs)
503 : : {
504 : 39160 : memset(X, 0, (A_limbs + B_limbs) * ciL);
505 : :
506 [ + + ]: 351660 : for (size_t i = 0; i < B_limbs; i++) {
507 : 312500 : (void) mbedtls_mpi_core_mla(X + i, A_limbs + 1, A, A_limbs, B[i]);
508 : : }
509 : 39160 : }
510 : :
511 : : /*
512 : : * Fast Montgomery initialization (thanks to Tom St Denis).
513 : : */
514 : 0 : mbedtls_mpi_uint mbedtls_mpi_core_montmul_init(const mbedtls_mpi_uint *N)
515 : : {
516 : 0 : mbedtls_mpi_uint x = N[0];
517 : :
518 : 0 : x += ((N[0] + 2) & 4) << 1;
519 : :
520 [ # # ]: 0 : for (unsigned int i = biL; i >= 8; i /= 2) {
521 : 0 : x *= (2 - (N[0] * x));
522 : : }
523 : :
524 : 0 : return ~x + 1;
525 : : }
526 : :
527 : 0 : void mbedtls_mpi_core_montmul(mbedtls_mpi_uint *X,
528 : : const mbedtls_mpi_uint *A,
529 : : const mbedtls_mpi_uint *B,
530 : : size_t B_limbs,
531 : : const mbedtls_mpi_uint *N,
532 : : size_t AN_limbs,
533 : : mbedtls_mpi_uint mm,
534 : : mbedtls_mpi_uint *T)
535 : : {
536 : 0 : memset(T, 0, (2 * AN_limbs + 1) * ciL);
537 : :
538 [ # # ]: 0 : for (size_t i = 0; i < AN_limbs; i++) {
539 : : /* T = (T + u0*B + u1*N) / 2^biL */
540 : 0 : mbedtls_mpi_uint u0 = A[i];
541 : 0 : mbedtls_mpi_uint u1 = (T[0] + u0 * B[0]) * mm;
542 : :
543 : 0 : (void) mbedtls_mpi_core_mla(T, AN_limbs + 2, B, B_limbs, u0);
544 : 0 : (void) mbedtls_mpi_core_mla(T, AN_limbs + 2, N, AN_limbs, u1);
545 : :
546 : 0 : T++;
547 : : }
548 : :
549 : : /*
550 : : * The result we want is (T >= N) ? T - N : T.
551 : : *
552 : : * For better constant-time properties in this function, we always do the
553 : : * subtraction, with the result in X.
554 : : *
555 : : * We also look to see if there was any carry in the final additions in the
556 : : * loop above.
557 : : */
558 : :
559 : 0 : mbedtls_mpi_uint carry = T[AN_limbs];
560 : 0 : mbedtls_mpi_uint borrow = mbedtls_mpi_core_sub(X, T, N, AN_limbs);
561 : :
562 : : /*
563 : : * Using R as the Montgomery radix (auxiliary modulus) i.e. 2^(biL*AN_limbs):
564 : : *
565 : : * T can be in one of 3 ranges:
566 : : *
567 : : * 1) T < N : (carry, borrow) = (0, 1): we want T
568 : : * 2) N <= T < R : (carry, borrow) = (0, 0): we want X
569 : : * 3) T >= R : (carry, borrow) = (1, 1): we want X
570 : : *
571 : : * and (carry, borrow) = (1, 0) can't happen.
572 : : *
573 : : * So the correct return value is already in X if (carry ^ borrow) = 0,
574 : : * but is in (the lower AN_limbs limbs of) T if (carry ^ borrow) = 1.
575 : : */
576 : 0 : mbedtls_ct_memcpy_if(mbedtls_ct_bool(carry ^ borrow),
577 : : (unsigned char *) X,
578 : : (unsigned char *) T,
579 : : NULL,
580 : : AN_limbs * sizeof(mbedtls_mpi_uint));
581 : 0 : }
582 : :
583 : 0 : int mbedtls_mpi_core_get_mont_r2_unsafe(mbedtls_mpi *X,
584 : : const mbedtls_mpi *N)
585 : : {
586 : 0 : int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
587 : :
588 [ # # ]: 0 : MBEDTLS_MPI_CHK(mbedtls_mpi_lset(X, 1));
589 [ # # ]: 0 : MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(X, N->n * 2 * biL));
590 [ # # ]: 0 : MBEDTLS_MPI_CHK(mbedtls_mpi_mod_mpi(X, X, N));
591 : 0 : MBEDTLS_MPI_CHK(mbedtls_mpi_shrink(X, N->n));
592 : :
593 : 0 : cleanup:
594 : 0 : return ret;
595 : : }
596 : :
597 : : MBEDTLS_STATIC_TESTABLE
598 : 0 : void mbedtls_mpi_core_ct_uint_table_lookup(mbedtls_mpi_uint *dest,
599 : : const mbedtls_mpi_uint *table,
600 : : size_t limbs,
601 : : size_t count,
602 : : size_t index)
603 : : {
604 [ # # ]: 0 : for (size_t i = 0; i < count; i++, table += limbs) {
605 : 0 : mbedtls_ct_condition_t assign = mbedtls_ct_uint_eq(i, index);
606 : 0 : mbedtls_mpi_core_cond_assign(dest, table, limbs, assign);
607 : : }
608 : 0 : }
609 : :
610 : : /* Fill X with n_bytes random bytes.
611 : : * X must already have room for those bytes.
612 : : * The ordering of the bytes returned from the RNG is suitable for
613 : : * deterministic ECDSA (see RFC 6979 §3.3 and the specification of
614 : : * mbedtls_mpi_core_random()).
615 : : */
616 : 12 : int mbedtls_mpi_core_fill_random(
617 : : mbedtls_mpi_uint *X, size_t X_limbs,
618 : : size_t n_bytes,
619 : : int (*f_rng)(void *, unsigned char *, size_t), void *p_rng)
620 : : {
621 : 12 : int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
622 : 12 : const size_t limbs = CHARS_TO_LIMBS(n_bytes);
623 : 12 : const size_t overhead = (limbs * ciL) - n_bytes;
624 : :
625 [ + - ]: 12 : if (X_limbs < limbs) {
626 : : return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
627 : : }
628 : :
629 : 12 : memset(X, 0, overhead);
630 : 12 : memset((unsigned char *) X + limbs * ciL, 0, (X_limbs - limbs) * ciL);
631 [ - + ]: 12 : MBEDTLS_MPI_CHK(f_rng(p_rng, (unsigned char *) X + overhead, n_bytes));
632 : 12 : mbedtls_mpi_core_bigendian_to_host(X, limbs);
633 : :
634 : : cleanup:
635 : : return ret;
636 : : }
637 : :
638 : 12 : int mbedtls_mpi_core_random(mbedtls_mpi_uint *X,
639 : : mbedtls_mpi_uint min,
640 : : const mbedtls_mpi_uint *N,
641 : : size_t limbs,
642 : : int (*f_rng)(void *, unsigned char *, size_t),
643 : : void *p_rng)
644 : : {
645 : 12 : mbedtls_ct_condition_t ge_lower = MBEDTLS_CT_TRUE, lt_upper = MBEDTLS_CT_FALSE;
646 : 12 : size_t n_bits = mbedtls_mpi_core_bitlen(N, limbs);
647 : 12 : size_t n_bytes = (n_bits + 7) / 8;
648 : 12 : int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
649 : :
650 : : /*
651 : : * When min == 0, each try has at worst a probability 1/2 of failing
652 : : * (the msb has a probability 1/2 of being 0, and then the result will
653 : : * be < N), so after 30 tries failure probability is a most 2**(-30).
654 : : *
655 : : * When N is just below a power of 2, as is the case when generating
656 : : * a random scalar on most elliptic curves, 1 try is enough with
657 : : * overwhelming probability. When N is just above a power of 2,
658 : : * as when generating a random scalar on secp224k1, each try has
659 : : * a probability of failing that is almost 1/2.
660 : : *
661 : : * The probabilities are almost the same if min is nonzero but negligible
662 : : * compared to N. This is always the case when N is crypto-sized, but
663 : : * it's convenient to support small N for testing purposes. When N
664 : : * is small, use a higher repeat count, otherwise the probability of
665 : : * failure is macroscopic.
666 : : */
667 [ - + ]: 12 : int count = (n_bytes > 4 ? 30 : 250);
668 : :
669 : : /*
670 : : * Match the procedure given in RFC 6979 §3.3 (deterministic ECDSA)
671 : : * when f_rng is a suitably parametrized instance of HMAC_DRBG:
672 : : * - use the same byte ordering;
673 : : * - keep the leftmost n_bits bits of the generated octet string;
674 : : * - try until result is in the desired range.
675 : : * This also avoids any bias, which is especially important for ECDSA.
676 : : */
677 : 12 : do {
678 [ - + ]: 12 : MBEDTLS_MPI_CHK(mbedtls_mpi_core_fill_random(X, limbs,
679 : : n_bytes,
680 : : f_rng, p_rng));
681 : 12 : mbedtls_mpi_core_shift_r(X, limbs, 8 * n_bytes - n_bits);
682 : :
683 [ - + ]: 12 : if (--count == 0) {
684 : 0 : ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
685 : 0 : goto cleanup;
686 : : }
687 : :
688 : 12 : ge_lower = mbedtls_mpi_core_uint_le_mpi(min, X, limbs);
689 : 12 : lt_upper = mbedtls_mpi_core_lt_ct(X, N, limbs);
690 [ - + ]: 12 : } while (mbedtls_ct_bool_and(ge_lower, lt_upper) == MBEDTLS_CT_FALSE);
691 : :
692 : 12 : cleanup:
693 : 12 : return ret;
694 : : }
695 : :
696 : 0 : static size_t exp_mod_get_window_size(size_t Ebits)
697 : : {
698 : : #if MBEDTLS_MPI_WINDOW_SIZE >= 6
699 : : return (Ebits > 671) ? 6 : (Ebits > 239) ? 5 : (Ebits > 79) ? 4 : 1;
700 : : #elif MBEDTLS_MPI_WINDOW_SIZE == 5
701 : : return (Ebits > 239) ? 5 : (Ebits > 79) ? 4 : 1;
702 : : #elif MBEDTLS_MPI_WINDOW_SIZE > 1
703 [ # # ]: 0 : return (Ebits > 79) ? MBEDTLS_MPI_WINDOW_SIZE : 1;
704 : : #else
705 : : (void) Ebits;
706 : : return 1;
707 : : #endif
708 : : }
709 : :
710 : 0 : size_t mbedtls_mpi_core_exp_mod_working_limbs(size_t AN_limbs, size_t E_limbs)
711 : : {
712 : 0 : const size_t wsize = exp_mod_get_window_size(E_limbs * biL);
713 : 0 : const size_t welem = ((size_t) 1) << wsize;
714 : :
715 : : /* How big does each part of the working memory pool need to be? */
716 : 0 : const size_t table_limbs = welem * AN_limbs;
717 : 0 : const size_t select_limbs = AN_limbs;
718 : 0 : const size_t temp_limbs = 2 * AN_limbs + 1;
719 : :
720 : 0 : return table_limbs + select_limbs + temp_limbs;
721 : : }
722 : :
723 : 0 : static void exp_mod_precompute_window(const mbedtls_mpi_uint *A,
724 : : const mbedtls_mpi_uint *N,
725 : : size_t AN_limbs,
726 : : mbedtls_mpi_uint mm,
727 : : const mbedtls_mpi_uint *RR,
728 : : size_t welem,
729 : : mbedtls_mpi_uint *Wtable,
730 : : mbedtls_mpi_uint *temp)
731 : : {
732 : : /* W[0] = 1 (in Montgomery presentation) */
733 : 0 : memset(Wtable, 0, AN_limbs * ciL);
734 : 0 : Wtable[0] = 1;
735 : 0 : mbedtls_mpi_core_montmul(Wtable, Wtable, RR, AN_limbs, N, AN_limbs, mm, temp);
736 : :
737 : : /* W[1] = A (already in Montgomery presentation) */
738 : 0 : mbedtls_mpi_uint *W1 = Wtable + AN_limbs;
739 : 0 : memcpy(W1, A, AN_limbs * ciL);
740 : :
741 : : /* W[i+1] = W[i] * W[1], i >= 2 */
742 : 0 : mbedtls_mpi_uint *Wprev = W1;
743 [ # # ]: 0 : for (size_t i = 2; i < welem; i++) {
744 : 0 : mbedtls_mpi_uint *Wcur = Wprev + AN_limbs;
745 : 0 : mbedtls_mpi_core_montmul(Wcur, Wprev, W1, AN_limbs, N, AN_limbs, mm, temp);
746 : 0 : Wprev = Wcur;
747 : : }
748 : 0 : }
749 : :
750 : : #if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
751 : : void (*mbedtls_safe_codepath_hook)(void) = NULL;
752 : : void (*mbedtls_unsafe_codepath_hook)(void) = NULL;
753 : : #endif
754 : :
755 : : /*
756 : : * This function calculates the indices of the exponent where the exponentiation algorithm should
757 : : * start processing.
758 : : *
759 : : * Warning! If the parameter E_public has MBEDTLS_MPI_IS_PUBLIC as its value,
760 : : * this function is not constant time with respect to the exponent (parameter E).
761 : : */
762 : 0 : static inline void exp_mod_calc_first_bit_optionally_safe(const mbedtls_mpi_uint *E,
763 : : size_t E_limbs,
764 : : int E_public,
765 : : size_t *E_limb_index,
766 : : size_t *E_bit_index)
767 : : {
768 [ # # ]: 0 : if (E_public == MBEDTLS_MPI_IS_PUBLIC) {
769 : : /*
770 : : * Skip leading zero bits.
771 : : */
772 : 0 : size_t E_bits = mbedtls_mpi_core_bitlen(E, E_limbs);
773 : 0 : if (E_bits == 0) {
774 : : /*
775 : : * If E is 0 mbedtls_mpi_core_bitlen() returns 0. Even if that is the case, we will want
776 : : * to represent it as a single 0 bit and as such the bitlength will be 1.
777 : : */
778 : : E_bits = 1;
779 : : }
780 : :
781 : 0 : *E_limb_index = E_bits / biL;
782 : 0 : *E_bit_index = E_bits % biL;
783 : :
784 : : #if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
785 : : if (mbedtls_unsafe_codepath_hook != NULL) {
786 : : mbedtls_unsafe_codepath_hook();
787 : : }
788 : : #endif
789 : : } else {
790 : : /*
791 : : * Here we need to be constant time with respect to E and can't do anything better than
792 : : * start at the first allocated bit.
793 : : */
794 : 0 : *E_limb_index = E_limbs;
795 : 0 : *E_bit_index = 0;
796 : : #if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
797 : : if (mbedtls_safe_codepath_hook != NULL) {
798 : : mbedtls_safe_codepath_hook();
799 : : }
800 : : #endif
801 : : }
802 : 0 : }
803 : :
804 : : /*
805 : : * Warning! If the parameter window_public has MBEDTLS_MPI_IS_PUBLIC as its value, this function is
806 : : * not constant time with respect to the window parameter and consequently the exponent of the
807 : : * exponentiation (parameter E of mbedtls_mpi_core_exp_mod_optionally_safe).
808 : : */
809 : 0 : static inline void exp_mod_table_lookup_optionally_safe(mbedtls_mpi_uint *Wselect,
810 : : mbedtls_mpi_uint *Wtable,
811 : : size_t AN_limbs, size_t welem,
812 : : mbedtls_mpi_uint window,
813 : : int window_public)
814 : : {
815 [ # # ]: 0 : if (window_public == MBEDTLS_MPI_IS_PUBLIC) {
816 : 0 : memcpy(Wselect, Wtable + window * AN_limbs, AN_limbs * ciL);
817 : : #if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
818 : : if (mbedtls_unsafe_codepath_hook != NULL) {
819 : : mbedtls_unsafe_codepath_hook();
820 : : }
821 : : #endif
822 : : } else {
823 : : /* Select Wtable[window] without leaking window through
824 : : * memory access patterns. */
825 : 0 : mbedtls_mpi_core_ct_uint_table_lookup(Wselect, Wtable,
826 : : AN_limbs, welem, window);
827 : : #if defined(MBEDTLS_TEST_HOOKS) && !defined(MBEDTLS_THREADING_C)
828 : : if (mbedtls_safe_codepath_hook != NULL) {
829 : : mbedtls_safe_codepath_hook();
830 : : }
831 : : #endif
832 : : }
833 : 0 : }
834 : :
835 : : /* Exponentiation: X := A^E mod N.
836 : : *
837 : : * Warning! If the parameter E_public has MBEDTLS_MPI_IS_PUBLIC as its value,
838 : : * this function is not constant time with respect to the exponent (parameter E).
839 : : *
840 : : * A must already be in Montgomery form.
841 : : *
842 : : * As in other bignum functions, assume that AN_limbs and E_limbs are nonzero.
843 : : *
844 : : * RR must contain 2^{2*biL} mod N.
845 : : *
846 : : * The algorithm is a variant of Left-to-right k-ary exponentiation: HAC 14.82
847 : : * (The difference is that the body in our loop processes a single bit instead
848 : : * of a full window.)
849 : : */
850 : 0 : static void mbedtls_mpi_core_exp_mod_optionally_safe(mbedtls_mpi_uint *X,
851 : : const mbedtls_mpi_uint *A,
852 : : const mbedtls_mpi_uint *N,
853 : : size_t AN_limbs,
854 : : const mbedtls_mpi_uint *E,
855 : : size_t E_limbs,
856 : : int E_public,
857 : : const mbedtls_mpi_uint *RR,
858 : : mbedtls_mpi_uint *T)
859 : : {
860 : : /* We'll process the bits of E from most significant
861 : : * (limb_index=E_limbs-1, E_bit_index=biL-1) to least significant
862 : : * (limb_index=0, E_bit_index=0). */
863 : 0 : size_t E_limb_index = E_limbs;
864 : 0 : size_t E_bit_index = 0;
865 : 0 : exp_mod_calc_first_bit_optionally_safe(E, E_limbs, E_public,
866 : : &E_limb_index, &E_bit_index);
867 : :
868 : 0 : const size_t wsize = exp_mod_get_window_size(E_limb_index * biL);
869 : 0 : const size_t welem = ((size_t) 1) << wsize;
870 : :
871 : : /* This is how we will use the temporary storage T, which must have space
872 : : * for table_limbs, select_limbs and (2 * AN_limbs + 1) for montmul. */
873 : 0 : const size_t table_limbs = welem * AN_limbs;
874 : 0 : const size_t select_limbs = AN_limbs;
875 : :
876 : : /* Pointers to specific parts of the temporary working memory pool */
877 : 0 : mbedtls_mpi_uint *const Wtable = T;
878 : 0 : mbedtls_mpi_uint *const Wselect = Wtable + table_limbs;
879 : 0 : mbedtls_mpi_uint *const temp = Wselect + select_limbs;
880 : :
881 : : /*
882 : : * Window precomputation
883 : : */
884 : :
885 : 0 : const mbedtls_mpi_uint mm = mbedtls_mpi_core_montmul_init(N);
886 : :
887 : : /* Set Wtable[i] = A^i (in Montgomery representation) */
888 : 0 : exp_mod_precompute_window(A, N, AN_limbs,
889 : : mm, RR,
890 : : welem, Wtable, temp);
891 : :
892 : : /*
893 : : * Fixed window exponentiation
894 : : */
895 : :
896 : : /* X = 1 (in Montgomery presentation) initially */
897 : 0 : memcpy(X, Wtable, AN_limbs * ciL);
898 : :
899 : : /* At any given time, window contains window_bits bits from E.
900 : : * window_bits can go up to wsize. */
901 : 0 : size_t window_bits = 0;
902 : 0 : mbedtls_mpi_uint window = 0;
903 : :
904 : 0 : do {
905 : : /* Square */
906 : 0 : mbedtls_mpi_core_montmul(X, X, X, AN_limbs, N, AN_limbs, mm, temp);
907 : :
908 : : /* Move to the next bit of the exponent */
909 [ # # ]: 0 : if (E_bit_index == 0) {
910 : 0 : --E_limb_index;
911 : 0 : E_bit_index = biL - 1;
912 : : } else {
913 : 0 : --E_bit_index;
914 : : }
915 : : /* Insert next exponent bit into window */
916 : 0 : ++window_bits;
917 : 0 : window <<= 1;
918 : 0 : window |= (E[E_limb_index] >> E_bit_index) & 1;
919 : :
920 : : /* Clear window if it's full. Also clear the window at the end,
921 : : * when we've finished processing the exponent. */
922 [ # # # # ]: 0 : if (window_bits == wsize ||
923 [ # # ]: 0 : (E_bit_index == 0 && E_limb_index == 0)) {
924 : :
925 : 0 : exp_mod_table_lookup_optionally_safe(Wselect, Wtable, AN_limbs, welem,
926 : : window, E_public);
927 : : /* Multiply X by the selected element. */
928 : 0 : mbedtls_mpi_core_montmul(X, X, Wselect, AN_limbs, N, AN_limbs, mm,
929 : : temp);
930 : 0 : window = 0;
931 : 0 : window_bits = 0;
932 : : }
933 [ # # # # ]: 0 : } while (!(E_bit_index == 0 && E_limb_index == 0));
934 : 0 : }
935 : :
936 : 0 : void mbedtls_mpi_core_exp_mod(mbedtls_mpi_uint *X,
937 : : const mbedtls_mpi_uint *A,
938 : : const mbedtls_mpi_uint *N, size_t AN_limbs,
939 : : const mbedtls_mpi_uint *E, size_t E_limbs,
940 : : const mbedtls_mpi_uint *RR,
941 : : mbedtls_mpi_uint *T)
942 : : {
943 : 0 : mbedtls_mpi_core_exp_mod_optionally_safe(X,
944 : : A,
945 : : N,
946 : : AN_limbs,
947 : : E,
948 : : E_limbs,
949 : : MBEDTLS_MPI_IS_SECRET,
950 : : RR,
951 : : T);
952 : 0 : }
953 : :
954 : 0 : void mbedtls_mpi_core_exp_mod_unsafe(mbedtls_mpi_uint *X,
955 : : const mbedtls_mpi_uint *A,
956 : : const mbedtls_mpi_uint *N, size_t AN_limbs,
957 : : const mbedtls_mpi_uint *E, size_t E_limbs,
958 : : const mbedtls_mpi_uint *RR,
959 : : mbedtls_mpi_uint *T)
960 : : {
961 : 0 : mbedtls_mpi_core_exp_mod_optionally_safe(X,
962 : : A,
963 : : N,
964 : : AN_limbs,
965 : : E,
966 : : E_limbs,
967 : : MBEDTLS_MPI_IS_PUBLIC,
968 : : RR,
969 : : T);
970 : 0 : }
971 : :
972 : 11384 : mbedtls_mpi_uint mbedtls_mpi_core_sub_int(mbedtls_mpi_uint *X,
973 : : const mbedtls_mpi_uint *A,
974 : : mbedtls_mpi_uint c, /* doubles as carry */
975 : : size_t limbs)
976 : : {
977 [ + + ]: 88932 : for (size_t i = 0; i < limbs; i++) {
978 : 77548 : mbedtls_mpi_uint s = A[i];
979 : 77548 : mbedtls_mpi_uint t = s - c;
980 : 77548 : c = (t > s);
981 : 77548 : X[i] = t;
982 : : }
983 : :
984 : 11384 : return c;
985 : : }
986 : :
987 : 0 : mbedtls_ct_condition_t mbedtls_mpi_core_check_zero_ct(const mbedtls_mpi_uint *A,
988 : : size_t limbs)
989 : : {
990 : 0 : volatile const mbedtls_mpi_uint *force_read_A = A;
991 : 0 : mbedtls_mpi_uint bits = 0;
992 : :
993 [ # # ]: 0 : for (size_t i = 0; i < limbs; i++) {
994 : 0 : bits |= force_read_A[i];
995 : : }
996 : :
997 : 0 : return mbedtls_ct_bool(bits);
998 : : }
999 : :
1000 : 0 : void mbedtls_mpi_core_to_mont_rep(mbedtls_mpi_uint *X,
1001 : : const mbedtls_mpi_uint *A,
1002 : : const mbedtls_mpi_uint *N,
1003 : : size_t AN_limbs,
1004 : : mbedtls_mpi_uint mm,
1005 : : const mbedtls_mpi_uint *rr,
1006 : : mbedtls_mpi_uint *T)
1007 : : {
1008 : 0 : mbedtls_mpi_core_montmul(X, A, rr, AN_limbs, N, AN_limbs, mm, T);
1009 : 0 : }
1010 : :
1011 : 0 : void mbedtls_mpi_core_from_mont_rep(mbedtls_mpi_uint *X,
1012 : : const mbedtls_mpi_uint *A,
1013 : : const mbedtls_mpi_uint *N,
1014 : : size_t AN_limbs,
1015 : : mbedtls_mpi_uint mm,
1016 : : mbedtls_mpi_uint *T)
1017 : : {
1018 : 0 : const mbedtls_mpi_uint Rinv = 1; /* 1/R in Mont. rep => 1 */
1019 : :
1020 : 0 : mbedtls_mpi_core_montmul(X, A, &Rinv, 1, N, AN_limbs, mm, T);
1021 : 0 : }
1022 : :
1023 : : /*
1024 : : * Compute X = A - B mod N.
1025 : : * Both A and B must be in [0, N) and so will the output.
1026 : : */
1027 : 28672 : static void mpi_core_sub_mod(mbedtls_mpi_uint *X,
1028 : : const mbedtls_mpi_uint *A,
1029 : : const mbedtls_mpi_uint *B,
1030 : : const mbedtls_mpi_uint *N,
1031 : : size_t limbs)
1032 : : {
1033 : 28672 : mbedtls_mpi_uint c = mbedtls_mpi_core_sub(X, A, B, limbs);
1034 : 28672 : (void) mbedtls_mpi_core_add_if(X, N, limbs, (unsigned) c);
1035 : 28672 : }
1036 : :
1037 : : /*
1038 : : * Divide X by 2 mod N in place, assuming N is odd.
1039 : : * The input must be in [0, N) and so will the output.
1040 : : */
1041 : : MBEDTLS_STATIC_TESTABLE
1042 : 28672 : void mbedtls_mpi_core_div2_mod_odd(mbedtls_mpi_uint *X,
1043 : : const mbedtls_mpi_uint *N,
1044 : : size_t limbs)
1045 : : {
1046 : : /* If X is odd, add N to make it even before shifting. */
1047 : 28672 : unsigned odd = (unsigned) X[0] & 1;
1048 : 28672 : mbedtls_mpi_uint c = mbedtls_mpi_core_add_if(X, N, limbs, odd);
1049 : 28672 : mbedtls_mpi_core_shift_r(X, limbs, 1);
1050 : 28672 : X[limbs - 1] |= c << (biL - 1);
1051 : 28672 : }
1052 : :
1053 : : /*
1054 : : * Constant-time GCD and modular inversion - odd modulus.
1055 : : *
1056 : : * Pre-conditions: see public documentation.
1057 : : *
1058 : : * See https://www.jstage.jst.go.jp/article/transinf/E106.D/9/E106.D_2022ICP0009/_pdf
1059 : : *
1060 : : * The paper gives two computationally equivalent algorithms: Alg 7 (readable)
1061 : : * and Alg 8 (constant-time). We use a third version that's hopefully both:
1062 : : *
1063 : : * u, v = A, N # N is called p in the paper but doesn't have to be prime
1064 : : * q, r = 0, 1
1065 : : * repeat bits(A_limbs + N_limbs) times:
1066 : : * d = v - u # t1 in Alg 7
1067 : : * t1 = (u and v both odd) ? u : d # t1 in Alg 8
1068 : : * t2 = (u and v both odd) ? d : (u odd) ? v : u # t2 in Alg 8
1069 : : * t2 >>= 1
1070 : : * swap = t1 > t2 # similar to s, z in Alg 8
1071 : : * u, v = (swap) ? t2, t1 : t1, t2
1072 : : *
1073 : : * d = r - q mod N # t2 in Alg 7
1074 : : * t1 = (u and v both odd) ? q : d # t3 in Alg 8
1075 : : * t2 = (u and v both odd) ? d : (u odd) ? r : q # t4 Alg 8
1076 : : * t2 /= 2 mod N # see below (pre_com)
1077 : : * q, r = (swap) ? t2, t1 : t1, t2
1078 : : * return v, q # v: GCD, see Alg 6; q: no mult by pre_com, see below
1079 : : *
1080 : : * The ternary operators in the above pseudo-code need to be realised in a
1081 : : * constant-time fashion. We use conditional assign for t1, t2 and conditional
1082 : : * swap for the final update. (Note: the similarity between branches of Alg 7
1083 : : * are highlighted in tables 2 and 3 and the surrounding text.)
1084 : : *
1085 : : * Also, we re-order operations, grouping things related to the inverse, which
1086 : : * facilitates making its computation optional, and requires fewer temporaries.
1087 : : *
1088 : : * The only actual change from the paper is dropping the trick with pre_com,
1089 : : * which I think complicates things for no benefit.
1090 : : * See the comment on the big I != NULL block below for details.
1091 : : */
1092 : 56 : void mbedtls_mpi_core_gcd_modinv_odd(mbedtls_mpi_uint *G,
1093 : : mbedtls_mpi_uint *I,
1094 : : const mbedtls_mpi_uint *A,
1095 : : size_t A_limbs,
1096 : : const mbedtls_mpi_uint *N,
1097 : : size_t N_limbs,
1098 : : mbedtls_mpi_uint *T)
1099 : : {
1100 : : /* GCD and modinv, names common to Alg 7 and Alg 8 */
1101 : 56 : mbedtls_mpi_uint *u = T + 0 * N_limbs;
1102 : 56 : mbedtls_mpi_uint *v = G;
1103 : :
1104 : : /* GCD and modinv, my name (t1, t2 from Alg 7) */
1105 : 56 : mbedtls_mpi_uint *d = T + 1 * N_limbs;
1106 : :
1107 : : /* GCD and modinv, names from Alg 8 (note: t1, t2 from Alg 7 are d above) */
1108 : 56 : mbedtls_mpi_uint *t1 = T + 2 * N_limbs;
1109 : 56 : mbedtls_mpi_uint *t2 = T + 3 * N_limbs;
1110 : :
1111 : : /* modinv only, names common to Alg 7 and Alg 8 */
1112 : 56 : mbedtls_mpi_uint *q = I;
1113 [ + - ]: 56 : mbedtls_mpi_uint *r = I != NULL ? T + 4 * N_limbs : NULL;
1114 : :
1115 : : /*
1116 : : * Initial values:
1117 : : * u, v = A, N
1118 : : * q, r = 0, 1
1119 : : *
1120 : : * We only write to G (aka v) after reading from inputs (A and N), which
1121 : : * allows aliasing, except with N when I != NULL, as then we'll be operating
1122 : : * mod N on q and r later - see the public documentation.
1123 : : */
1124 : 56 : if (A_limbs > N_limbs) {
1125 : : /* Violating this precondition should not result in memory errors. */
1126 : : A_limbs = N_limbs;
1127 : : }
1128 [ + - ]: 56 : memcpy(u, A, A_limbs * ciL);
1129 [ + - ]: 56 : memset((char *) u + A_limbs * ciL, 0, (N_limbs - A_limbs) * ciL);
1130 : :
1131 : : /* Avoid possible UB with memcpy when src == dst. */
1132 [ + - ]: 56 : if (v != N) {
1133 : 56 : memcpy(v, N, N_limbs * ciL);
1134 : : }
1135 : :
1136 [ + - ]: 56 : if (I != NULL) {
1137 : 56 : memset(q, 0, N_limbs * ciL);
1138 : :
1139 : 56 : memset(r, 0, N_limbs * ciL);
1140 : 56 : r[0] = 1;
1141 : : }
1142 : :
1143 : : /*
1144 : : * At each step, out of u, v, v - u we keep one, shift another, and discard
1145 : : * the third, then update (u, v) with the ordered result.
1146 : : * Then we mirror those actions with q, r, r - q mod N.
1147 : : *
1148 : : * Loop invariants:
1149 : : * u <= v (on entry: A <= N)
1150 : : * GCD(u, v) == GCD(A, N) (on entry: trivial)
1151 : : * v = A * q mod N (on entry: N = A * 0 mod N)
1152 : : * u = A * r mod N (on entry: A = A * 1 mod N)
1153 : : * q, r in [0, N) (on entry: 0, 1)
1154 : : *
1155 : : * On exit:
1156 : : * u = 0
1157 : : * v = GCD(A, N) = A * q mod N
1158 : : * if v == 1 then 1 = A * q mod N ie q is A's inverse mod N
1159 : : * r = 0
1160 : : *
1161 : : * The exit state is a fixed point of the loop's body.
1162 : : * Alg 7 and Alg 8 use 2 * bitlen(N) iterations but Theorem 2 (above in the
1163 : : * paper) says bitlen(A) + bitlen(N) is actually enough.
1164 : : */
1165 [ + + ]: 28728 : for (size_t i = 0; i < (A_limbs + N_limbs) * biL; i++) {
1166 : : /* s, z in Alg 8 - use meaningful names instead */
1167 : 28672 : mbedtls_ct_condition_t u_odd = mbedtls_ct_bool(u[0] & 1);
1168 : 28672 : mbedtls_ct_condition_t v_odd = mbedtls_ct_bool(v[0] & 1);
1169 : :
1170 : : /* Other conditions that will be useful below */
1171 : 28672 : mbedtls_ct_condition_t u_odd_v_odd = mbedtls_ct_bool_and(u_odd, v_odd);
1172 : 28672 : mbedtls_ct_condition_t v_even = mbedtls_ct_bool_not(v_odd);
1173 : 28672 : mbedtls_ct_condition_t u_odd_v_even = mbedtls_ct_bool_and(u_odd, v_even);
1174 : :
1175 : : /* This is called t1 in Alg 7 (no name in Alg 8).
1176 : : * We know that u <= v so there is no carry */
1177 : 28672 : (void) mbedtls_mpi_core_sub(d, v, u, N_limbs);
1178 : :
1179 : : /* t1 (the thing that's kept) can be d (default) or u (if t2 is d) */
1180 : 28672 : memcpy(t1, d, N_limbs * ciL);
1181 : 28672 : mbedtls_mpi_core_cond_assign(t1, u, N_limbs, u_odd_v_odd);
1182 : :
1183 : : /* t2 (the thing that's shifted) can be u (if even), or v (if even),
1184 : : * or d (which is even if both u and v were odd) */
1185 : 28672 : memcpy(t2, u, N_limbs * ciL);
1186 : 28672 : mbedtls_mpi_core_cond_assign(t2, v, N_limbs, u_odd_v_even);
1187 : 28672 : mbedtls_mpi_core_cond_assign(t2, d, N_limbs, u_odd_v_odd);
1188 : :
1189 : 28672 : mbedtls_mpi_core_shift_r(t2, N_limbs, 1); // t2 is even
1190 : :
1191 : : /* Update u, v and re-order them if needed */
1192 : 28672 : memcpy(u, t1, N_limbs * ciL);
1193 : 28672 : memcpy(v, t2, N_limbs * ciL);
1194 : 28672 : mbedtls_ct_condition_t swap = mbedtls_mpi_core_lt_ct(v, u, N_limbs);
1195 : 28672 : mbedtls_mpi_core_cond_swap(u, v, N_limbs, swap);
1196 : :
1197 : : /* Now, if modinv was requested, do the same with q, r, but:
1198 : : * - decisions still based on u and v (their initial values);
1199 : : * - operations are now mod N;
1200 : : * - we re-use t1, t2 for what the paper calls t3, t4 in Alg 8.
1201 : : *
1202 : : * Here we slightly diverge from the paper and instead do the obvious
1203 : : * thing that preserves the invariants involving q and r: mirror
1204 : : * operations on u and v, ie also divide by 2 here (mod N).
1205 : : *
1206 : : * The paper uses a trick where it replaces division by 2 with
1207 : : * multiplication by 2 here, and compensates in the end by multiplying
1208 : : * by pre_com, which is probably intended as an optimisation.
1209 : : *
1210 : : * However I believe it's not actually an optimisation, since
1211 : : * constant-time modular multiplication by 2 (left-shift + conditional
1212 : : * subtract) is just as costly as constant-time modular division by 2
1213 : : * (conditional add + right-shift). So, skip it and keep things simple.
1214 : : */
1215 [ + - ]: 28672 : if (I != NULL) {
1216 : : /* This is called t2 in Alg 7 (no name in Alg 8). */
1217 : 28672 : mpi_core_sub_mod(d, q, r, N, N_limbs);
1218 : :
1219 : : /* t3 (the thing that's kept) */
1220 : 28672 : memcpy(t1, d, N_limbs * ciL);
1221 : 28672 : mbedtls_mpi_core_cond_assign(t1, r, N_limbs, u_odd_v_odd);
1222 : :
1223 : : /* t4 (the thing that's shifted) */
1224 : 28672 : memcpy(t2, r, N_limbs * ciL);
1225 : 28672 : mbedtls_mpi_core_cond_assign(t2, q, N_limbs, u_odd_v_even);
1226 : 28672 : mbedtls_mpi_core_cond_assign(t2, d, N_limbs, u_odd_v_odd);
1227 : :
1228 : 28672 : mbedtls_mpi_core_div2_mod_odd(t2, N, N_limbs);
1229 : :
1230 : : /* Update and possibly swap */
1231 : 28672 : memcpy(r, t1, N_limbs * ciL);
1232 : 28672 : memcpy(q, t2, N_limbs * ciL);
1233 : 28672 : mbedtls_mpi_core_cond_swap(r, q, N_limbs, swap);
1234 : : }
1235 : : }
1236 : :
1237 : : /* G and I already hold the correct values by virtue of being aliased */
1238 : 56 : }
1239 : :
1240 : : #endif /* MBEDTLS_BIGNUM_C */
|
