xref: /third_party/mbedtls/library/bignum.c (revision a8e1175b)
1/*
2 *  Multi-precision integer library
3 *
4 *  Copyright The Mbed TLS Contributors
5 *  SPDX-License-Identifier: Apache-2.0 OR GPL-2.0-or-later
6 */
7
8/*
9 *  The following sources were referenced in the design of this Multi-precision
10 *  Integer library:
11 *
12 *  [1] Handbook of Applied Cryptography - 1997
13 *      Menezes, van Oorschot and Vanstone
14 *
15 *  [2] Multi-Precision Math
16 *      Tom St Denis
17 *      https://github.com/libtom/libtommath/blob/develop/tommath.pdf
18 *
19 *  [3] GNU Multi-Precision Arithmetic Library
20 *      https://gmplib.org/manual/index.html
21 *
22 */
23
24#include "common.h"
25
26#if defined(MBEDTLS_BIGNUM_C)
27
28#include "mbedtls/bignum.h"
29#include "bignum_core.h"
30#include "bn_mul.h"
31#include "mbedtls/platform_util.h"
32#include "mbedtls/error.h"
33#include "constant_time_internal.h"
34
35#include <limits.h>
36#include <string.h>
37
38#include "mbedtls/platform.h"
39
40
41
42/*
43 * Conditionally select an MPI sign in constant time.
44 * (MPI sign is the field s in mbedtls_mpi. It is unsigned short and only 1 and -1 are valid
45 * values.)
46 */
47static inline signed short mbedtls_ct_mpi_sign_if(mbedtls_ct_condition_t cond,
48                                                  signed short sign1, signed short sign2)
49{
50    return (signed short) mbedtls_ct_uint_if(cond, sign1 + 1, sign2 + 1) - 1;
51}
52
53/*
54 * Compare signed values in constant time
55 */
56int mbedtls_mpi_lt_mpi_ct(const mbedtls_mpi *X,
57                          const mbedtls_mpi *Y,
58                          unsigned *ret)
59{
60    mbedtls_ct_condition_t different_sign, X_is_negative, Y_is_negative, result;
61
62    if (X->n != Y->n) {
63        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
64    }
65
66    /*
67     * Set N_is_negative to MBEDTLS_CT_FALSE if N >= 0, MBEDTLS_CT_TRUE if N < 0.
68     * We know that N->s == 1 if N >= 0 and N->s == -1 if N < 0.
69     */
70    X_is_negative = mbedtls_ct_bool((X->s & 2) >> 1);
71    Y_is_negative = mbedtls_ct_bool((Y->s & 2) >> 1);
72
73    /*
74     * If the signs are different, then the positive operand is the bigger.
75     * That is if X is negative (X_is_negative == 1), then X < Y is true and it
76     * is false if X is positive (X_is_negative == 0).
77     */
78    different_sign = mbedtls_ct_bool_ne(X_is_negative, Y_is_negative); // true if different sign
79    result = mbedtls_ct_bool_and(different_sign, X_is_negative);
80
81    /*
82     * Assuming signs are the same, compare X and Y. We switch the comparison
83     * order if they are negative so that we get the right result, regardles of
84     * sign.
85     */
86
87    /* This array is used to conditionally swap the pointers in const time */
88    void * const p[2] = { X->p, Y->p };
89    size_t i = mbedtls_ct_size_if_else_0(X_is_negative, 1);
90    mbedtls_ct_condition_t lt = mbedtls_mpi_core_lt_ct(p[i], p[i ^ 1], X->n);
91
92    /*
93     * Store in result iff the signs are the same (i.e., iff different_sign == false). If
94     * the signs differ, result has already been set, so we don't change it.
95     */
96    result = mbedtls_ct_bool_or(result,
97                                mbedtls_ct_bool_and(mbedtls_ct_bool_not(different_sign), lt));
98
99    *ret = mbedtls_ct_uint_if_else_0(result, 1);
100
101    return 0;
102}
103
104/*
105 * Conditionally assign X = Y, without leaking information
106 * about whether the assignment was made or not.
107 * (Leaking information about the respective sizes of X and Y is ok however.)
108 */
109#if defined(_MSC_VER) && defined(MBEDTLS_PLATFORM_IS_WINDOWS_ON_ARM64) && \
110    (_MSC_FULL_VER < 193131103)
111/*
112 * MSVC miscompiles this function if it's inlined prior to Visual Studio 2022 version 17.1. See:
113 * https://developercommunity.visualstudio.com/t/c-compiler-miscompiles-part-of-mbedtls-library-on/1646989
114 */
115__declspec(noinline)
116#endif
117int mbedtls_mpi_safe_cond_assign(mbedtls_mpi *X,
118                                 const mbedtls_mpi *Y,
119                                 unsigned char assign)
120{
121    int ret = 0;
122
123    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, Y->n));
124
125    {
126        mbedtls_ct_condition_t do_assign = mbedtls_ct_bool(assign);
127
128        X->s = mbedtls_ct_mpi_sign_if(do_assign, Y->s, X->s);
129
130        mbedtls_mpi_core_cond_assign(X->p, Y->p, Y->n, do_assign);
131
132        mbedtls_ct_condition_t do_not_assign = mbedtls_ct_bool_not(do_assign);
133        for (size_t i = Y->n; i < X->n; i++) {
134            X->p[i] = mbedtls_ct_mpi_uint_if_else_0(do_not_assign, X->p[i]);
135        }
136    }
137
138cleanup:
139    return ret;
140}
141
142/*
143 * Conditionally swap X and Y, without leaking information
144 * about whether the swap was made or not.
145 * Here it is not ok to simply swap the pointers, which would lead to
146 * different memory access patterns when X and Y are used afterwards.
147 */
148int mbedtls_mpi_safe_cond_swap(mbedtls_mpi *X,
149                               mbedtls_mpi *Y,
150                               unsigned char swap)
151{
152    int ret = 0;
153    int s;
154
155    if (X == Y) {
156        return 0;
157    }
158
159    mbedtls_ct_condition_t do_swap = mbedtls_ct_bool(swap);
160
161    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, Y->n));
162    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(Y, X->n));
163
164    s = X->s;
165    X->s = mbedtls_ct_mpi_sign_if(do_swap, Y->s, X->s);
166    Y->s = mbedtls_ct_mpi_sign_if(do_swap, s, Y->s);
167
168    mbedtls_mpi_core_cond_swap(X->p, Y->p, X->n, do_swap);
169
170cleanup:
171    return ret;
172}
173
174/* Implementation that should never be optimized out by the compiler */
175#define mbedtls_mpi_zeroize_and_free(v, n) mbedtls_zeroize_and_free(v, ciL * (n))
176
177/*
178 * Initialize one MPI
179 */
180void mbedtls_mpi_init(mbedtls_mpi *X)
181{
182    X->s = 1;
183    X->n = 0;
184    X->p = NULL;
185}
186
187/*
188 * Unallocate one MPI
189 */
190void mbedtls_mpi_free(mbedtls_mpi *X)
191{
192    if (X == NULL) {
193        return;
194    }
195
196    if (X->p != NULL) {
197        mbedtls_mpi_zeroize_and_free(X->p, X->n);
198    }
199
200    X->s = 1;
201    X->n = 0;
202    X->p = NULL;
203}
204
205/*
206 * Enlarge to the specified number of limbs
207 */
208int mbedtls_mpi_grow(mbedtls_mpi *X, size_t nblimbs)
209{
210    mbedtls_mpi_uint *p;
211
212    if (nblimbs > MBEDTLS_MPI_MAX_LIMBS) {
213        return MBEDTLS_ERR_MPI_ALLOC_FAILED;
214    }
215
216    if (X->n < nblimbs) {
217        if ((p = (mbedtls_mpi_uint *) mbedtls_calloc(nblimbs, ciL)) == NULL) {
218            return MBEDTLS_ERR_MPI_ALLOC_FAILED;
219        }
220
221        if (X->p != NULL) {
222            memcpy(p, X->p, X->n * ciL);
223            mbedtls_mpi_zeroize_and_free(X->p, X->n);
224        }
225
226        /* nblimbs fits in n because we ensure that MBEDTLS_MPI_MAX_LIMBS
227         * fits, and we've checked that nblimbs <= MBEDTLS_MPI_MAX_LIMBS. */
228        X->n = (unsigned short) nblimbs;
229        X->p = p;
230    }
231
232    return 0;
233}
234
235/*
236 * Resize down as much as possible,
237 * while keeping at least the specified number of limbs
238 */
239int mbedtls_mpi_shrink(mbedtls_mpi *X, size_t nblimbs)
240{
241    mbedtls_mpi_uint *p;
242    size_t i;
243
244    if (nblimbs > MBEDTLS_MPI_MAX_LIMBS) {
245        return MBEDTLS_ERR_MPI_ALLOC_FAILED;
246    }
247
248    /* Actually resize up if there are currently fewer than nblimbs limbs. */
249    if (X->n <= nblimbs) {
250        return mbedtls_mpi_grow(X, nblimbs);
251    }
252    /* After this point, then X->n > nblimbs and in particular X->n > 0. */
253
254    for (i = X->n - 1; i > 0; i--) {
255        if (X->p[i] != 0) {
256            break;
257        }
258    }
259    i++;
260
261    if (i < nblimbs) {
262        i = nblimbs;
263    }
264
265    if ((p = (mbedtls_mpi_uint *) mbedtls_calloc(i, ciL)) == NULL) {
266        return MBEDTLS_ERR_MPI_ALLOC_FAILED;
267    }
268
269    if (X->p != NULL) {
270        memcpy(p, X->p, i * ciL);
271        mbedtls_mpi_zeroize_and_free(X->p, X->n);
272    }
273
274    /* i fits in n because we ensure that MBEDTLS_MPI_MAX_LIMBS
275     * fits, and we've checked that i <= nblimbs <= MBEDTLS_MPI_MAX_LIMBS. */
276    X->n = (unsigned short) i;
277    X->p = p;
278
279    return 0;
280}
281
282/* Resize X to have exactly n limbs and set it to 0. */
283static int mbedtls_mpi_resize_clear(mbedtls_mpi *X, size_t limbs)
284{
285    if (limbs == 0) {
286        mbedtls_mpi_free(X);
287        return 0;
288    } else if (X->n == limbs) {
289        memset(X->p, 0, limbs * ciL);
290        X->s = 1;
291        return 0;
292    } else {
293        mbedtls_mpi_free(X);
294        return mbedtls_mpi_grow(X, limbs);
295    }
296}
297
298/*
299 * Copy the contents of Y into X.
300 *
301 * This function is not constant-time. Leading zeros in Y may be removed.
302 *
303 * Ensure that X does not shrink. This is not guaranteed by the public API,
304 * but some code in the bignum module might still rely on this property.
305 */
306int mbedtls_mpi_copy(mbedtls_mpi *X, const mbedtls_mpi *Y)
307{
308    int ret = 0;
309    size_t i;
310
311    if (X == Y) {
312        return 0;
313    }
314
315    if (Y->n == 0) {
316        if (X->n != 0) {
317            X->s = 1;
318            memset(X->p, 0, X->n * ciL);
319        }
320        return 0;
321    }
322
323    for (i = Y->n - 1; i > 0; i--) {
324        if (Y->p[i] != 0) {
325            break;
326        }
327    }
328    i++;
329
330    X->s = Y->s;
331
332    if (X->n < i) {
333        MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, i));
334    } else {
335        memset(X->p + i, 0, (X->n - i) * ciL);
336    }
337
338    memcpy(X->p, Y->p, i * ciL);
339
340cleanup:
341
342    return ret;
343}
344
345/*
346 * Swap the contents of X and Y
347 */
348void mbedtls_mpi_swap(mbedtls_mpi *X, mbedtls_mpi *Y)
349{
350    mbedtls_mpi T;
351
352    memcpy(&T,  X, sizeof(mbedtls_mpi));
353    memcpy(X,  Y, sizeof(mbedtls_mpi));
354    memcpy(Y, &T, sizeof(mbedtls_mpi));
355}
356
357static inline mbedtls_mpi_uint mpi_sint_abs(mbedtls_mpi_sint z)
358{
359    if (z >= 0) {
360        return z;
361    }
362    /* Take care to handle the most negative value (-2^(biL-1)) correctly.
363     * A naive -z would have undefined behavior.
364     * Write this in a way that makes popular compilers happy (GCC, Clang,
365     * MSVC). */
366    return (mbedtls_mpi_uint) 0 - (mbedtls_mpi_uint) z;
367}
368
369/* Convert x to a sign, i.e. to 1, if x is positive, or -1, if x is negative.
370 * This looks awkward but generates smaller code than (x < 0 ? -1 : 1) */
371#define TO_SIGN(x) ((mbedtls_mpi_sint) (((mbedtls_mpi_uint) x) >> (biL - 1)) * -2 + 1)
372
373/*
374 * Set value from integer
375 */
376int mbedtls_mpi_lset(mbedtls_mpi *X, mbedtls_mpi_sint z)
377{
378    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
379
380    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, 1));
381    memset(X->p, 0, X->n * ciL);
382
383    X->p[0] = mpi_sint_abs(z);
384    X->s    = TO_SIGN(z);
385
386cleanup:
387
388    return ret;
389}
390
391/*
392 * Get a specific bit
393 */
394int mbedtls_mpi_get_bit(const mbedtls_mpi *X, size_t pos)
395{
396    if (X->n * biL <= pos) {
397        return 0;
398    }
399
400    return (X->p[pos / biL] >> (pos % biL)) & 0x01;
401}
402
403/*
404 * Set a bit to a specific value of 0 or 1
405 */
406int mbedtls_mpi_set_bit(mbedtls_mpi *X, size_t pos, unsigned char val)
407{
408    int ret = 0;
409    size_t off = pos / biL;
410    size_t idx = pos % biL;
411
412    if (val != 0 && val != 1) {
413        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
414    }
415
416    if (X->n * biL <= pos) {
417        if (val == 0) {
418            return 0;
419        }
420
421        MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, off + 1));
422    }
423
424    X->p[off] &= ~((mbedtls_mpi_uint) 0x01 << idx);
425    X->p[off] |= (mbedtls_mpi_uint) val << idx;
426
427cleanup:
428
429    return ret;
430}
431
432/*
433 * Return the number of less significant zero-bits
434 */
435size_t mbedtls_mpi_lsb(const mbedtls_mpi *X)
436{
437    size_t i;
438
439#if defined(__has_builtin)
440#if (MBEDTLS_MPI_UINT_MAX == UINT_MAX) && __has_builtin(__builtin_ctz)
441    #define mbedtls_mpi_uint_ctz __builtin_ctz
442#elif (MBEDTLS_MPI_UINT_MAX == ULONG_MAX) && __has_builtin(__builtin_ctzl)
443    #define mbedtls_mpi_uint_ctz __builtin_ctzl
444#elif (MBEDTLS_MPI_UINT_MAX == ULLONG_MAX) && __has_builtin(__builtin_ctzll)
445    #define mbedtls_mpi_uint_ctz __builtin_ctzll
446#endif
447#endif
448
449#if defined(mbedtls_mpi_uint_ctz)
450    for (i = 0; i < X->n; i++) {
451        if (X->p[i] != 0) {
452            return i * biL + mbedtls_mpi_uint_ctz(X->p[i]);
453        }
454    }
455#else
456    size_t count = 0;
457    for (i = 0; i < X->n; i++) {
458        for (size_t j = 0; j < biL; j++, count++) {
459            if (((X->p[i] >> j) & 1) != 0) {
460                return count;
461            }
462        }
463    }
464#endif
465
466    return 0;
467}
468
469/*
470 * Return the number of bits
471 */
472size_t mbedtls_mpi_bitlen(const mbedtls_mpi *X)
473{
474    return mbedtls_mpi_core_bitlen(X->p, X->n);
475}
476
477/*
478 * Return the total size in bytes
479 */
480size_t mbedtls_mpi_size(const mbedtls_mpi *X)
481{
482    return (mbedtls_mpi_bitlen(X) + 7) >> 3;
483}
484
485/*
486 * Convert an ASCII character to digit value
487 */
488static int mpi_get_digit(mbedtls_mpi_uint *d, int radix, char c)
489{
490    *d = 255;
491
492    if (c >= 0x30 && c <= 0x39) {
493        *d = c - 0x30;
494    }
495    if (c >= 0x41 && c <= 0x46) {
496        *d = c - 0x37;
497    }
498    if (c >= 0x61 && c <= 0x66) {
499        *d = c - 0x57;
500    }
501
502    if (*d >= (mbedtls_mpi_uint) radix) {
503        return MBEDTLS_ERR_MPI_INVALID_CHARACTER;
504    }
505
506    return 0;
507}
508
509/*
510 * Import from an ASCII string
511 */
512int mbedtls_mpi_read_string(mbedtls_mpi *X, int radix, const char *s)
513{
514    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
515    size_t i, j, slen, n;
516    int sign = 1;
517    mbedtls_mpi_uint d;
518    mbedtls_mpi T;
519
520    if (radix < 2 || radix > 16) {
521        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
522    }
523
524    mbedtls_mpi_init(&T);
525
526    if (s[0] == 0) {
527        mbedtls_mpi_free(X);
528        return 0;
529    }
530
531    if (s[0] == '-') {
532        ++s;
533        sign = -1;
534    }
535
536    slen = strlen(s);
537
538    if (radix == 16) {
539        if (slen > SIZE_MAX >> 2) {
540            return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
541        }
542
543        n = BITS_TO_LIMBS(slen << 2);
544
545        MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, n));
546        MBEDTLS_MPI_CHK(mbedtls_mpi_lset(X, 0));
547
548        for (i = slen, j = 0; i > 0; i--, j++) {
549            MBEDTLS_MPI_CHK(mpi_get_digit(&d, radix, s[i - 1]));
550            X->p[j / (2 * ciL)] |= d << ((j % (2 * ciL)) << 2);
551        }
552    } else {
553        MBEDTLS_MPI_CHK(mbedtls_mpi_lset(X, 0));
554
555        for (i = 0; i < slen; i++) {
556            MBEDTLS_MPI_CHK(mpi_get_digit(&d, radix, s[i]));
557            MBEDTLS_MPI_CHK(mbedtls_mpi_mul_int(&T, X, radix));
558            MBEDTLS_MPI_CHK(mbedtls_mpi_add_int(X, &T, d));
559        }
560    }
561
562    if (sign < 0 && mbedtls_mpi_bitlen(X) != 0) {
563        X->s = -1;
564    }
565
566cleanup:
567
568    mbedtls_mpi_free(&T);
569
570    return ret;
571}
572
573/*
574 * Helper to write the digits high-order first.
575 */
576static int mpi_write_hlp(mbedtls_mpi *X, int radix,
577                         char **p, const size_t buflen)
578{
579    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
580    mbedtls_mpi_uint r;
581    size_t length = 0;
582    char *p_end = *p + buflen;
583
584    do {
585        if (length >= buflen) {
586            return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
587        }
588
589        MBEDTLS_MPI_CHK(mbedtls_mpi_mod_int(&r, X, radix));
590        MBEDTLS_MPI_CHK(mbedtls_mpi_div_int(X, NULL, X, radix));
591        /*
592         * Write the residue in the current position, as an ASCII character.
593         */
594        if (r < 0xA) {
595            *(--p_end) = (char) ('0' + r);
596        } else {
597            *(--p_end) = (char) ('A' + (r - 0xA));
598        }
599
600        length++;
601    } while (mbedtls_mpi_cmp_int(X, 0) != 0);
602
603    memmove(*p, p_end, length);
604    *p += length;
605
606cleanup:
607
608    return ret;
609}
610
611/*
612 * Export into an ASCII string
613 */
614int mbedtls_mpi_write_string(const mbedtls_mpi *X, int radix,
615                             char *buf, size_t buflen, size_t *olen)
616{
617    int ret = 0;
618    size_t n;
619    char *p;
620    mbedtls_mpi T;
621
622    if (radix < 2 || radix > 16) {
623        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
624    }
625
626    n = mbedtls_mpi_bitlen(X);   /* Number of bits necessary to present `n`. */
627    if (radix >=  4) {
628        n >>= 1;                 /* Number of 4-adic digits necessary to present
629                                  * `n`. If radix > 4, this might be a strict
630                                  * overapproximation of the number of
631                                  * radix-adic digits needed to present `n`. */
632    }
633    if (radix >= 16) {
634        n >>= 1;                 /* Number of hexadecimal digits necessary to
635                                  * present `n`. */
636
637    }
638    n += 1; /* Terminating null byte */
639    n += 1; /* Compensate for the divisions above, which round down `n`
640             * in case it's not even. */
641    n += 1; /* Potential '-'-sign. */
642    n += (n & 1);   /* Make n even to have enough space for hexadecimal writing,
643                     * which always uses an even number of hex-digits. */
644
645    if (buflen < n) {
646        *olen = n;
647        return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
648    }
649
650    p = buf;
651    mbedtls_mpi_init(&T);
652
653    if (X->s == -1) {
654        *p++ = '-';
655        buflen--;
656    }
657
658    if (radix == 16) {
659        int c;
660        size_t i, j, k;
661
662        for (i = X->n, k = 0; i > 0; i--) {
663            for (j = ciL; j > 0; j--) {
664                c = (X->p[i - 1] >> ((j - 1) << 3)) & 0xFF;
665
666                if (c == 0 && k == 0 && (i + j) != 2) {
667                    continue;
668                }
669
670                *(p++) = "0123456789ABCDEF" [c / 16];
671                *(p++) = "0123456789ABCDEF" [c % 16];
672                k = 1;
673            }
674        }
675    } else {
676        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&T, X));
677
678        if (T.s == -1) {
679            T.s = 1;
680        }
681
682        MBEDTLS_MPI_CHK(mpi_write_hlp(&T, radix, &p, buflen));
683    }
684
685    *p++ = '\0';
686    *olen = (size_t) (p - buf);
687
688cleanup:
689
690    mbedtls_mpi_free(&T);
691
692    return ret;
693}
694
695#if defined(MBEDTLS_FS_IO)
696/*
697 * Read X from an opened file
698 */
699int mbedtls_mpi_read_file(mbedtls_mpi *X, int radix, FILE *fin)
700{
701    mbedtls_mpi_uint d;
702    size_t slen;
703    char *p;
704    /*
705     * Buffer should have space for (short) label and decimal formatted MPI,
706     * newline characters and '\0'
707     */
708    char s[MBEDTLS_MPI_RW_BUFFER_SIZE];
709
710    if (radix < 2 || radix > 16) {
711        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
712    }
713
714    memset(s, 0, sizeof(s));
715    if (fgets(s, sizeof(s) - 1, fin) == NULL) {
716        return MBEDTLS_ERR_MPI_FILE_IO_ERROR;
717    }
718
719    slen = strlen(s);
720    if (slen == sizeof(s) - 2) {
721        return MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL;
722    }
723
724    if (slen > 0 && s[slen - 1] == '\n') {
725        slen--; s[slen] = '\0';
726    }
727    if (slen > 0 && s[slen - 1] == '\r') {
728        slen--; s[slen] = '\0';
729    }
730
731    p = s + slen;
732    while (p-- > s) {
733        if (mpi_get_digit(&d, radix, *p) != 0) {
734            break;
735        }
736    }
737
738    return mbedtls_mpi_read_string(X, radix, p + 1);
739}
740
741/*
742 * Write X into an opened file (or stdout if fout == NULL)
743 */
744int mbedtls_mpi_write_file(const char *p, const mbedtls_mpi *X, int radix, FILE *fout)
745{
746    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
747    size_t n, slen, plen;
748    /*
749     * Buffer should have space for (short) label and decimal formatted MPI,
750     * newline characters and '\0'
751     */
752    char s[MBEDTLS_MPI_RW_BUFFER_SIZE];
753
754    if (radix < 2 || radix > 16) {
755        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
756    }
757
758    memset(s, 0, sizeof(s));
759
760    MBEDTLS_MPI_CHK(mbedtls_mpi_write_string(X, radix, s, sizeof(s) - 2, &n));
761
762    if (p == NULL) {
763        p = "";
764    }
765
766    plen = strlen(p);
767    slen = strlen(s);
768    s[slen++] = '\r';
769    s[slen++] = '\n';
770
771    if (fout != NULL) {
772        if (fwrite(p, 1, plen, fout) != plen ||
773            fwrite(s, 1, slen, fout) != slen) {
774            return MBEDTLS_ERR_MPI_FILE_IO_ERROR;
775        }
776    } else {
777        mbedtls_printf("%s%s", p, s);
778    }
779
780cleanup:
781
782    return ret;
783}
784#endif /* MBEDTLS_FS_IO */
785
786/*
787 * Import X from unsigned binary data, little endian
788 *
789 * This function is guaranteed to return an MPI with exactly the necessary
790 * number of limbs (in particular, it does not skip 0s in the input).
791 */
792int mbedtls_mpi_read_binary_le(mbedtls_mpi *X,
793                               const unsigned char *buf, size_t buflen)
794{
795    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
796    const size_t limbs = CHARS_TO_LIMBS(buflen);
797
798    /* Ensure that target MPI has exactly the necessary number of limbs */
799    MBEDTLS_MPI_CHK(mbedtls_mpi_resize_clear(X, limbs));
800
801    MBEDTLS_MPI_CHK(mbedtls_mpi_core_read_le(X->p, X->n, buf, buflen));
802
803cleanup:
804
805    /*
806     * This function is also used to import keys. However, wiping the buffers
807     * upon failure is not necessary because failure only can happen before any
808     * input is copied.
809     */
810    return ret;
811}
812
813/*
814 * Import X from unsigned binary data, big endian
815 *
816 * This function is guaranteed to return an MPI with exactly the necessary
817 * number of limbs (in particular, it does not skip 0s in the input).
818 */
819int mbedtls_mpi_read_binary(mbedtls_mpi *X, const unsigned char *buf, size_t buflen)
820{
821    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
822    const size_t limbs = CHARS_TO_LIMBS(buflen);
823
824    /* Ensure that target MPI has exactly the necessary number of limbs */
825    MBEDTLS_MPI_CHK(mbedtls_mpi_resize_clear(X, limbs));
826
827    MBEDTLS_MPI_CHK(mbedtls_mpi_core_read_be(X->p, X->n, buf, buflen));
828
829cleanup:
830
831    /*
832     * This function is also used to import keys. However, wiping the buffers
833     * upon failure is not necessary because failure only can happen before any
834     * input is copied.
835     */
836    return ret;
837}
838
839/*
840 * Export X into unsigned binary data, little endian
841 */
842int mbedtls_mpi_write_binary_le(const mbedtls_mpi *X,
843                                unsigned char *buf, size_t buflen)
844{
845    return mbedtls_mpi_core_write_le(X->p, X->n, buf, buflen);
846}
847
848/*
849 * Export X into unsigned binary data, big endian
850 */
851int mbedtls_mpi_write_binary(const mbedtls_mpi *X,
852                             unsigned char *buf, size_t buflen)
853{
854    return mbedtls_mpi_core_write_be(X->p, X->n, buf, buflen);
855}
856
857/*
858 * Left-shift: X <<= count
859 */
860int mbedtls_mpi_shift_l(mbedtls_mpi *X, size_t count)
861{
862    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
863    size_t i;
864
865    i = mbedtls_mpi_bitlen(X) + count;
866
867    if (X->n * biL < i) {
868        MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, BITS_TO_LIMBS(i)));
869    }
870
871    ret = 0;
872
873    mbedtls_mpi_core_shift_l(X->p, X->n, count);
874cleanup:
875
876    return ret;
877}
878
879/*
880 * Right-shift: X >>= count
881 */
882int mbedtls_mpi_shift_r(mbedtls_mpi *X, size_t count)
883{
884    if (X->n != 0) {
885        mbedtls_mpi_core_shift_r(X->p, X->n, count);
886    }
887    return 0;
888}
889
890/*
891 * Compare unsigned values
892 */
893int mbedtls_mpi_cmp_abs(const mbedtls_mpi *X, const mbedtls_mpi *Y)
894{
895    size_t i, j;
896
897    for (i = X->n; i > 0; i--) {
898        if (X->p[i - 1] != 0) {
899            break;
900        }
901    }
902
903    for (j = Y->n; j > 0; j--) {
904        if (Y->p[j - 1] != 0) {
905            break;
906        }
907    }
908
909    /* If i == j == 0, i.e. abs(X) == abs(Y),
910     * we end up returning 0 at the end of the function. */
911
912    if (i > j) {
913        return 1;
914    }
915    if (j > i) {
916        return -1;
917    }
918
919    for (; i > 0; i--) {
920        if (X->p[i - 1] > Y->p[i - 1]) {
921            return 1;
922        }
923        if (X->p[i - 1] < Y->p[i - 1]) {
924            return -1;
925        }
926    }
927
928    return 0;
929}
930
931/*
932 * Compare signed values
933 */
934int mbedtls_mpi_cmp_mpi(const mbedtls_mpi *X, const mbedtls_mpi *Y)
935{
936    size_t i, j;
937
938    for (i = X->n; i > 0; i--) {
939        if (X->p[i - 1] != 0) {
940            break;
941        }
942    }
943
944    for (j = Y->n; j > 0; j--) {
945        if (Y->p[j - 1] != 0) {
946            break;
947        }
948    }
949
950    if (i == 0 && j == 0) {
951        return 0;
952    }
953
954    if (i > j) {
955        return X->s;
956    }
957    if (j > i) {
958        return -Y->s;
959    }
960
961    if (X->s > 0 && Y->s < 0) {
962        return 1;
963    }
964    if (Y->s > 0 && X->s < 0) {
965        return -1;
966    }
967
968    for (; i > 0; i--) {
969        if (X->p[i - 1] > Y->p[i - 1]) {
970            return X->s;
971        }
972        if (X->p[i - 1] < Y->p[i - 1]) {
973            return -X->s;
974        }
975    }
976
977    return 0;
978}
979
980/*
981 * Compare signed values
982 */
983int mbedtls_mpi_cmp_int(const mbedtls_mpi *X, mbedtls_mpi_sint z)
984{
985    mbedtls_mpi Y;
986    mbedtls_mpi_uint p[1];
987
988    *p  = mpi_sint_abs(z);
989    Y.s = TO_SIGN(z);
990    Y.n = 1;
991    Y.p = p;
992
993    return mbedtls_mpi_cmp_mpi(X, &Y);
994}
995
996/*
997 * Unsigned addition: X = |A| + |B|  (HAC 14.7)
998 */
999int mbedtls_mpi_add_abs(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B)
1000{
1001    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1002    size_t j;
1003    mbedtls_mpi_uint *p;
1004    mbedtls_mpi_uint c;
1005
1006    if (X == B) {
1007        const mbedtls_mpi *T = A; A = X; B = T;
1008    }
1009
1010    if (X != A) {
1011        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(X, A));
1012    }
1013
1014    /*
1015     * X must always be positive as a result of unsigned additions.
1016     */
1017    X->s = 1;
1018
1019    for (j = B->n; j > 0; j--) {
1020        if (B->p[j - 1] != 0) {
1021            break;
1022        }
1023    }
1024
1025    /* Exit early to avoid undefined behavior on NULL+0 when X->n == 0
1026     * and B is 0 (of any size). */
1027    if (j == 0) {
1028        return 0;
1029    }
1030
1031    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, j));
1032
1033    /* j is the number of non-zero limbs of B. Add those to X. */
1034
1035    p = X->p;
1036
1037    c = mbedtls_mpi_core_add(p, p, B->p, j);
1038
1039    p += j;
1040
1041    /* Now propagate any carry */
1042
1043    while (c != 0) {
1044        if (j >= X->n) {
1045            MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, j + 1));
1046            p = X->p + j;
1047        }
1048
1049        *p += c; c = (*p < c); j++; p++;
1050    }
1051
1052cleanup:
1053
1054    return ret;
1055}
1056
1057/*
1058 * Unsigned subtraction: X = |A| - |B|  (HAC 14.9, 14.10)
1059 */
1060int mbedtls_mpi_sub_abs(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B)
1061{
1062    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1063    size_t n;
1064    mbedtls_mpi_uint carry;
1065
1066    for (n = B->n; n > 0; n--) {
1067        if (B->p[n - 1] != 0) {
1068            break;
1069        }
1070    }
1071    if (n > A->n) {
1072        /* B >= (2^ciL)^n > A */
1073        ret = MBEDTLS_ERR_MPI_NEGATIVE_VALUE;
1074        goto cleanup;
1075    }
1076
1077    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, A->n));
1078
1079    /* Set the high limbs of X to match A. Don't touch the lower limbs
1080     * because X might be aliased to B, and we must not overwrite the
1081     * significant digits of B. */
1082    if (A->n > n && A != X) {
1083        memcpy(X->p + n, A->p + n, (A->n - n) * ciL);
1084    }
1085    if (X->n > A->n) {
1086        memset(X->p + A->n, 0, (X->n - A->n) * ciL);
1087    }
1088
1089    carry = mbedtls_mpi_core_sub(X->p, A->p, B->p, n);
1090    if (carry != 0) {
1091        /* Propagate the carry through the rest of X. */
1092        carry = mbedtls_mpi_core_sub_int(X->p + n, X->p + n, carry, X->n - n);
1093
1094        /* If we have further carry/borrow, the result is negative. */
1095        if (carry != 0) {
1096            ret = MBEDTLS_ERR_MPI_NEGATIVE_VALUE;
1097            goto cleanup;
1098        }
1099    }
1100
1101    /* X should always be positive as a result of unsigned subtractions. */
1102    X->s = 1;
1103
1104cleanup:
1105    return ret;
1106}
1107
1108/* Common function for signed addition and subtraction.
1109 * Calculate A + B * flip_B where flip_B is 1 or -1.
1110 */
1111static int add_sub_mpi(mbedtls_mpi *X,
1112                       const mbedtls_mpi *A, const mbedtls_mpi *B,
1113                       int flip_B)
1114{
1115    int ret, s;
1116
1117    s = A->s;
1118    if (A->s * B->s * flip_B < 0) {
1119        int cmp = mbedtls_mpi_cmp_abs(A, B);
1120        if (cmp >= 0) {
1121            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(X, A, B));
1122            /* If |A| = |B|, the result is 0 and we must set the sign bit
1123             * to +1 regardless of which of A or B was negative. Otherwise,
1124             * since |A| > |B|, the sign is the sign of A. */
1125            X->s = cmp == 0 ? 1 : s;
1126        } else {
1127            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(X, B, A));
1128            /* Since |A| < |B|, the sign is the opposite of A. */
1129            X->s = -s;
1130        }
1131    } else {
1132        MBEDTLS_MPI_CHK(mbedtls_mpi_add_abs(X, A, B));
1133        X->s = s;
1134    }
1135
1136cleanup:
1137
1138    return ret;
1139}
1140
1141/*
1142 * Signed addition: X = A + B
1143 */
1144int mbedtls_mpi_add_mpi(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B)
1145{
1146    return add_sub_mpi(X, A, B, 1);
1147}
1148
1149/*
1150 * Signed subtraction: X = A - B
1151 */
1152int mbedtls_mpi_sub_mpi(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B)
1153{
1154    return add_sub_mpi(X, A, B, -1);
1155}
1156
1157/*
1158 * Signed addition: X = A + b
1159 */
1160int mbedtls_mpi_add_int(mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b)
1161{
1162    mbedtls_mpi B;
1163    mbedtls_mpi_uint p[1];
1164
1165    p[0] = mpi_sint_abs(b);
1166    B.s = TO_SIGN(b);
1167    B.n = 1;
1168    B.p = p;
1169
1170    return mbedtls_mpi_add_mpi(X, A, &B);
1171}
1172
1173/*
1174 * Signed subtraction: X = A - b
1175 */
1176int mbedtls_mpi_sub_int(mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b)
1177{
1178    mbedtls_mpi B;
1179    mbedtls_mpi_uint p[1];
1180
1181    p[0] = mpi_sint_abs(b);
1182    B.s = TO_SIGN(b);
1183    B.n = 1;
1184    B.p = p;
1185
1186    return mbedtls_mpi_sub_mpi(X, A, &B);
1187}
1188
1189/*
1190 * Baseline multiplication: X = A * B  (HAC 14.12)
1191 */
1192int mbedtls_mpi_mul_mpi(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B)
1193{
1194    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1195    size_t i, j;
1196    mbedtls_mpi TA, TB;
1197    int result_is_zero = 0;
1198
1199    mbedtls_mpi_init(&TA);
1200    mbedtls_mpi_init(&TB);
1201
1202    if (X == A) {
1203        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TA, A)); A = &TA;
1204    }
1205    if (X == B) {
1206        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TB, B)); B = &TB;
1207    }
1208
1209    for (i = A->n; i > 0; i--) {
1210        if (A->p[i - 1] != 0) {
1211            break;
1212        }
1213    }
1214    if (i == 0) {
1215        result_is_zero = 1;
1216    }
1217
1218    for (j = B->n; j > 0; j--) {
1219        if (B->p[j - 1] != 0) {
1220            break;
1221        }
1222    }
1223    if (j == 0) {
1224        result_is_zero = 1;
1225    }
1226
1227    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, i + j));
1228    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(X, 0));
1229
1230    mbedtls_mpi_core_mul(X->p, A->p, i, B->p, j);
1231
1232    /* If the result is 0, we don't shortcut the operation, which reduces
1233     * but does not eliminate side channels leaking the zero-ness. We do
1234     * need to take care to set the sign bit properly since the library does
1235     * not fully support an MPI object with a value of 0 and s == -1. */
1236    if (result_is_zero) {
1237        X->s = 1;
1238    } else {
1239        X->s = A->s * B->s;
1240    }
1241
1242cleanup:
1243
1244    mbedtls_mpi_free(&TB); mbedtls_mpi_free(&TA);
1245
1246    return ret;
1247}
1248
1249/*
1250 * Baseline multiplication: X = A * b
1251 */
1252int mbedtls_mpi_mul_int(mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_uint b)
1253{
1254    size_t n = A->n;
1255    while (n > 0 && A->p[n - 1] == 0) {
1256        --n;
1257    }
1258
1259    /* The general method below doesn't work if b==0. */
1260    if (b == 0 || n == 0) {
1261        return mbedtls_mpi_lset(X, 0);
1262    }
1263
1264    /* Calculate A*b as A + A*(b-1) to take advantage of mbedtls_mpi_core_mla */
1265    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1266    /* In general, A * b requires 1 limb more than b. If
1267     * A->p[n - 1] * b / b == A->p[n - 1], then A * b fits in the same
1268     * number of limbs as A and the call to grow() is not required since
1269     * copy() will take care of the growth if needed. However, experimentally,
1270     * making the call to grow() unconditional causes slightly fewer
1271     * calls to calloc() in ECP code, presumably because it reuses the
1272     * same mpi for a while and this way the mpi is more likely to directly
1273     * grow to its final size.
1274     *
1275     * Note that calculating A*b as 0 + A*b doesn't work as-is because
1276     * A,X can be the same. */
1277    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, n + 1));
1278    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(X, A));
1279    mbedtls_mpi_core_mla(X->p, X->n, A->p, n, b - 1);
1280
1281cleanup:
1282    return ret;
1283}
1284
1285/*
1286 * Unsigned integer divide - double mbedtls_mpi_uint dividend, u1/u0, and
1287 * mbedtls_mpi_uint divisor, d
1288 */
1289static mbedtls_mpi_uint mbedtls_int_div_int(mbedtls_mpi_uint u1,
1290                                            mbedtls_mpi_uint u0,
1291                                            mbedtls_mpi_uint d,
1292                                            mbedtls_mpi_uint *r)
1293{
1294#if defined(MBEDTLS_HAVE_UDBL)
1295    mbedtls_t_udbl dividend, quotient;
1296#else
1297    const mbedtls_mpi_uint radix = (mbedtls_mpi_uint) 1 << biH;
1298    const mbedtls_mpi_uint uint_halfword_mask = ((mbedtls_mpi_uint) 1 << biH) - 1;
1299    mbedtls_mpi_uint d0, d1, q0, q1, rAX, r0, quotient;
1300    mbedtls_mpi_uint u0_msw, u0_lsw;
1301    size_t s;
1302#endif
1303
1304    /*
1305     * Check for overflow
1306     */
1307    if (0 == d || u1 >= d) {
1308        if (r != NULL) {
1309            *r = ~(mbedtls_mpi_uint) 0u;
1310        }
1311
1312        return ~(mbedtls_mpi_uint) 0u;
1313    }
1314
1315#if defined(MBEDTLS_HAVE_UDBL)
1316    dividend  = (mbedtls_t_udbl) u1 << biL;
1317    dividend |= (mbedtls_t_udbl) u0;
1318    quotient = dividend / d;
1319    if (quotient > ((mbedtls_t_udbl) 1 << biL) - 1) {
1320        quotient = ((mbedtls_t_udbl) 1 << biL) - 1;
1321    }
1322
1323    if (r != NULL) {
1324        *r = (mbedtls_mpi_uint) (dividend - (quotient * d));
1325    }
1326
1327    return (mbedtls_mpi_uint) quotient;
1328#else
1329
1330    /*
1331     * Algorithm D, Section 4.3.1 - The Art of Computer Programming
1332     *   Vol. 2 - Seminumerical Algorithms, Knuth
1333     */
1334
1335    /*
1336     * Normalize the divisor, d, and dividend, u0, u1
1337     */
1338    s = mbedtls_mpi_core_clz(d);
1339    d = d << s;
1340
1341    u1 = u1 << s;
1342    u1 |= (u0 >> (biL - s)) & (-(mbedtls_mpi_sint) s >> (biL - 1));
1343    u0 =  u0 << s;
1344
1345    d1 = d >> biH;
1346    d0 = d & uint_halfword_mask;
1347
1348    u0_msw = u0 >> biH;
1349    u0_lsw = u0 & uint_halfword_mask;
1350
1351    /*
1352     * Find the first quotient and remainder
1353     */
1354    q1 = u1 / d1;
1355    r0 = u1 - d1 * q1;
1356
1357    while (q1 >= radix || (q1 * d0 > radix * r0 + u0_msw)) {
1358        q1 -= 1;
1359        r0 += d1;
1360
1361        if (r0 >= radix) {
1362            break;
1363        }
1364    }
1365
1366    rAX = (u1 * radix) + (u0_msw - q1 * d);
1367    q0 = rAX / d1;
1368    r0 = rAX - q0 * d1;
1369
1370    while (q0 >= radix || (q0 * d0 > radix * r0 + u0_lsw)) {
1371        q0 -= 1;
1372        r0 += d1;
1373
1374        if (r0 >= radix) {
1375            break;
1376        }
1377    }
1378
1379    if (r != NULL) {
1380        *r = (rAX * radix + u0_lsw - q0 * d) >> s;
1381    }
1382
1383    quotient = q1 * radix + q0;
1384
1385    return quotient;
1386#endif
1387}
1388
1389/*
1390 * Division by mbedtls_mpi: A = Q * B + R  (HAC 14.20)
1391 */
1392int mbedtls_mpi_div_mpi(mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A,
1393                        const mbedtls_mpi *B)
1394{
1395    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1396    size_t i, n, t, k;
1397    mbedtls_mpi X, Y, Z, T1, T2;
1398    mbedtls_mpi_uint TP2[3];
1399
1400    if (mbedtls_mpi_cmp_int(B, 0) == 0) {
1401        return MBEDTLS_ERR_MPI_DIVISION_BY_ZERO;
1402    }
1403
1404    mbedtls_mpi_init(&X); mbedtls_mpi_init(&Y); mbedtls_mpi_init(&Z);
1405    mbedtls_mpi_init(&T1);
1406    /*
1407     * Avoid dynamic memory allocations for constant-size T2.
1408     *
1409     * T2 is used for comparison only and the 3 limbs are assigned explicitly,
1410     * so nobody increase the size of the MPI and we're safe to use an on-stack
1411     * buffer.
1412     */
1413    T2.s = 1;
1414    T2.n = sizeof(TP2) / sizeof(*TP2);
1415    T2.p = TP2;
1416
1417    if (mbedtls_mpi_cmp_abs(A, B) < 0) {
1418        if (Q != NULL) {
1419            MBEDTLS_MPI_CHK(mbedtls_mpi_lset(Q, 0));
1420        }
1421        if (R != NULL) {
1422            MBEDTLS_MPI_CHK(mbedtls_mpi_copy(R, A));
1423        }
1424        return 0;
1425    }
1426
1427    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&X, A));
1428    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&Y, B));
1429    X.s = Y.s = 1;
1430
1431    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(&Z, A->n + 2));
1432    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&Z,  0));
1433    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(&T1, A->n + 2));
1434
1435    k = mbedtls_mpi_bitlen(&Y) % biL;
1436    if (k < biL - 1) {
1437        k = biL - 1 - k;
1438        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&X, k));
1439        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&Y, k));
1440    } else {
1441        k = 0;
1442    }
1443
1444    n = X.n - 1;
1445    t = Y.n - 1;
1446    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&Y, biL * (n - t)));
1447
1448    while (mbedtls_mpi_cmp_mpi(&X, &Y) >= 0) {
1449        Z.p[n - t]++;
1450        MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&X, &X, &Y));
1451    }
1452    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&Y, biL * (n - t)));
1453
1454    for (i = n; i > t; i--) {
1455        if (X.p[i] >= Y.p[t]) {
1456            Z.p[i - t - 1] = ~(mbedtls_mpi_uint) 0u;
1457        } else {
1458            Z.p[i - t - 1] = mbedtls_int_div_int(X.p[i], X.p[i - 1],
1459                                                 Y.p[t], NULL);
1460        }
1461
1462        T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1463        T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1464        T2.p[2] = X.p[i];
1465
1466        Z.p[i - t - 1]++;
1467        do {
1468            Z.p[i - t - 1]--;
1469
1470            MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&T1, 0));
1471            T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1472            T1.p[1] = Y.p[t];
1473            MBEDTLS_MPI_CHK(mbedtls_mpi_mul_int(&T1, &T1, Z.p[i - t - 1]));
1474        } while (mbedtls_mpi_cmp_mpi(&T1, &T2) > 0);
1475
1476        MBEDTLS_MPI_CHK(mbedtls_mpi_mul_int(&T1, &Y, Z.p[i - t - 1]));
1477        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&T1,  biL * (i - t - 1)));
1478        MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&X, &X, &T1));
1479
1480        if (mbedtls_mpi_cmp_int(&X, 0) < 0) {
1481            MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&T1, &Y));
1482            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&T1, biL * (i - t - 1)));
1483            MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(&X, &X, &T1));
1484            Z.p[i - t - 1]--;
1485        }
1486    }
1487
1488    if (Q != NULL) {
1489        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(Q, &Z));
1490        Q->s = A->s * B->s;
1491    }
1492
1493    if (R != NULL) {
1494        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&X, k));
1495        X.s = A->s;
1496        MBEDTLS_MPI_CHK(mbedtls_mpi_copy(R, &X));
1497
1498        if (mbedtls_mpi_cmp_int(R, 0) == 0) {
1499            R->s = 1;
1500        }
1501    }
1502
1503cleanup:
1504
1505    mbedtls_mpi_free(&X); mbedtls_mpi_free(&Y); mbedtls_mpi_free(&Z);
1506    mbedtls_mpi_free(&T1);
1507    mbedtls_platform_zeroize(TP2, sizeof(TP2));
1508
1509    return ret;
1510}
1511
1512/*
1513 * Division by int: A = Q * b + R
1514 */
1515int mbedtls_mpi_div_int(mbedtls_mpi *Q, mbedtls_mpi *R,
1516                        const mbedtls_mpi *A,
1517                        mbedtls_mpi_sint b)
1518{
1519    mbedtls_mpi B;
1520    mbedtls_mpi_uint p[1];
1521
1522    p[0] = mpi_sint_abs(b);
1523    B.s = TO_SIGN(b);
1524    B.n = 1;
1525    B.p = p;
1526
1527    return mbedtls_mpi_div_mpi(Q, R, A, &B);
1528}
1529
1530/*
1531 * Modulo: R = A mod B
1532 */
1533int mbedtls_mpi_mod_mpi(mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B)
1534{
1535    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1536
1537    if (mbedtls_mpi_cmp_int(B, 0) < 0) {
1538        return MBEDTLS_ERR_MPI_NEGATIVE_VALUE;
1539    }
1540
1541    MBEDTLS_MPI_CHK(mbedtls_mpi_div_mpi(NULL, R, A, B));
1542
1543    while (mbedtls_mpi_cmp_int(R, 0) < 0) {
1544        MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(R, R, B));
1545    }
1546
1547    while (mbedtls_mpi_cmp_mpi(R, B) >= 0) {
1548        MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(R, R, B));
1549    }
1550
1551cleanup:
1552
1553    return ret;
1554}
1555
1556/*
1557 * Modulo: r = A mod b
1558 */
1559int mbedtls_mpi_mod_int(mbedtls_mpi_uint *r, const mbedtls_mpi *A, mbedtls_mpi_sint b)
1560{
1561    size_t i;
1562    mbedtls_mpi_uint x, y, z;
1563
1564    if (b == 0) {
1565        return MBEDTLS_ERR_MPI_DIVISION_BY_ZERO;
1566    }
1567
1568    if (b < 0) {
1569        return MBEDTLS_ERR_MPI_NEGATIVE_VALUE;
1570    }
1571
1572    /*
1573     * handle trivial cases
1574     */
1575    if (b == 1 || A->n == 0) {
1576        *r = 0;
1577        return 0;
1578    }
1579
1580    if (b == 2) {
1581        *r = A->p[0] & 1;
1582        return 0;
1583    }
1584
1585    /*
1586     * general case
1587     */
1588    for (i = A->n, y = 0; i > 0; i--) {
1589        x  = A->p[i - 1];
1590        y  = (y << biH) | (x >> biH);
1591        z  = y / b;
1592        y -= z * b;
1593
1594        x <<= biH;
1595        y  = (y << biH) | (x >> biH);
1596        z  = y / b;
1597        y -= z * b;
1598    }
1599
1600    /*
1601     * If A is negative, then the current y represents a negative value.
1602     * Flipping it to the positive side.
1603     */
1604    if (A->s < 0 && y != 0) {
1605        y = b - y;
1606    }
1607
1608    *r = y;
1609
1610    return 0;
1611}
1612
1613int mbedtls_mpi_exp_mod(mbedtls_mpi *X, const mbedtls_mpi *A,
1614                        const mbedtls_mpi *E, const mbedtls_mpi *N,
1615                        mbedtls_mpi *prec_RR)
1616{
1617    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1618
1619    if (mbedtls_mpi_cmp_int(N, 0) <= 0 || (N->p[0] & 1) == 0) {
1620        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1621    }
1622
1623    if (mbedtls_mpi_cmp_int(E, 0) < 0) {
1624        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1625    }
1626
1627    if (mbedtls_mpi_bitlen(E) > MBEDTLS_MPI_MAX_BITS ||
1628        mbedtls_mpi_bitlen(N) > MBEDTLS_MPI_MAX_BITS) {
1629        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1630    }
1631
1632    /*
1633     * Ensure that the exponent that we are passing to the core is not NULL.
1634     */
1635    if (E->n == 0) {
1636        ret = mbedtls_mpi_lset(X, 1);
1637        return ret;
1638    }
1639
1640    /*
1641     * Allocate working memory for mbedtls_mpi_core_exp_mod()
1642     */
1643    size_t T_limbs = mbedtls_mpi_core_exp_mod_working_limbs(N->n, E->n);
1644    mbedtls_mpi_uint *T = (mbedtls_mpi_uint *) mbedtls_calloc(T_limbs, sizeof(mbedtls_mpi_uint));
1645    if (T == NULL) {
1646        return MBEDTLS_ERR_MPI_ALLOC_FAILED;
1647    }
1648
1649    mbedtls_mpi RR;
1650    mbedtls_mpi_init(&RR);
1651
1652    /*
1653     * If 1st call, pre-compute R^2 mod N
1654     */
1655    if (prec_RR == NULL || prec_RR->p == NULL) {
1656        MBEDTLS_MPI_CHK(mbedtls_mpi_core_get_mont_r2_unsafe(&RR, N));
1657
1658        if (prec_RR != NULL) {
1659            *prec_RR = RR;
1660        }
1661    } else {
1662        MBEDTLS_MPI_CHK(mbedtls_mpi_grow(prec_RR, N->n));
1663        RR = *prec_RR;
1664    }
1665
1666    /*
1667     * To preserve constness we need to make a copy of A. Using X for this to
1668     * save memory.
1669     */
1670    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(X, A));
1671
1672    /*
1673     * Compensate for negative A (and correct at the end).
1674     */
1675    X->s = 1;
1676
1677    /*
1678     * Make sure that X is in a form that is safe for consumption by
1679     * the core functions.
1680     *
1681     * - The core functions will not touch the limbs of X above N->n. The
1682     *   result will be correct if those limbs are 0, which the mod call
1683     *   ensures.
1684     * - Also, X must have at least as many limbs as N for the calls to the
1685     *   core functions.
1686     */
1687    if (mbedtls_mpi_cmp_mpi(X, N) >= 0) {
1688        MBEDTLS_MPI_CHK(mbedtls_mpi_mod_mpi(X, X, N));
1689    }
1690    MBEDTLS_MPI_CHK(mbedtls_mpi_grow(X, N->n));
1691
1692    /*
1693     * Convert to and from Montgomery around mbedtls_mpi_core_exp_mod().
1694     */
1695    {
1696        mbedtls_mpi_uint mm = mbedtls_mpi_core_montmul_init(N->p);
1697        mbedtls_mpi_core_to_mont_rep(X->p, X->p, N->p, N->n, mm, RR.p, T);
1698        mbedtls_mpi_core_exp_mod(X->p, X->p, N->p, N->n, E->p, E->n, RR.p, T);
1699        mbedtls_mpi_core_from_mont_rep(X->p, X->p, N->p, N->n, mm, T);
1700    }
1701
1702    /*
1703     * Correct for negative A.
1704     */
1705    if (A->s == -1 && (E->p[0] & 1) != 0) {
1706        mbedtls_ct_condition_t is_x_non_zero = mbedtls_mpi_core_check_zero_ct(X->p, X->n);
1707        X->s = mbedtls_ct_mpi_sign_if(is_x_non_zero, -1, 1);
1708
1709        MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(X, N, X));
1710    }
1711
1712cleanup:
1713
1714    mbedtls_mpi_zeroize_and_free(T, T_limbs);
1715
1716    if (prec_RR == NULL || prec_RR->p == NULL) {
1717        mbedtls_mpi_free(&RR);
1718    }
1719
1720    return ret;
1721}
1722
1723/*
1724 * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1725 */
1726int mbedtls_mpi_gcd(mbedtls_mpi *G, const mbedtls_mpi *A, const mbedtls_mpi *B)
1727{
1728    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1729    size_t lz, lzt;
1730    mbedtls_mpi TA, TB;
1731
1732    mbedtls_mpi_init(&TA); mbedtls_mpi_init(&TB);
1733
1734    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TA, A));
1735    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TB, B));
1736
1737    lz = mbedtls_mpi_lsb(&TA);
1738    lzt = mbedtls_mpi_lsb(&TB);
1739
1740    /* The loop below gives the correct result when A==0 but not when B==0.
1741     * So have a special case for B==0. Leverage the fact that we just
1742     * calculated the lsb and lsb(B)==0 iff B is odd or 0 to make the test
1743     * slightly more efficient than cmp_int(). */
1744    if (lzt == 0 && mbedtls_mpi_get_bit(&TB, 0) == 0) {
1745        ret = mbedtls_mpi_copy(G, A);
1746        goto cleanup;
1747    }
1748
1749    if (lzt < lz) {
1750        lz = lzt;
1751    }
1752
1753    TA.s = TB.s = 1;
1754
1755    /* We mostly follow the procedure described in HAC 14.54, but with some
1756     * minor differences:
1757     * - Sequences of multiplications or divisions by 2 are grouped into a
1758     *   single shift operation.
1759     * - The procedure in HAC assumes that 0 < TB <= TA.
1760     *     - The condition TB <= TA is not actually necessary for correctness.
1761     *       TA and TB have symmetric roles except for the loop termination
1762     *       condition, and the shifts at the beginning of the loop body
1763     *       remove any significance from the ordering of TA vs TB before
1764     *       the shifts.
1765     *     - If TA = 0, the loop goes through 0 iterations and the result is
1766     *       correctly TB.
1767     *     - The case TB = 0 was short-circuited above.
1768     *
1769     * For the correctness proof below, decompose the original values of
1770     * A and B as
1771     *   A = sa * 2^a * A' with A'=0 or A' odd, and sa = +-1
1772     *   B = sb * 2^b * B' with B'=0 or B' odd, and sb = +-1
1773     * Then gcd(A, B) = 2^{min(a,b)} * gcd(A',B'),
1774     * and gcd(A',B') is odd or 0.
1775     *
1776     * At the beginning, we have TA = |A| and TB = |B| so gcd(A,B) = gcd(TA,TB).
1777     * The code maintains the following invariant:
1778     *     gcd(A,B) = 2^k * gcd(TA,TB) for some k   (I)
1779     */
1780
1781    /* Proof that the loop terminates:
1782     * At each iteration, either the right-shift by 1 is made on a nonzero
1783     * value and the nonnegative integer bitlen(TA) + bitlen(TB) decreases
1784     * by at least 1, or the right-shift by 1 is made on zero and then
1785     * TA becomes 0 which ends the loop (TB cannot be 0 if it is right-shifted
1786     * since in that case TB is calculated from TB-TA with the condition TB>TA).
1787     */
1788    while (mbedtls_mpi_cmp_int(&TA, 0) != 0) {
1789        /* Divisions by 2 preserve the invariant (I). */
1790        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TA, mbedtls_mpi_lsb(&TA)));
1791        MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TB, mbedtls_mpi_lsb(&TB)));
1792
1793        /* Set either TA or TB to |TA-TB|/2. Since TA and TB are both odd,
1794         * TA-TB is even so the division by 2 has an integer result.
1795         * Invariant (I) is preserved since any odd divisor of both TA and TB
1796         * also divides |TA-TB|/2, and any odd divisor of both TA and |TA-TB|/2
1797         * also divides TB, and any odd divisor of both TB and |TA-TB|/2 also
1798         * divides TA.
1799         */
1800        if (mbedtls_mpi_cmp_mpi(&TA, &TB) >= 0) {
1801            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(&TA, &TA, &TB));
1802            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TA, 1));
1803        } else {
1804            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_abs(&TB, &TB, &TA));
1805            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TB, 1));
1806        }
1807        /* Note that one of TA or TB is still odd. */
1808    }
1809
1810    /* By invariant (I), gcd(A,B) = 2^k * gcd(TA,TB) for some k.
1811     * At the loop exit, TA = 0, so gcd(TA,TB) = TB.
1812     * - If there was at least one loop iteration, then one of TA or TB is odd,
1813     *   and TA = 0, so TB is odd and gcd(TA,TB) = gcd(A',B'). In this case,
1814     *   lz = min(a,b) so gcd(A,B) = 2^lz * TB.
1815     * - If there was no loop iteration, then A was 0, and gcd(A,B) = B.
1816     *   In this case, lz = 0 and B = TB so gcd(A,B) = B = 2^lz * TB as well.
1817     */
1818
1819    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_l(&TB, lz));
1820    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(G, &TB));
1821
1822cleanup:
1823
1824    mbedtls_mpi_free(&TA); mbedtls_mpi_free(&TB);
1825
1826    return ret;
1827}
1828
1829/*
1830 * Fill X with size bytes of random.
1831 * The bytes returned from the RNG are used in a specific order which
1832 * is suitable for deterministic ECDSA (see the specification of
1833 * mbedtls_mpi_random() and the implementation in mbedtls_mpi_fill_random()).
1834 */
1835int mbedtls_mpi_fill_random(mbedtls_mpi *X, size_t size,
1836                            int (*f_rng)(void *, unsigned char *, size_t),
1837                            void *p_rng)
1838{
1839    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1840    const size_t limbs = CHARS_TO_LIMBS(size);
1841
1842    /* Ensure that target MPI has exactly the necessary number of limbs */
1843    MBEDTLS_MPI_CHK(mbedtls_mpi_resize_clear(X, limbs));
1844    if (size == 0) {
1845        return 0;
1846    }
1847
1848    ret = mbedtls_mpi_core_fill_random(X->p, X->n, size, f_rng, p_rng);
1849
1850cleanup:
1851    return ret;
1852}
1853
1854int mbedtls_mpi_random(mbedtls_mpi *X,
1855                       mbedtls_mpi_sint min,
1856                       const mbedtls_mpi *N,
1857                       int (*f_rng)(void *, unsigned char *, size_t),
1858                       void *p_rng)
1859{
1860    if (min < 0) {
1861        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1862    }
1863    if (mbedtls_mpi_cmp_int(N, min) <= 0) {
1864        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1865    }
1866
1867    /* Ensure that target MPI has exactly the same number of limbs
1868     * as the upper bound, even if the upper bound has leading zeros.
1869     * This is necessary for mbedtls_mpi_core_random. */
1870    int ret = mbedtls_mpi_resize_clear(X, N->n);
1871    if (ret != 0) {
1872        return ret;
1873    }
1874
1875    return mbedtls_mpi_core_random(X->p, min, N->p, X->n, f_rng, p_rng);
1876}
1877
1878/*
1879 * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1880 */
1881int mbedtls_mpi_inv_mod(mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *N)
1882{
1883    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
1884    mbedtls_mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1885
1886    if (mbedtls_mpi_cmp_int(N, 1) <= 0) {
1887        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1888    }
1889
1890    mbedtls_mpi_init(&TA); mbedtls_mpi_init(&TU); mbedtls_mpi_init(&U1); mbedtls_mpi_init(&U2);
1891    mbedtls_mpi_init(&G); mbedtls_mpi_init(&TB); mbedtls_mpi_init(&TV);
1892    mbedtls_mpi_init(&V1); mbedtls_mpi_init(&V2);
1893
1894    MBEDTLS_MPI_CHK(mbedtls_mpi_gcd(&G, A, N));
1895
1896    if (mbedtls_mpi_cmp_int(&G, 1) != 0) {
1897        ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
1898        goto cleanup;
1899    }
1900
1901    MBEDTLS_MPI_CHK(mbedtls_mpi_mod_mpi(&TA, A, N));
1902    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TU, &TA));
1903    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TB, N));
1904    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&TV, N));
1905
1906    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&U1, 1));
1907    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&U2, 0));
1908    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&V1, 0));
1909    MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&V2, 1));
1910
1911    do {
1912        while ((TU.p[0] & 1) == 0) {
1913            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TU, 1));
1914
1915            if ((U1.p[0] & 1) != 0 || (U2.p[0] & 1) != 0) {
1916                MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(&U1, &U1, &TB));
1917                MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&U2, &U2, &TA));
1918            }
1919
1920            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&U1, 1));
1921            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&U2, 1));
1922        }
1923
1924        while ((TV.p[0] & 1) == 0) {
1925            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&TV, 1));
1926
1927            if ((V1.p[0] & 1) != 0 || (V2.p[0] & 1) != 0) {
1928                MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(&V1, &V1, &TB));
1929                MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&V2, &V2, &TA));
1930            }
1931
1932            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&V1, 1));
1933            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&V2, 1));
1934        }
1935
1936        if (mbedtls_mpi_cmp_mpi(&TU, &TV) >= 0) {
1937            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&TU, &TU, &TV));
1938            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&U1, &U1, &V1));
1939            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&U2, &U2, &V2));
1940        } else {
1941            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&TV, &TV, &TU));
1942            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&V1, &V1, &U1));
1943            MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&V2, &V2, &U2));
1944        }
1945    } while (mbedtls_mpi_cmp_int(&TU, 0) != 0);
1946
1947    while (mbedtls_mpi_cmp_int(&V1, 0) < 0) {
1948        MBEDTLS_MPI_CHK(mbedtls_mpi_add_mpi(&V1, &V1, N));
1949    }
1950
1951    while (mbedtls_mpi_cmp_mpi(&V1, N) >= 0) {
1952        MBEDTLS_MPI_CHK(mbedtls_mpi_sub_mpi(&V1, &V1, N));
1953    }
1954
1955    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(X, &V1));
1956
1957cleanup:
1958
1959    mbedtls_mpi_free(&TA); mbedtls_mpi_free(&TU); mbedtls_mpi_free(&U1); mbedtls_mpi_free(&U2);
1960    mbedtls_mpi_free(&G); mbedtls_mpi_free(&TB); mbedtls_mpi_free(&TV);
1961    mbedtls_mpi_free(&V1); mbedtls_mpi_free(&V2);
1962
1963    return ret;
1964}
1965
1966#if defined(MBEDTLS_GENPRIME)
1967
1968/* Gaps between primes, starting at 3. https://oeis.org/A001223 */
1969static const unsigned char small_prime_gaps[] = {
1970    2, 2, 4, 2, 4, 2, 4, 6,
1971    2, 6, 4, 2, 4, 6, 6, 2,
1972    6, 4, 2, 6, 4, 6, 8, 4,
1973    2, 4, 2, 4, 14, 4, 6, 2,
1974    10, 2, 6, 6, 4, 6, 6, 2,
1975    10, 2, 4, 2, 12, 12, 4, 2,
1976    4, 6, 2, 10, 6, 6, 6, 2,
1977    6, 4, 2, 10, 14, 4, 2, 4,
1978    14, 6, 10, 2, 4, 6, 8, 6,
1979    6, 4, 6, 8, 4, 8, 10, 2,
1980    10, 2, 6, 4, 6, 8, 4, 2,
1981    4, 12, 8, 4, 8, 4, 6, 12,
1982    2, 18, 6, 10, 6, 6, 2, 6,
1983    10, 6, 6, 2, 6, 6, 4, 2,
1984    12, 10, 2, 4, 6, 6, 2, 12,
1985    4, 6, 8, 10, 8, 10, 8, 6,
1986    6, 4, 8, 6, 4, 8, 4, 14,
1987    10, 12, 2, 10, 2, 4, 2, 10,
1988    14, 4, 2, 4, 14, 4, 2, 4,
1989    20, 4, 8, 10, 8, 4, 6, 6,
1990    14, 4, 6, 6, 8, 6, /*reaches 997*/
1991    0 /* the last entry is effectively unused */
1992};
1993
1994/*
1995 * Small divisors test (X must be positive)
1996 *
1997 * Return values:
1998 * 0: no small factor (possible prime, more tests needed)
1999 * 1: certain prime
2000 * MBEDTLS_ERR_MPI_NOT_ACCEPTABLE: certain non-prime
2001 * other negative: error
2002 */
2003static int mpi_check_small_factors(const mbedtls_mpi *X)
2004{
2005    int ret = 0;
2006    size_t i;
2007    mbedtls_mpi_uint r;
2008    unsigned p = 3; /* The first odd prime */
2009
2010    if ((X->p[0] & 1) == 0) {
2011        return MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2012    }
2013
2014    for (i = 0; i < sizeof(small_prime_gaps); p += small_prime_gaps[i], i++) {
2015        MBEDTLS_MPI_CHK(mbedtls_mpi_mod_int(&r, X, p));
2016        if (r == 0) {
2017            if (mbedtls_mpi_cmp_int(X, p) == 0) {
2018                return 1;
2019            } else {
2020                return MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2021            }
2022        }
2023    }
2024
2025cleanup:
2026    return ret;
2027}
2028
2029/*
2030 * Miller-Rabin pseudo-primality test  (HAC 4.24)
2031 */
2032static int mpi_miller_rabin(const mbedtls_mpi *X, size_t rounds,
2033                            int (*f_rng)(void *, unsigned char *, size_t),
2034                            void *p_rng)
2035{
2036    int ret, count;
2037    size_t i, j, k, s;
2038    mbedtls_mpi W, R, T, A, RR;
2039
2040    mbedtls_mpi_init(&W); mbedtls_mpi_init(&R);
2041    mbedtls_mpi_init(&T); mbedtls_mpi_init(&A);
2042    mbedtls_mpi_init(&RR);
2043
2044    /*
2045     * W = |X| - 1
2046     * R = W >> lsb( W )
2047     */
2048    MBEDTLS_MPI_CHK(mbedtls_mpi_sub_int(&W, X, 1));
2049    s = mbedtls_mpi_lsb(&W);
2050    MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&R, &W));
2051    MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&R, s));
2052
2053    for (i = 0; i < rounds; i++) {
2054        /*
2055         * pick a random A, 1 < A < |X| - 1
2056         */
2057        count = 0;
2058        do {
2059            MBEDTLS_MPI_CHK(mbedtls_mpi_fill_random(&A, X->n * ciL, f_rng, p_rng));
2060
2061            j = mbedtls_mpi_bitlen(&A);
2062            k = mbedtls_mpi_bitlen(&W);
2063            if (j > k) {
2064                A.p[A.n - 1] &= ((mbedtls_mpi_uint) 1 << (k - (A.n - 1) * biL - 1)) - 1;
2065            }
2066
2067            if (count++ > 30) {
2068                ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2069                goto cleanup;
2070            }
2071
2072        } while (mbedtls_mpi_cmp_mpi(&A, &W) >= 0 ||
2073                 mbedtls_mpi_cmp_int(&A, 1)  <= 0);
2074
2075        /*
2076         * A = A^R mod |X|
2077         */
2078        MBEDTLS_MPI_CHK(mbedtls_mpi_exp_mod(&A, &A, &R, X, &RR));
2079
2080        if (mbedtls_mpi_cmp_mpi(&A, &W) == 0 ||
2081            mbedtls_mpi_cmp_int(&A,  1) == 0) {
2082            continue;
2083        }
2084
2085        j = 1;
2086        while (j < s && mbedtls_mpi_cmp_mpi(&A, &W) != 0) {
2087            /*
2088             * A = A * A mod |X|
2089             */
2090            MBEDTLS_MPI_CHK(mbedtls_mpi_mul_mpi(&T, &A, &A));
2091            MBEDTLS_MPI_CHK(mbedtls_mpi_mod_mpi(&A, &T, X));
2092
2093            if (mbedtls_mpi_cmp_int(&A, 1) == 0) {
2094                break;
2095            }
2096
2097            j++;
2098        }
2099
2100        /*
2101         * not prime if A != |X| - 1 or A == 1
2102         */
2103        if (mbedtls_mpi_cmp_mpi(&A, &W) != 0 ||
2104            mbedtls_mpi_cmp_int(&A,  1) == 0) {
2105            ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2106            break;
2107        }
2108    }
2109
2110cleanup:
2111    mbedtls_mpi_free(&W); mbedtls_mpi_free(&R);
2112    mbedtls_mpi_free(&T); mbedtls_mpi_free(&A);
2113    mbedtls_mpi_free(&RR);
2114
2115    return ret;
2116}
2117
2118/*
2119 * Pseudo-primality test: small factors, then Miller-Rabin
2120 */
2121int mbedtls_mpi_is_prime_ext(const mbedtls_mpi *X, int rounds,
2122                             int (*f_rng)(void *, unsigned char *, size_t),
2123                             void *p_rng)
2124{
2125    int ret = MBEDTLS_ERR_ERROR_CORRUPTION_DETECTED;
2126    mbedtls_mpi XX;
2127
2128    XX.s = 1;
2129    XX.n = X->n;
2130    XX.p = X->p;
2131
2132    if (mbedtls_mpi_cmp_int(&XX, 0) == 0 ||
2133        mbedtls_mpi_cmp_int(&XX, 1) == 0) {
2134        return MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2135    }
2136
2137    if (mbedtls_mpi_cmp_int(&XX, 2) == 0) {
2138        return 0;
2139    }
2140
2141    if ((ret = mpi_check_small_factors(&XX)) != 0) {
2142        if (ret == 1) {
2143            return 0;
2144        }
2145
2146        return ret;
2147    }
2148
2149    return mpi_miller_rabin(&XX, rounds, f_rng, p_rng);
2150}
2151
2152/*
2153 * Prime number generation
2154 *
2155 * To generate an RSA key in a way recommended by FIPS 186-4, both primes must
2156 * be either 1024 bits or 1536 bits long, and flags must contain
2157 * MBEDTLS_MPI_GEN_PRIME_FLAG_LOW_ERR.
2158 */
2159int mbedtls_mpi_gen_prime(mbedtls_mpi *X, size_t nbits, int flags,
2160                          int (*f_rng)(void *, unsigned char *, size_t),
2161                          void *p_rng)
2162{
2163#ifdef MBEDTLS_HAVE_INT64
2164// ceil(2^63.5)
2165#define CEIL_MAXUINT_DIV_SQRT2 0xb504f333f9de6485ULL
2166#else
2167// ceil(2^31.5)
2168#define CEIL_MAXUINT_DIV_SQRT2 0xb504f334U
2169#endif
2170    int ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2171    size_t k, n;
2172    int rounds;
2173    mbedtls_mpi_uint r;
2174    mbedtls_mpi Y;
2175
2176    if (nbits < 3 || nbits > MBEDTLS_MPI_MAX_BITS) {
2177        return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
2178    }
2179
2180    mbedtls_mpi_init(&Y);
2181
2182    n = BITS_TO_LIMBS(nbits);
2183
2184    if ((flags & MBEDTLS_MPI_GEN_PRIME_FLAG_LOW_ERR) == 0) {
2185        /*
2186         * 2^-80 error probability, number of rounds chosen per HAC, table 4.4
2187         */
2188        rounds = ((nbits >= 1300) ?  2 : (nbits >=  850) ?  3 :
2189                  (nbits >=  650) ?  4 : (nbits >=  350) ?  8 :
2190                  (nbits >=  250) ? 12 : (nbits >=  150) ? 18 : 27);
2191    } else {
2192        /*
2193         * 2^-100 error probability, number of rounds computed based on HAC,
2194         * fact 4.48
2195         */
2196        rounds = ((nbits >= 1450) ?  4 : (nbits >=  1150) ?  5 :
2197                  (nbits >= 1000) ?  6 : (nbits >=   850) ?  7 :
2198                  (nbits >=  750) ?  8 : (nbits >=   500) ? 13 :
2199                  (nbits >=  250) ? 28 : (nbits >=   150) ? 40 : 51);
2200    }
2201
2202    while (1) {
2203        MBEDTLS_MPI_CHK(mbedtls_mpi_fill_random(X, n * ciL, f_rng, p_rng));
2204        /* make sure generated number is at least (nbits-1)+0.5 bits (FIPS 186-4 �B.3.3 steps 4.4, 5.5) */
2205        if (X->p[n-1] < CEIL_MAXUINT_DIV_SQRT2) {
2206            continue;
2207        }
2208
2209        k = n * biL;
2210        if (k > nbits) {
2211            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(X, k - nbits));
2212        }
2213        X->p[0] |= 1;
2214
2215        if ((flags & MBEDTLS_MPI_GEN_PRIME_FLAG_DH) == 0) {
2216            ret = mbedtls_mpi_is_prime_ext(X, rounds, f_rng, p_rng);
2217
2218            if (ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE) {
2219                goto cleanup;
2220            }
2221        } else {
2222            /*
2223             * A necessary condition for Y and X = 2Y + 1 to be prime
2224             * is X = 2 mod 3 (which is equivalent to Y = 2 mod 3).
2225             * Make sure it is satisfied, while keeping X = 3 mod 4
2226             */
2227
2228            X->p[0] |= 2;
2229
2230            MBEDTLS_MPI_CHK(mbedtls_mpi_mod_int(&r, X, 3));
2231            if (r == 0) {
2232                MBEDTLS_MPI_CHK(mbedtls_mpi_add_int(X, X, 8));
2233            } else if (r == 1) {
2234                MBEDTLS_MPI_CHK(mbedtls_mpi_add_int(X, X, 4));
2235            }
2236
2237            /* Set Y = (X-1) / 2, which is X / 2 because X is odd */
2238            MBEDTLS_MPI_CHK(mbedtls_mpi_copy(&Y, X));
2239            MBEDTLS_MPI_CHK(mbedtls_mpi_shift_r(&Y, 1));
2240
2241            while (1) {
2242                /*
2243                 * First, check small factors for X and Y
2244                 * before doing Miller-Rabin on any of them
2245                 */
2246                if ((ret = mpi_check_small_factors(X)) == 0 &&
2247                    (ret = mpi_check_small_factors(&Y)) == 0 &&
2248                    (ret = mpi_miller_rabin(X, rounds, f_rng, p_rng))
2249                    == 0 &&
2250                    (ret = mpi_miller_rabin(&Y, rounds, f_rng, p_rng))
2251                    == 0) {
2252                    goto cleanup;
2253                }
2254
2255                if (ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE) {
2256                    goto cleanup;
2257                }
2258
2259                /*
2260                 * Next candidates. We want to preserve Y = (X-1) / 2 and
2261                 * Y = 1 mod 2 and Y = 2 mod 3 (eq X = 3 mod 4 and X = 2 mod 3)
2262                 * so up Y by 6 and X by 12.
2263                 */
2264                MBEDTLS_MPI_CHK(mbedtls_mpi_add_int(X,  X, 12));
2265                MBEDTLS_MPI_CHK(mbedtls_mpi_add_int(&Y, &Y, 6));
2266            }
2267        }
2268    }
2269
2270cleanup:
2271
2272    mbedtls_mpi_free(&Y);
2273
2274    return ret;
2275}
2276
2277#endif /* MBEDTLS_GENPRIME */
2278
2279#if defined(MBEDTLS_SELF_TEST)
2280
2281#define GCD_PAIR_COUNT  3
2282
2283static const int gcd_pairs[GCD_PAIR_COUNT][3] =
2284{
2285    { 693, 609, 21 },
2286    { 1764, 868, 28 },
2287    { 768454923, 542167814, 1 }
2288};
2289
2290/*
2291 * Checkup routine
2292 */
2293int mbedtls_mpi_self_test(int verbose)
2294{
2295    int ret, i;
2296    mbedtls_mpi A, E, N, X, Y, U, V;
2297
2298    mbedtls_mpi_init(&A); mbedtls_mpi_init(&E); mbedtls_mpi_init(&N); mbedtls_mpi_init(&X);
2299    mbedtls_mpi_init(&Y); mbedtls_mpi_init(&U); mbedtls_mpi_init(&V);
2300
2301    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&A, 16,
2302                                            "EFE021C2645FD1DC586E69184AF4A31E" \
2303                                            "D5F53E93B5F123FA41680867BA110131" \
2304                                            "944FE7952E2517337780CB0DB80E61AA" \
2305                                            "E7C8DDC6C5C6AADEB34EB38A2F40D5E6"));
2306
2307    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&E, 16,
2308                                            "B2E7EFD37075B9F03FF989C7C5051C20" \
2309                                            "34D2A323810251127E7BF8625A4F49A5" \
2310                                            "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2311                                            "5B5C25763222FEFCCFC38B832366C29E"));
2312
2313    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&N, 16,
2314                                            "0066A198186C18C10B2F5ED9B522752A" \
2315                                            "9830B69916E535C8F047518A889A43A5" \
2316                                            "94B6BED27A168D31D4A52F88925AA8F5"));
2317
2318    MBEDTLS_MPI_CHK(mbedtls_mpi_mul_mpi(&X, &A, &N));
2319
2320    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&U, 16,
2321                                            "602AB7ECA597A3D6B56FF9829A5E8B85" \
2322                                            "9E857EA95A03512E2BAE7391688D264A" \
2323                                            "A5663B0341DB9CCFD2C4C5F421FEC814" \
2324                                            "8001B72E848A38CAE1C65F78E56ABDEF" \
2325                                            "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2326                                            "ECF677152EF804370C1A305CAF3B5BF1" \
2327                                            "30879B56C61DE584A0F53A2447A51E"));
2328
2329    if (verbose != 0) {
2330        mbedtls_printf("  MPI test #1 (mul_mpi): ");
2331    }
2332
2333    if (mbedtls_mpi_cmp_mpi(&X, &U) != 0) {
2334        if (verbose != 0) {
2335            mbedtls_printf("failed\n");
2336        }
2337
2338        ret = 1;
2339        goto cleanup;
2340    }
2341
2342    if (verbose != 0) {
2343        mbedtls_printf("passed\n");
2344    }
2345
2346    MBEDTLS_MPI_CHK(mbedtls_mpi_div_mpi(&X, &Y, &A, &N));
2347
2348    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&U, 16,
2349                                            "256567336059E52CAE22925474705F39A94"));
2350
2351    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&V, 16,
2352                                            "6613F26162223DF488E9CD48CC132C7A" \
2353                                            "0AC93C701B001B092E4E5B9F73BCD27B" \
2354                                            "9EE50D0657C77F374E903CDFA4C642"));
2355
2356    if (verbose != 0) {
2357        mbedtls_printf("  MPI test #2 (div_mpi): ");
2358    }
2359
2360    if (mbedtls_mpi_cmp_mpi(&X, &U) != 0 ||
2361        mbedtls_mpi_cmp_mpi(&Y, &V) != 0) {
2362        if (verbose != 0) {
2363            mbedtls_printf("failed\n");
2364        }
2365
2366        ret = 1;
2367        goto cleanup;
2368    }
2369
2370    if (verbose != 0) {
2371        mbedtls_printf("passed\n");
2372    }
2373
2374    MBEDTLS_MPI_CHK(mbedtls_mpi_exp_mod(&X, &A, &E, &N, NULL));
2375
2376    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&U, 16,
2377                                            "36E139AEA55215609D2816998ED020BB" \
2378                                            "BD96C37890F65171D948E9BC7CBAA4D9" \
2379                                            "325D24D6A3C12710F10A09FA08AB87"));
2380
2381    if (verbose != 0) {
2382        mbedtls_printf("  MPI test #3 (exp_mod): ");
2383    }
2384
2385    if (mbedtls_mpi_cmp_mpi(&X, &U) != 0) {
2386        if (verbose != 0) {
2387            mbedtls_printf("failed\n");
2388        }
2389
2390        ret = 1;
2391        goto cleanup;
2392    }
2393
2394    if (verbose != 0) {
2395        mbedtls_printf("passed\n");
2396    }
2397
2398    MBEDTLS_MPI_CHK(mbedtls_mpi_inv_mod(&X, &A, &N));
2399
2400    MBEDTLS_MPI_CHK(mbedtls_mpi_read_string(&U, 16,
2401                                            "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2402                                            "C3DBA76456363A10869622EAC2DD84EC" \
2403                                            "C5B8A74DAC4D09E03B5E0BE779F2DF61"));
2404
2405    if (verbose != 0) {
2406        mbedtls_printf("  MPI test #4 (inv_mod): ");
2407    }
2408
2409    if (mbedtls_mpi_cmp_mpi(&X, &U) != 0) {
2410        if (verbose != 0) {
2411            mbedtls_printf("failed\n");
2412        }
2413
2414        ret = 1;
2415        goto cleanup;
2416    }
2417
2418    if (verbose != 0) {
2419        mbedtls_printf("passed\n");
2420    }
2421
2422    if (verbose != 0) {
2423        mbedtls_printf("  MPI test #5 (simple gcd): ");
2424    }
2425
2426    for (i = 0; i < GCD_PAIR_COUNT; i++) {
2427        MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&X, gcd_pairs[i][0]));
2428        MBEDTLS_MPI_CHK(mbedtls_mpi_lset(&Y, gcd_pairs[i][1]));
2429
2430        MBEDTLS_MPI_CHK(mbedtls_mpi_gcd(&A, &X, &Y));
2431
2432        if (mbedtls_mpi_cmp_int(&A, gcd_pairs[i][2]) != 0) {
2433            if (verbose != 0) {
2434                mbedtls_printf("failed at %d\n", i);
2435            }
2436
2437            ret = 1;
2438            goto cleanup;
2439        }
2440    }
2441
2442    if (verbose != 0) {
2443        mbedtls_printf("passed\n");
2444    }
2445
2446cleanup:
2447
2448    if (ret != 0 && verbose != 0) {
2449        mbedtls_printf("Unexpected error, return code = %08X\n", (unsigned int) ret);
2450    }
2451
2452    mbedtls_mpi_free(&A); mbedtls_mpi_free(&E); mbedtls_mpi_free(&N); mbedtls_mpi_free(&X);
2453    mbedtls_mpi_free(&Y); mbedtls_mpi_free(&U); mbedtls_mpi_free(&V);
2454
2455    if (verbose != 0) {
2456        mbedtls_printf("\n");
2457    }
2458
2459    return ret;
2460}
2461
2462#endif /* MBEDTLS_SELF_TEST */
2463
2464#endif /* MBEDTLS_BIGNUM_C */
2465