xref: /reactos/dll/3rdparty/mbedtls/bignum.c (revision 34593d93)
1 /*
2  *  Multi-precision integer library
3  *
4  *  Copyright (C) 2006-2015, ARM Limited, All Rights Reserved
5  *  SPDX-License-Identifier: Apache-2.0 OR GPL-2.0-or-later
6  *
7  *  This file is provided under the Apache License 2.0, or the
8  *  GNU General Public License v2.0 or later.
9  *
10  *  **********
11  *  Apache License 2.0:
12  *
13  *  Licensed under the Apache License, Version 2.0 (the "License"); you may
14  *  not use this file except in compliance with the License.
15  *  You may obtain a copy of the License at
16  *
17  *  http://www.apache.org/licenses/LICENSE-2.0
18  *
19  *  Unless required by applicable law or agreed to in writing, software
20  *  distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
21  *  WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
22  *  See the License for the specific language governing permissions and
23  *  limitations under the License.
24  *
25  *  **********
26  *
27  *  **********
28  *  GNU General Public License v2.0 or later:
29  *
30  *  This program is free software; you can redistribute it and/or modify
31  *  it under the terms of the GNU General Public License as published by
32  *  the Free Software Foundation; either version 2 of the License, or
33  *  (at your option) any later version.
34  *
35  *  This program is distributed in the hope that it will be useful,
36  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
37  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
38  *  GNU General Public License for more details.
39  *
40  *  You should have received a copy of the GNU General Public License along
41  *  with this program; if not, write to the Free Software Foundation, Inc.,
42  *  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
43  *
44  *  **********
45  *
46  *  This file is part of mbed TLS (https://tls.mbed.org)
47  */
48 
49 /*
50  *  The following sources were referenced in the design of this Multi-precision
51  *  Integer library:
52  *
53  *  [1] Handbook of Applied Cryptography - 1997
54  *      Menezes, van Oorschot and Vanstone
55  *
56  *  [2] Multi-Precision Math
57  *      Tom St Denis
58  *      https://github.com/libtom/libtommath/blob/develop/tommath.pdf
59  *
60  *  [3] GNU Multi-Precision Arithmetic Library
61  *      https://gmplib.org/manual/index.html
62  *
63  */
64 
65 #if !defined(MBEDTLS_CONFIG_FILE)
66 #include "mbedtls/config.h"
67 #else
68 #include MBEDTLS_CONFIG_FILE
69 #endif
70 
71 #if defined(MBEDTLS_BIGNUM_C)
72 
73 #include "mbedtls/bignum.h"
74 #include "mbedtls/bn_mul.h"
75 
76 #include <string.h>
77 
78 #if defined(MBEDTLS_PLATFORM_C)
79 #include "mbedtls/platform.h"
80 #else
81 #include <stdio.h>
82 #include <stdlib.h>
83 #define mbedtls_printf     printf
84 #define mbedtls_calloc    calloc
85 #define mbedtls_free       free
86 #endif
87 
88 /* Implementation that should never be optimized out by the compiler */
89 static void mbedtls_mpi_zeroize( mbedtls_mpi_uint *v, size_t n ) {
90     volatile mbedtls_mpi_uint *p = v; while( n-- ) *p++ = 0;
91 }
92 
93 /* Implementation that should never be optimized out by the compiler */
94 static void mbedtls_zeroize( void *v, size_t n ) {
95     volatile unsigned char *p = v; while( n-- ) *p++ = 0;
96 }
97 
98 #define ciL    (sizeof(mbedtls_mpi_uint))         /* chars in limb  */
99 #define biL    (ciL << 3)               /* bits  in limb  */
100 #define biH    (ciL << 2)               /* half limb size */
101 
102 #define MPI_SIZE_T_MAX  ( (size_t) -1 ) /* SIZE_T_MAX is not standard */
103 
104 /*
105  * Convert between bits/chars and number of limbs
106  * Divide first in order to avoid potential overflows
107  */
108 #define BITS_TO_LIMBS(i)  ( (i) / biL + ( (i) % biL != 0 ) )
109 #define CHARS_TO_LIMBS(i) ( (i) / ciL + ( (i) % ciL != 0 ) )
110 
111 /*
112  * Initialize one MPI
113  */
114 void mbedtls_mpi_init( mbedtls_mpi *X )
115 {
116     if( X == NULL )
117         return;
118 
119     X->s = 1;
120     X->n = 0;
121     X->p = NULL;
122 }
123 
124 /*
125  * Unallocate one MPI
126  */
127 void mbedtls_mpi_free( mbedtls_mpi *X )
128 {
129     if( X == NULL )
130         return;
131 
132     if( X->p != NULL )
133     {
134         mbedtls_mpi_zeroize( X->p, X->n );
135         mbedtls_free( X->p );
136     }
137 
138     X->s = 1;
139     X->n = 0;
140     X->p = NULL;
141 }
142 
143 /*
144  * Enlarge to the specified number of limbs
145  */
146 int mbedtls_mpi_grow( mbedtls_mpi *X, size_t nblimbs )
147 {
148     mbedtls_mpi_uint *p;
149 
150     if( nblimbs > MBEDTLS_MPI_MAX_LIMBS )
151         return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
152 
153     if( X->n < nblimbs )
154     {
155         if( ( p = (mbedtls_mpi_uint*)mbedtls_calloc( nblimbs, ciL ) ) == NULL )
156             return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
157 
158         if( X->p != NULL )
159         {
160             memcpy( p, X->p, X->n * ciL );
161             mbedtls_mpi_zeroize( X->p, X->n );
162             mbedtls_free( X->p );
163         }
164 
165         X->n = nblimbs;
166         X->p = p;
167     }
168 
169     return( 0 );
170 }
171 
172 /*
173  * Resize down as much as possible,
174  * while keeping at least the specified number of limbs
175  */
176 int mbedtls_mpi_shrink( mbedtls_mpi *X, size_t nblimbs )
177 {
178     mbedtls_mpi_uint *p;
179     size_t i;
180 
181     /* Actually resize up if there are currently fewer than nblimbs limbs. */
182     if( X->n <= nblimbs )
183         return( mbedtls_mpi_grow( X, nblimbs ) );
184     /* After this point, then X->n > nblimbs and in particular X->n > 0. */
185 
186     for( i = X->n - 1; i > 0; i-- )
187         if( X->p[i] != 0 )
188             break;
189     i++;
190 
191     if( i < nblimbs )
192         i = nblimbs;
193 
194     if( ( p = (mbedtls_mpi_uint*)mbedtls_calloc( i, ciL ) ) == NULL )
195         return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
196 
197     if( X->p != NULL )
198     {
199         memcpy( p, X->p, i * ciL );
200         mbedtls_mpi_zeroize( X->p, X->n );
201         mbedtls_free( X->p );
202     }
203 
204     X->n = i;
205     X->p = p;
206 
207     return( 0 );
208 }
209 
210 /*
211  * Copy the contents of Y into X
212  */
213 int mbedtls_mpi_copy( mbedtls_mpi *X, const mbedtls_mpi *Y )
214 {
215     int ret;
216     size_t i;
217 
218     if( X == Y )
219         return( 0 );
220 
221     if( Y->n == 0 )
222     {
223         mbedtls_mpi_free( X );
224         return( 0 );
225     }
226 
227     for( i = Y->n - 1; i > 0; i-- )
228         if( Y->p[i] != 0 )
229             break;
230     i++;
231 
232     X->s = Y->s;
233 
234     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i ) );
235 
236     memset( X->p, 0, X->n * ciL );
237     memcpy( X->p, Y->p, i * ciL );
238 
239 cleanup:
240 
241     return( ret );
242 }
243 
244 /*
245  * Swap the contents of X and Y
246  */
247 void mbedtls_mpi_swap( mbedtls_mpi *X, mbedtls_mpi *Y )
248 {
249     mbedtls_mpi T;
250 
251     memcpy( &T,  X, sizeof( mbedtls_mpi ) );
252     memcpy(  X,  Y, sizeof( mbedtls_mpi ) );
253     memcpy(  Y, &T, sizeof( mbedtls_mpi ) );
254 }
255 
256 /*
257  * Conditionally assign dest = src, without leaking information
258  * about whether the assignment was made or not.
259  * dest and src must be arrays of limbs of size n.
260  * assign must be 0 or 1.
261  */
262 static void mpi_safe_cond_assign( size_t n,
263                                   mbedtls_mpi_uint *dest,
264                                   const mbedtls_mpi_uint *src,
265                                   unsigned char assign )
266 {
267     size_t i;
268     for( i = 0; i < n; i++ )
269         dest[i] = dest[i] * ( 1 - assign ) + src[i] * assign;
270 }
271 
272 /*
273  * Conditionally assign X = Y, without leaking information
274  * about whether the assignment was made or not.
275  * (Leaking information about the respective sizes of X and Y is ok however.)
276  */
277 int mbedtls_mpi_safe_cond_assign( mbedtls_mpi *X, const mbedtls_mpi *Y, unsigned char assign )
278 {
279     int ret = 0;
280     size_t i;
281 
282     /* make sure assign is 0 or 1 in a time-constant manner */
283     assign = (assign | (unsigned char)-assign) >> 7;
284 
285     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, Y->n ) );
286 
287     X->s = X->s * ( 1 - assign ) + Y->s * assign;
288 
289     mpi_safe_cond_assign( Y->n, X->p, Y->p, assign );
290 
291     for( i = Y->n; i < X->n; i++ )
292         X->p[i] *= ( 1 - assign );
293 
294 cleanup:
295     return( ret );
296 }
297 
298 /*
299  * Conditionally swap X and Y, without leaking information
300  * about whether the swap was made or not.
301  * Here it is not ok to simply swap the pointers, which whould lead to
302  * different memory access patterns when X and Y are used afterwards.
303  */
304 int mbedtls_mpi_safe_cond_swap( mbedtls_mpi *X, mbedtls_mpi *Y, unsigned char swap )
305 {
306     int ret, s;
307     size_t i;
308     mbedtls_mpi_uint tmp;
309 
310     if( X == Y )
311         return( 0 );
312 
313     /* make sure swap is 0 or 1 in a time-constant manner */
314     swap = (swap | (unsigned char)-swap) >> 7;
315 
316     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, Y->n ) );
317     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( Y, X->n ) );
318 
319     s = X->s;
320     X->s = X->s * ( 1 - swap ) + Y->s * swap;
321     Y->s = Y->s * ( 1 - swap ) +    s * swap;
322 
323 
324     for( i = 0; i < X->n; i++ )
325     {
326         tmp = X->p[i];
327         X->p[i] = X->p[i] * ( 1 - swap ) + Y->p[i] * swap;
328         Y->p[i] = Y->p[i] * ( 1 - swap ) +     tmp * swap;
329     }
330 
331 cleanup:
332     return( ret );
333 }
334 
335 /*
336  * Set value from integer
337  */
338 int mbedtls_mpi_lset( mbedtls_mpi *X, mbedtls_mpi_sint z )
339 {
340     int ret;
341 
342     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, 1 ) );
343     memset( X->p, 0, X->n * ciL );
344 
345     X->p[0] = ( z < 0 ) ? -z : z;
346     X->s    = ( z < 0 ) ? -1 : 1;
347 
348 cleanup:
349 
350     return( ret );
351 }
352 
353 /*
354  * Get a specific bit
355  */
356 int mbedtls_mpi_get_bit( const mbedtls_mpi *X, size_t pos )
357 {
358     if( X->n * biL <= pos )
359         return( 0 );
360 
361     return( ( X->p[pos / biL] >> ( pos % biL ) ) & 0x01 );
362 }
363 
364 /* Get a specific byte, without range checks. */
365 #define GET_BYTE( X, i )                                \
366     ( ( ( X )->p[( i ) / ciL] >> ( ( ( i ) % ciL ) * 8 ) ) & 0xff )
367 
368 /*
369  * Set a bit to a specific value of 0 or 1
370  */
371 int mbedtls_mpi_set_bit( mbedtls_mpi *X, size_t pos, unsigned char val )
372 {
373     int ret = 0;
374     size_t off = pos / biL;
375     size_t idx = pos % biL;
376 
377     if( val != 0 && val != 1 )
378         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
379 
380     if( X->n * biL <= pos )
381     {
382         if( val == 0 )
383             return( 0 );
384 
385         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, off + 1 ) );
386     }
387 
388     X->p[off] &= ~( (mbedtls_mpi_uint) 0x01 << idx );
389     X->p[off] |= (mbedtls_mpi_uint) val << idx;
390 
391 cleanup:
392 
393     return( ret );
394 }
395 
396 /*
397  * Return the number of less significant zero-bits
398  */
399 size_t mbedtls_mpi_lsb( const mbedtls_mpi *X )
400 {
401     size_t i, j, count = 0;
402 
403     for( i = 0; i < X->n; i++ )
404         for( j = 0; j < biL; j++, count++ )
405             if( ( ( X->p[i] >> j ) & 1 ) != 0 )
406                 return( count );
407 
408     return( 0 );
409 }
410 
411 /*
412  * Count leading zero bits in a given integer
413  */
414 static size_t mbedtls_clz( const mbedtls_mpi_uint x )
415 {
416     size_t j;
417     mbedtls_mpi_uint mask = (mbedtls_mpi_uint) 1 << (biL - 1);
418 
419     for( j = 0; j < biL; j++ )
420     {
421         if( x & mask ) break;
422 
423         mask >>= 1;
424     }
425 
426     return j;
427 }
428 
429 /*
430  * Return the number of bits
431  */
432 size_t mbedtls_mpi_bitlen( const mbedtls_mpi *X )
433 {
434     size_t i, j;
435 
436     if( X->n == 0 )
437         return( 0 );
438 
439     for( i = X->n - 1; i > 0; i-- )
440         if( X->p[i] != 0 )
441             break;
442 
443     j = biL - mbedtls_clz( X->p[i] );
444 
445     return( ( i * biL ) + j );
446 }
447 
448 /*
449  * Return the total size in bytes
450  */
451 size_t mbedtls_mpi_size( const mbedtls_mpi *X )
452 {
453     return( ( mbedtls_mpi_bitlen( X ) + 7 ) >> 3 );
454 }
455 
456 /*
457  * Convert an ASCII character to digit value
458  */
459 static int mpi_get_digit( mbedtls_mpi_uint *d, int radix, char c )
460 {
461     *d = 255;
462 
463     if( c >= 0x30 && c <= 0x39 ) *d = c - 0x30;
464     if( c >= 0x41 && c <= 0x46 ) *d = c - 0x37;
465     if( c >= 0x61 && c <= 0x66 ) *d = c - 0x57;
466 
467     if( *d >= (mbedtls_mpi_uint) radix )
468         return( MBEDTLS_ERR_MPI_INVALID_CHARACTER );
469 
470     return( 0 );
471 }
472 
473 /*
474  * Import from an ASCII string
475  */
476 int mbedtls_mpi_read_string( mbedtls_mpi *X, int radix, const char *s )
477 {
478     int ret;
479     size_t i, j, slen, n;
480     mbedtls_mpi_uint d;
481     mbedtls_mpi T;
482 
483     if( radix < 2 || radix > 16 )
484         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
485 
486     mbedtls_mpi_init( &T );
487 
488     slen = strlen( s );
489 
490     if( radix == 16 )
491     {
492         if( slen > MPI_SIZE_T_MAX >> 2 )
493             return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
494 
495         n = BITS_TO_LIMBS( slen << 2 );
496 
497         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, n ) );
498         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
499 
500         for( i = slen, j = 0; i > 0; i--, j++ )
501         {
502             if( i == 1 && s[i - 1] == '-' )
503             {
504                 X->s = -1;
505                 break;
506             }
507 
508             MBEDTLS_MPI_CHK( mpi_get_digit( &d, radix, s[i - 1] ) );
509             X->p[j / ( 2 * ciL )] |= d << ( ( j % ( 2 * ciL ) ) << 2 );
510         }
511     }
512     else
513     {
514         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
515 
516         for( i = 0; i < slen; i++ )
517         {
518             if( i == 0 && s[i] == '-' )
519             {
520                 X->s = -1;
521                 continue;
522             }
523 
524             MBEDTLS_MPI_CHK( mpi_get_digit( &d, radix, s[i] ) );
525             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T, X, radix ) );
526 
527             if( X->s == 1 )
528             {
529                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, &T, d ) );
530             }
531             else
532             {
533                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_int( X, &T, d ) );
534             }
535         }
536     }
537 
538 cleanup:
539 
540     mbedtls_mpi_free( &T );
541 
542     return( ret );
543 }
544 
545 /*
546  * Helper to write the digits high-order first.
547  */
548 static int mpi_write_hlp( mbedtls_mpi *X, int radix,
549                           char **p, const size_t buflen )
550 {
551     int ret;
552     mbedtls_mpi_uint r;
553     size_t length = 0;
554     char *p_end = *p + buflen;
555 
556     do
557     {
558         if( length >= buflen )
559         {
560             return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
561         }
562 
563         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, radix ) );
564         MBEDTLS_MPI_CHK( mbedtls_mpi_div_int( X, NULL, X, radix ) );
565         /*
566          * Write the residue in the current position, as an ASCII character.
567          */
568         if( r < 0xA )
569             *(--p_end) = (char)( '0' + r );
570         else
571             *(--p_end) = (char)( 'A' + ( r - 0xA ) );
572 
573         length++;
574     } while( mbedtls_mpi_cmp_int( X, 0 ) != 0 );
575 
576     memmove( *p, p_end, length );
577     *p += length;
578 
579 cleanup:
580 
581     return( ret );
582 }
583 
584 /*
585  * Export into an ASCII string
586  */
587 int mbedtls_mpi_write_string( const mbedtls_mpi *X, int radix,
588                               char *buf, size_t buflen, size_t *olen )
589 {
590     int ret = 0;
591     size_t n;
592     char *p;
593     mbedtls_mpi T;
594 
595     if( radix < 2 || radix > 16 )
596         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
597 
598     n = mbedtls_mpi_bitlen( X ); /* Number of bits necessary to present `n`. */
599     if( radix >=  4 ) n >>= 1;   /* Number of 4-adic digits necessary to present
600                                   * `n`. If radix > 4, this might be a strict
601                                   * overapproximation of the number of
602                                   * radix-adic digits needed to present `n`. */
603     if( radix >= 16 ) n >>= 1;   /* Number of hexadecimal digits necessary to
604                                   * present `n`. */
605 
606     n += 1; /* Terminating null byte */
607     n += 1; /* Compensate for the divisions above, which round down `n`
608              * in case it's not even. */
609     n += 1; /* Potential '-'-sign. */
610     n += ( n & 1 ); /* Make n even to have enough space for hexadecimal writing,
611                      * which always uses an even number of hex-digits. */
612 
613     if( buflen < n )
614     {
615         *olen = n;
616         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
617     }
618 
619     p = buf;
620     mbedtls_mpi_init( &T );
621 
622     if( X->s == -1 )
623     {
624         *p++ = '-';
625         buflen--;
626     }
627 
628     if( radix == 16 )
629     {
630         int c;
631         size_t i, j, k;
632 
633         for( i = X->n, k = 0; i > 0; i-- )
634         {
635             for( j = ciL; j > 0; j-- )
636             {
637                 c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
638 
639                 if( c == 0 && k == 0 && ( i + j ) != 2 )
640                     continue;
641 
642                 *(p++) = "0123456789ABCDEF" [c / 16];
643                 *(p++) = "0123456789ABCDEF" [c % 16];
644                 k = 1;
645             }
646         }
647     }
648     else
649     {
650         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T, X ) );
651 
652         if( T.s == -1 )
653             T.s = 1;
654 
655         MBEDTLS_MPI_CHK( mpi_write_hlp( &T, radix, &p, buflen ) );
656     }
657 
658     *p++ = '\0';
659     *olen = p - buf;
660 
661 cleanup:
662 
663     mbedtls_mpi_free( &T );
664 
665     return( ret );
666 }
667 
668 #if defined(MBEDTLS_FS_IO)
669 /*
670  * Read X from an opened file
671  */
672 int mbedtls_mpi_read_file( mbedtls_mpi *X, int radix, FILE *fin )
673 {
674     mbedtls_mpi_uint d;
675     size_t slen;
676     char *p;
677     /*
678      * Buffer should have space for (short) label and decimal formatted MPI,
679      * newline characters and '\0'
680      */
681     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
682 
683     memset( s, 0, sizeof( s ) );
684     if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
685         return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
686 
687     slen = strlen( s );
688     if( slen == sizeof( s ) - 2 )
689         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
690 
691     if( slen > 0 && s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
692     if( slen > 0 && s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }
693 
694     p = s + slen;
695     while( p-- > s )
696         if( mpi_get_digit( &d, radix, *p ) != 0 )
697             break;
698 
699     return( mbedtls_mpi_read_string( X, radix, p + 1 ) );
700 }
701 
702 /*
703  * Write X into an opened file (or stdout if fout == NULL)
704  */
705 int mbedtls_mpi_write_file( const char *p, const mbedtls_mpi *X, int radix, FILE *fout )
706 {
707     int ret;
708     size_t n, slen, plen;
709     /*
710      * Buffer should have space for (short) label and decimal formatted MPI,
711      * newline characters and '\0'
712      */
713     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
714 
715     memset( s, 0, sizeof( s ) );
716 
717     MBEDTLS_MPI_CHK( mbedtls_mpi_write_string( X, radix, s, sizeof( s ) - 2, &n ) );
718 
719     if( p == NULL ) p = "";
720 
721     plen = strlen( p );
722     slen = strlen( s );
723     s[slen++] = '\r';
724     s[slen++] = '\n';
725 
726     if( fout != NULL )
727     {
728         if( fwrite( p, 1, plen, fout ) != plen ||
729             fwrite( s, 1, slen, fout ) != slen )
730             return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
731     }
732     else
733         mbedtls_printf( "%s%s", p, s );
734 
735 cleanup:
736 
737     return( ret );
738 }
739 #endif /* MBEDTLS_FS_IO */
740 
741 /*
742  * Import X from unsigned binary data, big endian
743  */
744 int mbedtls_mpi_read_binary( mbedtls_mpi *X, const unsigned char *buf, size_t buflen )
745 {
746     int ret;
747     size_t i, j;
748     size_t const limbs = CHARS_TO_LIMBS( buflen );
749 
750     /* Ensure that target MPI has exactly the necessary number of limbs */
751     if( X->n != limbs )
752     {
753         mbedtls_mpi_free( X );
754         mbedtls_mpi_init( X );
755         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, limbs ) );
756     }
757 
758     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
759 
760     for( i = buflen, j = 0; i > 0; i--, j++ )
761         X->p[j / ciL] |= ((mbedtls_mpi_uint) buf[i - 1]) << ((j % ciL) << 3);
762 
763 cleanup:
764 
765     return( ret );
766 }
767 
768 /*
769  * Export X into unsigned binary data, big endian
770  */
771 int mbedtls_mpi_write_binary( const mbedtls_mpi *X,
772                               unsigned char *buf, size_t buflen )
773 {
774     size_t stored_bytes = X->n * ciL;
775     size_t bytes_to_copy;
776     unsigned char *p;
777     size_t i;
778 
779     if( stored_bytes < buflen )
780     {
781         /* There is enough space in the output buffer. Write initial
782          * null bytes and record the position at which to start
783          * writing the significant bytes. In this case, the execution
784          * trace of this function does not depend on the value of the
785          * number. */
786         bytes_to_copy = stored_bytes;
787         p = buf + buflen - stored_bytes;
788         memset( buf, 0, buflen - stored_bytes );
789     }
790     else
791     {
792         /* The output buffer is smaller than the allocated size of X.
793          * However X may fit if its leading bytes are zero. */
794         bytes_to_copy = buflen;
795         p = buf;
796         for( i = bytes_to_copy; i < stored_bytes; i++ )
797         {
798             if( GET_BYTE( X, i ) != 0 )
799                 return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
800         }
801     }
802 
803     for( i = 0; i < bytes_to_copy; i++ )
804         p[bytes_to_copy - i - 1] = GET_BYTE( X, i );
805 
806     return( 0 );
807 }
808 
809 /*
810  * Left-shift: X <<= count
811  */
812 int mbedtls_mpi_shift_l( mbedtls_mpi *X, size_t count )
813 {
814     int ret;
815     size_t i, v0, t1;
816     mbedtls_mpi_uint r0 = 0, r1;
817 
818     v0 = count / (biL    );
819     t1 = count & (biL - 1);
820 
821     i = mbedtls_mpi_bitlen( X ) + count;
822 
823     if( X->n * biL < i )
824         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, BITS_TO_LIMBS( i ) ) );
825 
826     ret = 0;
827 
828     /*
829      * shift by count / limb_size
830      */
831     if( v0 > 0 )
832     {
833         for( i = X->n; i > v0; i-- )
834             X->p[i - 1] = X->p[i - v0 - 1];
835 
836         for( ; i > 0; i-- )
837             X->p[i - 1] = 0;
838     }
839 
840     /*
841      * shift by count % limb_size
842      */
843     if( t1 > 0 )
844     {
845         for( i = v0; i < X->n; i++ )
846         {
847             r1 = X->p[i] >> (biL - t1);
848             X->p[i] <<= t1;
849             X->p[i] |= r0;
850             r0 = r1;
851         }
852     }
853 
854 cleanup:
855 
856     return( ret );
857 }
858 
859 /*
860  * Right-shift: X >>= count
861  */
862 int mbedtls_mpi_shift_r( mbedtls_mpi *X, size_t count )
863 {
864     size_t i, v0, v1;
865     mbedtls_mpi_uint r0 = 0, r1;
866 
867     v0 = count /  biL;
868     v1 = count & (biL - 1);
869 
870     if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
871         return mbedtls_mpi_lset( X, 0 );
872 
873     /*
874      * shift by count / limb_size
875      */
876     if( v0 > 0 )
877     {
878         for( i = 0; i < X->n - v0; i++ )
879             X->p[i] = X->p[i + v0];
880 
881         for( ; i < X->n; i++ )
882             X->p[i] = 0;
883     }
884 
885     /*
886      * shift by count % limb_size
887      */
888     if( v1 > 0 )
889     {
890         for( i = X->n; i > 0; i-- )
891         {
892             r1 = X->p[i - 1] << (biL - v1);
893             X->p[i - 1] >>= v1;
894             X->p[i - 1] |= r0;
895             r0 = r1;
896         }
897     }
898 
899     return( 0 );
900 }
901 
902 /*
903  * Compare unsigned values
904  */
905 int mbedtls_mpi_cmp_abs( const mbedtls_mpi *X, const mbedtls_mpi *Y )
906 {
907     size_t i, j;
908 
909     for( i = X->n; i > 0; i-- )
910         if( X->p[i - 1] != 0 )
911             break;
912 
913     for( j = Y->n; j > 0; j-- )
914         if( Y->p[j - 1] != 0 )
915             break;
916 
917     if( i == 0 && j == 0 )
918         return( 0 );
919 
920     if( i > j ) return(  1 );
921     if( j > i ) return( -1 );
922 
923     for( ; i > 0; i-- )
924     {
925         if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
926         if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
927     }
928 
929     return( 0 );
930 }
931 
932 /*
933  * Compare signed values
934  */
935 int mbedtls_mpi_cmp_mpi( const mbedtls_mpi *X, const mbedtls_mpi *Y )
936 {
937     size_t i, j;
938 
939     for( i = X->n; i > 0; i-- )
940         if( X->p[i - 1] != 0 )
941             break;
942 
943     for( j = Y->n; j > 0; j-- )
944         if( Y->p[j - 1] != 0 )
945             break;
946 
947     if( i == 0 && j == 0 )
948         return( 0 );
949 
950     if( i > j ) return(  X->s );
951     if( j > i ) return( -Y->s );
952 
953     if( X->s > 0 && Y->s < 0 ) return(  1 );
954     if( Y->s > 0 && X->s < 0 ) return( -1 );
955 
956     for( ; i > 0; i-- )
957     {
958         if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
959         if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
960     }
961 
962     return( 0 );
963 }
964 
965 /** Decide if an integer is less than the other, without branches.
966  *
967  * \param x         First integer.
968  * \param y         Second integer.
969  *
970  * \return          1 if \p x is less than \p y, 0 otherwise
971  */
972 static unsigned ct_lt_mpi_uint( const mbedtls_mpi_uint x,
973         const mbedtls_mpi_uint y )
974 {
975     mbedtls_mpi_uint ret;
976     mbedtls_mpi_uint cond;
977 
978     /*
979      * Check if the most significant bits (MSB) of the operands are different.
980      */
981     cond = ( x ^ y );
982     /*
983      * If the MSB are the same then the difference x-y will be negative (and
984      * have its MSB set to 1 during conversion to unsigned) if and only if x<y.
985      */
986     ret = ( x - y ) & ~cond;
987     /*
988      * If the MSB are different, then the operand with the MSB of 1 is the
989      * bigger. (That is if y has MSB of 1, then x<y is true and it is false if
990      * the MSB of y is 0.)
991      */
992     ret |= y & cond;
993 
994 
995     ret = ret >> ( biL - 1 );
996 
997     return (unsigned) ret;
998 }
999 
1000 /*
1001  * Compare signed values in constant time
1002  */
1003 int mbedtls_mpi_lt_mpi_ct( const mbedtls_mpi *X, const mbedtls_mpi *Y,
1004         unsigned *ret )
1005 {
1006     size_t i;
1007     /* The value of any of these variables is either 0 or 1 at all times. */
1008     unsigned cond, done, X_is_negative, Y_is_negative;
1009 
1010     if( X->n != Y->n )
1011         return MBEDTLS_ERR_MPI_BAD_INPUT_DATA;
1012 
1013     /*
1014      * Set sign_N to 1 if N >= 0, 0 if N < 0.
1015      * We know that N->s == 1 if N >= 0 and N->s == -1 if N < 0.
1016      */
1017     X_is_negative = ( X->s & 2 ) >> 1;
1018     Y_is_negative = ( Y->s & 2 ) >> 1;
1019 
1020     /*
1021      * If the signs are different, then the positive operand is the bigger.
1022      * That is if X is negative (X_is_negative == 1), then X < Y is true and it
1023      * is false if X is positive (X_is_negative == 0).
1024      */
1025     cond = ( X_is_negative ^ Y_is_negative );
1026     *ret = cond & X_is_negative;
1027 
1028     /*
1029      * This is a constant-time function. We might have the result, but we still
1030      * need to go through the loop. Record if we have the result already.
1031      */
1032     done = cond;
1033 
1034     for( i = X->n; i > 0; i-- )
1035     {
1036         /*
1037          * If Y->p[i - 1] < X->p[i - 1] then X < Y is true if and only if both
1038          * X and Y are negative.
1039          *
1040          * Again even if we can make a decision, we just mark the result and
1041          * the fact that we are done and continue looping.
1042          */
1043         cond = ct_lt_mpi_uint( Y->p[i - 1], X->p[i - 1] );
1044         *ret |= cond & ( 1 - done ) & X_is_negative;
1045         done |= cond;
1046 
1047         /*
1048          * If X->p[i - 1] < Y->p[i - 1] then X < Y is true if and only if both
1049          * X and Y are positive.
1050          *
1051          * Again even if we can make a decision, we just mark the result and
1052          * the fact that we are done and continue looping.
1053          */
1054         cond = ct_lt_mpi_uint( X->p[i - 1], Y->p[i - 1] );
1055         *ret |= cond & ( 1 - done ) & ( 1 - X_is_negative );
1056         done |= cond;
1057     }
1058 
1059     return( 0 );
1060 }
1061 
1062 /*
1063  * Compare signed values
1064  */
1065 int mbedtls_mpi_cmp_int( const mbedtls_mpi *X, mbedtls_mpi_sint z )
1066 {
1067     mbedtls_mpi Y;
1068     mbedtls_mpi_uint p[1];
1069 
1070     *p  = ( z < 0 ) ? -z : z;
1071     Y.s = ( z < 0 ) ? -1 : 1;
1072     Y.n = 1;
1073     Y.p = p;
1074 
1075     return( mbedtls_mpi_cmp_mpi( X, &Y ) );
1076 }
1077 
1078 /*
1079  * Unsigned addition: X = |A| + |B|  (HAC 14.7)
1080  */
1081 int mbedtls_mpi_add_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1082 {
1083     int ret;
1084     size_t i, j;
1085     mbedtls_mpi_uint *o, *p, c, tmp;
1086 
1087     if( X == B )
1088     {
1089         const mbedtls_mpi *T = A; A = X; B = T;
1090     }
1091 
1092     if( X != A )
1093         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
1094 
1095     /*
1096      * X should always be positive as a result of unsigned additions.
1097      */
1098     X->s = 1;
1099 
1100     for( j = B->n; j > 0; j-- )
1101         if( B->p[j - 1] != 0 )
1102             break;
1103 
1104     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
1105 
1106     o = B->p; p = X->p; c = 0;
1107 
1108     /*
1109      * tmp is used because it might happen that p == o
1110      */
1111     for( i = 0; i < j; i++, o++, p++ )
1112     {
1113         tmp= *o;
1114         *p +=  c; c  = ( *p <  c );
1115         *p += tmp; c += ( *p < tmp );
1116     }
1117 
1118     while( c != 0 )
1119     {
1120         if( i >= X->n )
1121         {
1122             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + 1 ) );
1123             p = X->p + i;
1124         }
1125 
1126         *p += c; c = ( *p < c ); i++; p++;
1127     }
1128 
1129 cleanup:
1130 
1131     return( ret );
1132 }
1133 
1134 /**
1135  * Helper for mbedtls_mpi subtraction.
1136  *
1137  * Calculate d - s where d and s have the same size.
1138  * This function operates modulo (2^ciL)^n and returns the carry
1139  * (1 if there was a wraparound, i.e. if `d < s`, and 0 otherwise).
1140  *
1141  * \param n             Number of limbs of \p d and \p s.
1142  * \param[in,out] d     On input, the left operand.
1143  *                      On output, the result of the subtraction:
1144  * \param[in] s         The right operand.
1145  *
1146  * \return              1 if `d < s`.
1147  *                      0 if `d >= s`.
1148  */
1149 static mbedtls_mpi_uint mpi_sub_hlp( size_t n,
1150                                      mbedtls_mpi_uint *d,
1151                                      const mbedtls_mpi_uint *s )
1152 {
1153     size_t i;
1154     mbedtls_mpi_uint c, z;
1155 
1156     for( i = c = 0; i < n; i++, s++, d++ )
1157     {
1158         z = ( *d <  c );     *d -=  c;
1159         c = ( *d < *s ) + z; *d -= *s;
1160     }
1161 
1162     return( c );
1163 }
1164 
1165 /*
1166  * Unsigned subtraction: X = |A| - |B|  (HAC 14.9, 14.10)
1167  */
1168 int mbedtls_mpi_sub_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1169 {
1170     mbedtls_mpi TB;
1171     int ret;
1172     size_t n;
1173     mbedtls_mpi_uint carry;
1174 
1175     mbedtls_mpi_init( &TB );
1176 
1177     if( X == B )
1178     {
1179         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
1180         B = &TB;
1181     }
1182 
1183     if( X != A )
1184         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
1185 
1186     /*
1187      * X should always be positive as a result of unsigned subtractions.
1188      */
1189     X->s = 1;
1190 
1191     ret = 0;
1192 
1193     for( n = B->n; n > 0; n-- )
1194         if( B->p[n - 1] != 0 )
1195             break;
1196 
1197     carry = mpi_sub_hlp( n, X->p, B->p );
1198     if( carry != 0 )
1199     {
1200         /* Propagate the carry to the first nonzero limb of X. */
1201         for( ; n < X->n && X->p[n] == 0; n++ )
1202             --X->p[n];
1203         /* If we ran out of space for the carry, it means that the result
1204          * is negative. */
1205         if( n == X->n )
1206             return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1207         --X->p[n];
1208     }
1209 
1210 cleanup:
1211 
1212     mbedtls_mpi_free( &TB );
1213 
1214     return( ret );
1215 }
1216 
1217 /*
1218  * Signed addition: X = A + B
1219  */
1220 int mbedtls_mpi_add_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1221 {
1222     int ret, s = A->s;
1223 
1224     if( A->s * B->s < 0 )
1225     {
1226         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1227         {
1228             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1229             X->s =  s;
1230         }
1231         else
1232         {
1233             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1234             X->s = -s;
1235         }
1236     }
1237     else
1238     {
1239         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1240         X->s = s;
1241     }
1242 
1243 cleanup:
1244 
1245     return( ret );
1246 }
1247 
1248 /*
1249  * Signed subtraction: X = A - B
1250  */
1251 int mbedtls_mpi_sub_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1252 {
1253     int ret, s = A->s;
1254 
1255     if( A->s * B->s > 0 )
1256     {
1257         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1258         {
1259             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1260             X->s =  s;
1261         }
1262         else
1263         {
1264             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1265             X->s = -s;
1266         }
1267     }
1268     else
1269     {
1270         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1271         X->s = s;
1272     }
1273 
1274 cleanup:
1275 
1276     return( ret );
1277 }
1278 
1279 /*
1280  * Signed addition: X = A + b
1281  */
1282 int mbedtls_mpi_add_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1283 {
1284     mbedtls_mpi _B;
1285     mbedtls_mpi_uint p[1];
1286 
1287     p[0] = ( b < 0 ) ? -b : b;
1288     _B.s = ( b < 0 ) ? -1 : 1;
1289     _B.n = 1;
1290     _B.p = p;
1291 
1292     return( mbedtls_mpi_add_mpi( X, A, &_B ) );
1293 }
1294 
1295 /*
1296  * Signed subtraction: X = A - b
1297  */
1298 int mbedtls_mpi_sub_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1299 {
1300     mbedtls_mpi _B;
1301     mbedtls_mpi_uint p[1];
1302 
1303     p[0] = ( b < 0 ) ? -b : b;
1304     _B.s = ( b < 0 ) ? -1 : 1;
1305     _B.n = 1;
1306     _B.p = p;
1307 
1308     return( mbedtls_mpi_sub_mpi( X, A, &_B ) );
1309 }
1310 
1311 /*
1312  * Helper for mbedtls_mpi multiplication
1313  */
1314 static
1315 #if defined(__APPLE__) && defined(__arm__)
1316 /*
1317  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
1318  * appears to need this to prevent bad ARM code generation at -O3.
1319  */
1320 __attribute__ ((noinline))
1321 #endif
1322 void mpi_mul_hlp( size_t i, mbedtls_mpi_uint *s, mbedtls_mpi_uint *d, mbedtls_mpi_uint b )
1323 {
1324     mbedtls_mpi_uint c = 0, t = 0;
1325 
1326 #if defined(MULADDC_HUIT)
1327     for( ; i >= 8; i -= 8 )
1328     {
1329         MULADDC_INIT
1330         MULADDC_HUIT
1331         MULADDC_STOP
1332     }
1333 
1334     for( ; i > 0; i-- )
1335     {
1336         MULADDC_INIT
1337         MULADDC_CORE
1338         MULADDC_STOP
1339     }
1340 #else /* MULADDC_HUIT */
1341     for( ; i >= 16; i -= 16 )
1342     {
1343         MULADDC_INIT
1344         MULADDC_CORE   MULADDC_CORE
1345         MULADDC_CORE   MULADDC_CORE
1346         MULADDC_CORE   MULADDC_CORE
1347         MULADDC_CORE   MULADDC_CORE
1348 
1349         MULADDC_CORE   MULADDC_CORE
1350         MULADDC_CORE   MULADDC_CORE
1351         MULADDC_CORE   MULADDC_CORE
1352         MULADDC_CORE   MULADDC_CORE
1353         MULADDC_STOP
1354     }
1355 
1356     for( ; i >= 8; i -= 8 )
1357     {
1358         MULADDC_INIT
1359         MULADDC_CORE   MULADDC_CORE
1360         MULADDC_CORE   MULADDC_CORE
1361 
1362         MULADDC_CORE   MULADDC_CORE
1363         MULADDC_CORE   MULADDC_CORE
1364         MULADDC_STOP
1365     }
1366 
1367     for( ; i > 0; i-- )
1368     {
1369         MULADDC_INIT
1370         MULADDC_CORE
1371         MULADDC_STOP
1372     }
1373 #endif /* MULADDC_HUIT */
1374 
1375     t++;
1376 
1377     do {
1378         *d += c; c = ( *d < c ); d++;
1379     }
1380     while( c != 0 );
1381 }
1382 
1383 /*
1384  * Baseline multiplication: X = A * B  (HAC 14.12)
1385  */
1386 int mbedtls_mpi_mul_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1387 {
1388     int ret;
1389     size_t i, j;
1390     mbedtls_mpi TA, TB;
1391 
1392     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
1393 
1394     if( X == A ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) ); A = &TA; }
1395     if( X == B ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) ); B = &TB; }
1396 
1397     for( i = A->n; i > 0; i-- )
1398         if( A->p[i - 1] != 0 )
1399             break;
1400 
1401     for( j = B->n; j > 0; j-- )
1402         if( B->p[j - 1] != 0 )
1403             break;
1404 
1405     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + j ) );
1406     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
1407 
1408     for( i++; j > 0; j-- )
1409         mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1410 
1411     X->s = A->s * B->s;
1412 
1413 cleanup:
1414 
1415     mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TA );
1416 
1417     return( ret );
1418 }
1419 
1420 /*
1421  * Baseline multiplication: X = A * b
1422  */
1423 int mbedtls_mpi_mul_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_uint b )
1424 {
1425     mbedtls_mpi _B;
1426     mbedtls_mpi_uint p[1];
1427 
1428     _B.s = 1;
1429     _B.n = 1;
1430     _B.p = p;
1431     p[0] = b;
1432 
1433     return( mbedtls_mpi_mul_mpi( X, A, &_B ) );
1434 }
1435 
1436 /*
1437  * Unsigned integer divide - double mbedtls_mpi_uint dividend, u1/u0, and
1438  * mbedtls_mpi_uint divisor, d
1439  */
1440 static mbedtls_mpi_uint mbedtls_int_div_int( mbedtls_mpi_uint u1,
1441             mbedtls_mpi_uint u0, mbedtls_mpi_uint d, mbedtls_mpi_uint *r )
1442 {
1443 #if defined(MBEDTLS_HAVE_UDBL)
1444     mbedtls_t_udbl dividend, quotient;
1445 #else
1446     const mbedtls_mpi_uint radix = (mbedtls_mpi_uint) 1 << biH;
1447     const mbedtls_mpi_uint uint_halfword_mask = ( (mbedtls_mpi_uint) 1 << biH ) - 1;
1448     mbedtls_mpi_uint d0, d1, q0, q1, rAX, r0, quotient;
1449     mbedtls_mpi_uint u0_msw, u0_lsw;
1450     size_t s;
1451 #endif
1452 
1453     /*
1454      * Check for overflow
1455      */
1456     if( 0 == d || u1 >= d )
1457     {
1458         if (r != NULL) *r = ~0;
1459 
1460         return ( ~0 );
1461     }
1462 
1463 #if defined(MBEDTLS_HAVE_UDBL)
1464     dividend  = (mbedtls_t_udbl) u1 << biL;
1465     dividend |= (mbedtls_t_udbl) u0;
1466     quotient = dividend / d;
1467     if( quotient > ( (mbedtls_t_udbl) 1 << biL ) - 1 )
1468         quotient = ( (mbedtls_t_udbl) 1 << biL ) - 1;
1469 
1470     if( r != NULL )
1471         *r = (mbedtls_mpi_uint)( dividend - (quotient * d ) );
1472 
1473     return (mbedtls_mpi_uint) quotient;
1474 #else
1475 
1476     /*
1477      * Algorithm D, Section 4.3.1 - The Art of Computer Programming
1478      *   Vol. 2 - Seminumerical Algorithms, Knuth
1479      */
1480 
1481     /*
1482      * Normalize the divisor, d, and dividend, u0, u1
1483      */
1484     s = mbedtls_clz( d );
1485     d = d << s;
1486 
1487     u1 = u1 << s;
1488     u1 |= ( u0 >> ( biL - s ) ) & ( -(mbedtls_mpi_sint)s >> ( biL - 1 ) );
1489     u0 =  u0 << s;
1490 
1491     d1 = d >> biH;
1492     d0 = d & uint_halfword_mask;
1493 
1494     u0_msw = u0 >> biH;
1495     u0_lsw = u0 & uint_halfword_mask;
1496 
1497     /*
1498      * Find the first quotient and remainder
1499      */
1500     q1 = u1 / d1;
1501     r0 = u1 - d1 * q1;
1502 
1503     while( q1 >= radix || ( q1 * d0 > radix * r0 + u0_msw ) )
1504     {
1505         q1 -= 1;
1506         r0 += d1;
1507 
1508         if ( r0 >= radix ) break;
1509     }
1510 
1511     rAX = ( u1 * radix ) + ( u0_msw - q1 * d );
1512     q0 = rAX / d1;
1513     r0 = rAX - q0 * d1;
1514 
1515     while( q0 >= radix || ( q0 * d0 > radix * r0 + u0_lsw ) )
1516     {
1517         q0 -= 1;
1518         r0 += d1;
1519 
1520         if ( r0 >= radix ) break;
1521     }
1522 
1523     if (r != NULL)
1524         *r = ( rAX * radix + u0_lsw - q0 * d ) >> s;
1525 
1526     quotient = q1 * radix + q0;
1527 
1528     return quotient;
1529 #endif
1530 }
1531 
1532 /*
1533  * Division by mbedtls_mpi: A = Q * B + R  (HAC 14.20)
1534  */
1535 int mbedtls_mpi_div_mpi( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1536 {
1537     int ret;
1538     size_t i, n, t, k;
1539     mbedtls_mpi X, Y, Z, T1, T2;
1540 
1541     if( mbedtls_mpi_cmp_int( B, 0 ) == 0 )
1542         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1543 
1544     mbedtls_mpi_init( &X ); mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &Z );
1545     mbedtls_mpi_init( &T1 ); mbedtls_mpi_init( &T2 );
1546 
1547     if( mbedtls_mpi_cmp_abs( A, B ) < 0 )
1548     {
1549         if( Q != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_lset( Q, 0 ) );
1550         if( R != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, A ) );
1551         return( 0 );
1552     }
1553 
1554     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &X, A ) );
1555     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, B ) );
1556     X.s = Y.s = 1;
1557 
1558     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &Z, A->n + 2 ) );
1559     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Z,  0 ) );
1560     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T1, 2 ) );
1561     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T2, 3 ) );
1562 
1563     k = mbedtls_mpi_bitlen( &Y ) % biL;
1564     if( k < biL - 1 )
1565     {
1566         k = biL - 1 - k;
1567         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &X, k ) );
1568         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, k ) );
1569     }
1570     else k = 0;
1571 
1572     n = X.n - 1;
1573     t = Y.n - 1;
1574     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, biL * ( n - t ) ) );
1575 
1576     while( mbedtls_mpi_cmp_mpi( &X, &Y ) >= 0 )
1577     {
1578         Z.p[n - t]++;
1579         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &Y ) );
1580     }
1581     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, biL * ( n - t ) ) );
1582 
1583     for( i = n; i > t ; i-- )
1584     {
1585         if( X.p[i] >= Y.p[t] )
1586             Z.p[i - t - 1] = ~0;
1587         else
1588         {
1589             Z.p[i - t - 1] = mbedtls_int_div_int( X.p[i], X.p[i - 1],
1590                                                             Y.p[t], NULL);
1591         }
1592 
1593         Z.p[i - t - 1]++;
1594         do
1595         {
1596             Z.p[i - t - 1]--;
1597 
1598             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T1, 0 ) );
1599             T1.p[0] = ( t < 1 ) ? 0 : Y.p[t - 1];
1600             T1.p[1] = Y.p[t];
1601             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1602 
1603             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T2, 0 ) );
1604             T2.p[0] = ( i < 2 ) ? 0 : X.p[i - 2];
1605             T2.p[1] = ( i < 1 ) ? 0 : X.p[i - 1];
1606             T2.p[2] = X.p[i];
1607         }
1608         while( mbedtls_mpi_cmp_mpi( &T1, &T2 ) > 0 );
1609 
1610         MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1611         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1,  biL * ( i - t - 1 ) ) );
1612         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &T1 ) );
1613 
1614         if( mbedtls_mpi_cmp_int( &X, 0 ) < 0 )
1615         {
1616             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T1, &Y ) );
1617             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1, biL * ( i - t - 1 ) ) );
1618             MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &X, &X, &T1 ) );
1619             Z.p[i - t - 1]--;
1620         }
1621     }
1622 
1623     if( Q != NULL )
1624     {
1625         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( Q, &Z ) );
1626         Q->s = A->s * B->s;
1627     }
1628 
1629     if( R != NULL )
1630     {
1631         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &X, k ) );
1632         X.s = A->s;
1633         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, &X ) );
1634 
1635         if( mbedtls_mpi_cmp_int( R, 0 ) == 0 )
1636             R->s = 1;
1637     }
1638 
1639 cleanup:
1640 
1641     mbedtls_mpi_free( &X ); mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &Z );
1642     mbedtls_mpi_free( &T1 ); mbedtls_mpi_free( &T2 );
1643 
1644     return( ret );
1645 }
1646 
1647 /*
1648  * Division by int: A = Q * b + R
1649  */
1650 int mbedtls_mpi_div_int( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1651 {
1652     mbedtls_mpi _B;
1653     mbedtls_mpi_uint p[1];
1654 
1655     p[0] = ( b < 0 ) ? -b : b;
1656     _B.s = ( b < 0 ) ? -1 : 1;
1657     _B.n = 1;
1658     _B.p = p;
1659 
1660     return( mbedtls_mpi_div_mpi( Q, R, A, &_B ) );
1661 }
1662 
1663 /*
1664  * Modulo: R = A mod B
1665  */
1666 int mbedtls_mpi_mod_mpi( mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1667 {
1668     int ret;
1669 
1670     if( mbedtls_mpi_cmp_int( B, 0 ) < 0 )
1671         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1672 
1673     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( NULL, R, A, B ) );
1674 
1675     while( mbedtls_mpi_cmp_int( R, 0 ) < 0 )
1676       MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( R, R, B ) );
1677 
1678     while( mbedtls_mpi_cmp_mpi( R, B ) >= 0 )
1679       MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( R, R, B ) );
1680 
1681 cleanup:
1682 
1683     return( ret );
1684 }
1685 
1686 /*
1687  * Modulo: r = A mod b
1688  */
1689 int mbedtls_mpi_mod_int( mbedtls_mpi_uint *r, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1690 {
1691     size_t i;
1692     mbedtls_mpi_uint x, y, z;
1693 
1694     if( b == 0 )
1695         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1696 
1697     if( b < 0 )
1698         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1699 
1700     /*
1701      * handle trivial cases
1702      */
1703     if( b == 1 )
1704     {
1705         *r = 0;
1706         return( 0 );
1707     }
1708 
1709     if( b == 2 )
1710     {
1711         *r = A->p[0] & 1;
1712         return( 0 );
1713     }
1714 
1715     /*
1716      * general case
1717      */
1718     for( i = A->n, y = 0; i > 0; i-- )
1719     {
1720         x  = A->p[i - 1];
1721         y  = ( y << biH ) | ( x >> biH );
1722         z  = y / b;
1723         y -= z * b;
1724 
1725         x <<= biH;
1726         y  = ( y << biH ) | ( x >> biH );
1727         z  = y / b;
1728         y -= z * b;
1729     }
1730 
1731     /*
1732      * If A is negative, then the current y represents a negative value.
1733      * Flipping it to the positive side.
1734      */
1735     if( A->s < 0 && y != 0 )
1736         y = b - y;
1737 
1738     *r = y;
1739 
1740     return( 0 );
1741 }
1742 
1743 /*
1744  * Fast Montgomery initialization (thanks to Tom St Denis)
1745  */
1746 static void mpi_montg_init( mbedtls_mpi_uint *mm, const mbedtls_mpi *N )
1747 {
1748     mbedtls_mpi_uint x, m0 = N->p[0];
1749     unsigned int i;
1750 
1751     x  = m0;
1752     x += ( ( m0 + 2 ) & 4 ) << 1;
1753 
1754     for( i = biL; i >= 8; i /= 2 )
1755         x *= ( 2 - ( m0 * x ) );
1756 
1757     *mm = ~x + 1;
1758 }
1759 
1760 /** Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1761  *
1762  * \param[in,out]   A   One of the numbers to multiply.
1763  *                      It must have at least as many limbs as N
1764  *                      (A->n >= N->n), and any limbs beyond n are ignored.
1765  *                      On successful completion, A contains the result of
1766  *                      the multiplication A * B * R^-1 mod N where
1767  *                      R = (2^ciL)^n.
1768  * \param[in]       B   One of the numbers to multiply.
1769  *                      It must be nonzero and must not have more limbs than N
1770  *                      (B->n <= N->n).
1771  * \param[in]       N   The modulo. N must be odd.
1772  * \param           mm  The value calculated by `mpi_montg_init(&mm, N)`.
1773  *                      This is -N^-1 mod 2^ciL.
1774  * \param[in,out]   T   A bignum for temporary storage.
1775  *                      It must be at least twice the limb size of N plus 2
1776  *                      (T->n >= 2 * (N->n + 1)).
1777  *                      Its initial content is unused and
1778  *                      its final content is indeterminate.
1779  *                      Note that unlike the usual convention in the library
1780  *                      for `const mbedtls_mpi*`, the content of T can change.
1781  */
1782 static void mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,
1783                          const mbedtls_mpi *T )
1784 {
1785     size_t i, n, m;
1786     mbedtls_mpi_uint u0, u1, *d;
1787 
1788     memset( T->p, 0, T->n * ciL );
1789 
1790     d = T->p;
1791     n = N->n;
1792     m = ( B->n < n ) ? B->n : n;
1793 
1794     for( i = 0; i < n; i++ )
1795     {
1796         /*
1797          * T = (T + u0*B + u1*N) / 2^biL
1798          */
1799         u0 = A->p[i];
1800         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1801 
1802         mpi_mul_hlp( m, B->p, d, u0 );
1803         mpi_mul_hlp( n, N->p, d, u1 );
1804 
1805         *d++ = u0; d[n + 1] = 0;
1806     }
1807 
1808     /* At this point, d is either the desired result or the desired result
1809      * plus N. We now potentially subtract N, avoiding leaking whether the
1810      * subtraction is performed through side channels. */
1811 
1812     /* Copy the n least significant limbs of d to A, so that
1813      * A = d if d < N (recall that N has n limbs). */
1814     memcpy( A->p, d, n * ciL );
1815     /* If d >= N then we want to set A to d - N. To prevent timing attacks,
1816      * do the calculation without using conditional tests. */
1817     /* Set d to d0 + (2^biL)^n - N where d0 is the current value of d. */
1818     d[n] += 1;
1819     d[n] -= mpi_sub_hlp( n, d, N->p );
1820     /* If d0 < N then d < (2^biL)^n
1821      * so d[n] == 0 and we want to keep A as it is.
1822      * If d0 >= N then d >= (2^biL)^n, and d <= (2^biL)^n + N < 2 * (2^biL)^n
1823      * so d[n] == 1 and we want to set A to the result of the subtraction
1824      * which is d - (2^biL)^n, i.e. the n least significant limbs of d.
1825      * This exactly corresponds to a conditional assignment. */
1826     mpi_safe_cond_assign( n, A->p, d, (unsigned char) d[n] );
1827 }
1828 
1829 /*
1830  * Montgomery reduction: A = A * R^-1 mod N
1831  *
1832  * See mpi_montmul() regarding constraints and guarantees on the parameters.
1833  */
1834 static void mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N, mbedtls_mpi_uint mm, const mbedtls_mpi *T )
1835 {
1836     mbedtls_mpi_uint z = 1;
1837     mbedtls_mpi U;
1838 
1839     U.n = U.s = (int) z;
1840     U.p = &z;
1841 
1842     mpi_montmul( A, &U, N, mm, T );
1843 }
1844 
1845 /*
1846  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1847  */
1848 int mbedtls_mpi_exp_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *E, const mbedtls_mpi *N, mbedtls_mpi *_RR )
1849 {
1850     int ret;
1851     size_t wbits, wsize, one = 1;
1852     size_t i, j, nblimbs;
1853     size_t bufsize, nbits;
1854     mbedtls_mpi_uint ei, mm, state;
1855     mbedtls_mpi RR, T, W[ 2 << MBEDTLS_MPI_WINDOW_SIZE ], Apos;
1856     int neg;
1857 
1858     if( mbedtls_mpi_cmp_int( N, 0 ) <= 0 || ( N->p[0] & 1 ) == 0 )
1859         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1860 
1861     if( mbedtls_mpi_cmp_int( E, 0 ) < 0 )
1862         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1863 
1864     /*
1865      * Init temps and window size
1866      */
1867     mpi_montg_init( &mm, N );
1868     mbedtls_mpi_init( &RR ); mbedtls_mpi_init( &T );
1869     mbedtls_mpi_init( &Apos );
1870     memset( W, 0, sizeof( W ) );
1871 
1872     i = mbedtls_mpi_bitlen( E );
1873 
1874     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1875             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1876 
1877 #if( MBEDTLS_MPI_WINDOW_SIZE < 6 )
1878     if( wsize > MBEDTLS_MPI_WINDOW_SIZE )
1879         wsize = MBEDTLS_MPI_WINDOW_SIZE;
1880 #endif
1881 
1882     j = N->n + 1;
1883     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
1884     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[1],  j ) );
1885     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T, j * 2 ) );
1886 
1887     /*
1888      * Compensate for negative A (and correct at the end)
1889      */
1890     neg = ( A->s == -1 );
1891     if( neg )
1892     {
1893         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Apos, A ) );
1894         Apos.s = 1;
1895         A = &Apos;
1896     }
1897 
1898     /*
1899      * If 1st call, pre-compute R^2 mod N
1900      */
1901     if( _RR == NULL || _RR->p == NULL )
1902     {
1903         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &RR, 1 ) );
1904         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &RR, N->n * 2 * biL ) );
1905         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &RR, &RR, N ) );
1906 
1907         if( _RR != NULL )
1908             memcpy( _RR, &RR, sizeof( mbedtls_mpi ) );
1909     }
1910     else
1911         memcpy( &RR, _RR, sizeof( mbedtls_mpi ) );
1912 
1913     /*
1914      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1915      */
1916     if( mbedtls_mpi_cmp_mpi( A, N ) >= 0 )
1917         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &W[1], A, N ) );
1918     else
1919         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[1], A ) );
1920 
1921     mpi_montmul( &W[1], &RR, N, mm, &T );
1922 
1923     /*
1924      * X = R^2 * R^-1 mod N = R mod N
1925      */
1926     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &RR ) );
1927     mpi_montred( X, N, mm, &T );
1928 
1929     if( wsize > 1 )
1930     {
1931         /*
1932          * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1933          */
1934         j =  one << ( wsize - 1 );
1935 
1936         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[j], N->n + 1 ) );
1937         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[j], &W[1]    ) );
1938 
1939         for( i = 0; i < wsize - 1; i++ )
1940             mpi_montmul( &W[j], &W[j], N, mm, &T );
1941 
1942         /*
1943          * W[i] = W[i - 1] * W[1]
1944          */
1945         for( i = j + 1; i < ( one << wsize ); i++ )
1946         {
1947             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[i], N->n + 1 ) );
1948             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[i], &W[i - 1] ) );
1949 
1950             mpi_montmul( &W[i], &W[1], N, mm, &T );
1951         }
1952     }
1953 
1954     nblimbs = E->n;
1955     bufsize = 0;
1956     nbits   = 0;
1957     wbits   = 0;
1958     state   = 0;
1959 
1960     while( 1 )
1961     {
1962         if( bufsize == 0 )
1963         {
1964             if( nblimbs == 0 )
1965                 break;
1966 
1967             nblimbs--;
1968 
1969             bufsize = sizeof( mbedtls_mpi_uint ) << 3;
1970         }
1971 
1972         bufsize--;
1973 
1974         ei = (E->p[nblimbs] >> bufsize) & 1;
1975 
1976         /*
1977          * skip leading 0s
1978          */
1979         if( ei == 0 && state == 0 )
1980             continue;
1981 
1982         if( ei == 0 && state == 1 )
1983         {
1984             /*
1985              * out of window, square X
1986              */
1987             mpi_montmul( X, X, N, mm, &T );
1988             continue;
1989         }
1990 
1991         /*
1992          * add ei to current window
1993          */
1994         state = 2;
1995 
1996         nbits++;
1997         wbits |= ( ei << ( wsize - nbits ) );
1998 
1999         if( nbits == wsize )
2000         {
2001             /*
2002              * X = X^wsize R^-1 mod N
2003              */
2004             for( i = 0; i < wsize; i++ )
2005                 mpi_montmul( X, X, N, mm, &T );
2006 
2007             /*
2008              * X = X * W[wbits] R^-1 mod N
2009              */
2010             mpi_montmul( X, &W[wbits], N, mm, &T );
2011 
2012             state--;
2013             nbits = 0;
2014             wbits = 0;
2015         }
2016     }
2017 
2018     /*
2019      * process the remaining bits
2020      */
2021     for( i = 0; i < nbits; i++ )
2022     {
2023         mpi_montmul( X, X, N, mm, &T );
2024 
2025         wbits <<= 1;
2026 
2027         if( ( wbits & ( one << wsize ) ) != 0 )
2028             mpi_montmul( X, &W[1], N, mm, &T );
2029     }
2030 
2031     /*
2032      * X = A^E * R * R^-1 mod N = A^E mod N
2033      */
2034     mpi_montred( X, N, mm, &T );
2035 
2036     if( neg && E->n != 0 && ( E->p[0] & 1 ) != 0 )
2037     {
2038         X->s = -1;
2039         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( X, N, X ) );
2040     }
2041 
2042 cleanup:
2043 
2044     for( i = ( one << ( wsize - 1 ) ); i < ( one << wsize ); i++ )
2045         mbedtls_mpi_free( &W[i] );
2046 
2047     mbedtls_mpi_free( &W[1] ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &Apos );
2048 
2049     if( _RR == NULL || _RR->p == NULL )
2050         mbedtls_mpi_free( &RR );
2051 
2052     return( ret );
2053 }
2054 
2055 /*
2056  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
2057  */
2058 int mbedtls_mpi_gcd( mbedtls_mpi *G, const mbedtls_mpi *A, const mbedtls_mpi *B )
2059 {
2060     int ret;
2061     size_t lz, lzt;
2062     mbedtls_mpi TG, TA, TB;
2063 
2064     mbedtls_mpi_init( &TG ); mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
2065 
2066     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) );
2067     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
2068 
2069     lz = mbedtls_mpi_lsb( &TA );
2070     lzt = mbedtls_mpi_lsb( &TB );
2071 
2072     if( lzt < lz )
2073         lz = lzt;
2074 
2075     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, lz ) );
2076     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, lz ) );
2077 
2078     TA.s = TB.s = 1;
2079 
2080     while( mbedtls_mpi_cmp_int( &TA, 0 ) != 0 )
2081     {
2082         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, mbedtls_mpi_lsb( &TA ) ) );
2083         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, mbedtls_mpi_lsb( &TB ) ) );
2084 
2085         if( mbedtls_mpi_cmp_mpi( &TA, &TB ) >= 0 )
2086         {
2087             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TA, &TA, &TB ) );
2088             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, 1 ) );
2089         }
2090         else
2091         {
2092             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TB, &TB, &TA ) );
2093             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, 1 ) );
2094         }
2095     }
2096 
2097     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &TB, lz ) );
2098     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( G, &TB ) );
2099 
2100 cleanup:
2101 
2102     mbedtls_mpi_free( &TG ); mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TB );
2103 
2104     return( ret );
2105 }
2106 
2107 /*
2108  * Fill X with size bytes of random.
2109  *
2110  * Use a temporary bytes representation to make sure the result is the same
2111  * regardless of the platform endianness (useful when f_rng is actually
2112  * deterministic, eg for tests).
2113  */
2114 int mbedtls_mpi_fill_random( mbedtls_mpi *X, size_t size,
2115                      int (*f_rng)(void *, unsigned char *, size_t),
2116                      void *p_rng )
2117 {
2118     int ret;
2119     unsigned char buf[MBEDTLS_MPI_MAX_SIZE];
2120 
2121     if( size > MBEDTLS_MPI_MAX_SIZE )
2122         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
2123 
2124     MBEDTLS_MPI_CHK( f_rng( p_rng, buf, size ) );
2125     MBEDTLS_MPI_CHK( mbedtls_mpi_read_binary( X, buf, size ) );
2126 
2127 cleanup:
2128     mbedtls_zeroize( buf, sizeof( buf ) );
2129     return( ret );
2130 }
2131 
2132 /*
2133  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
2134  */
2135 int mbedtls_mpi_inv_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *N )
2136 {
2137     int ret;
2138     mbedtls_mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
2139 
2140     if( mbedtls_mpi_cmp_int( N, 1 ) <= 0 )
2141         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
2142 
2143     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TU ); mbedtls_mpi_init( &U1 ); mbedtls_mpi_init( &U2 );
2144     mbedtls_mpi_init( &G ); mbedtls_mpi_init( &TB ); mbedtls_mpi_init( &TV );
2145     mbedtls_mpi_init( &V1 ); mbedtls_mpi_init( &V2 );
2146 
2147     MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &G, A, N ) );
2148 
2149     if( mbedtls_mpi_cmp_int( &G, 1 ) != 0 )
2150     {
2151         ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2152         goto cleanup;
2153     }
2154 
2155     MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &TA, A, N ) );
2156     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TU, &TA ) );
2157     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, N ) );
2158     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TV, N ) );
2159 
2160     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U1, 1 ) );
2161     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U2, 0 ) );
2162     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V1, 0 ) );
2163     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V2, 1 ) );
2164 
2165     do
2166     {
2167         while( ( TU.p[0] & 1 ) == 0 )
2168         {
2169             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TU, 1 ) );
2170 
2171             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
2172             {
2173                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &U1, &U1, &TB ) );
2174                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &TA ) );
2175             }
2176 
2177             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U1, 1 ) );
2178             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U2, 1 ) );
2179         }
2180 
2181         while( ( TV.p[0] & 1 ) == 0 )
2182         {
2183             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TV, 1 ) );
2184 
2185             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
2186             {
2187                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, &TB ) );
2188                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &TA ) );
2189             }
2190 
2191             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V1, 1 ) );
2192             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V2, 1 ) );
2193         }
2194 
2195         if( mbedtls_mpi_cmp_mpi( &TU, &TV ) >= 0 )
2196         {
2197             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TU, &TU, &TV ) );
2198             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U1, &U1, &V1 ) );
2199             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &V2 ) );
2200         }
2201         else
2202         {
2203             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TV, &TV, &TU ) );
2204             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, &U1 ) );
2205             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &U2 ) );
2206         }
2207     }
2208     while( mbedtls_mpi_cmp_int( &TU, 0 ) != 0 );
2209 
2210     while( mbedtls_mpi_cmp_int( &V1, 0 ) < 0 )
2211         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, N ) );
2212 
2213     while( mbedtls_mpi_cmp_mpi( &V1, N ) >= 0 )
2214         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, N ) );
2215 
2216     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &V1 ) );
2217 
2218 cleanup:
2219 
2220     mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TU ); mbedtls_mpi_free( &U1 ); mbedtls_mpi_free( &U2 );
2221     mbedtls_mpi_free( &G ); mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TV );
2222     mbedtls_mpi_free( &V1 ); mbedtls_mpi_free( &V2 );
2223 
2224     return( ret );
2225 }
2226 
2227 #if defined(MBEDTLS_GENPRIME)
2228 
2229 static const int small_prime[] =
2230 {
2231         3,    5,    7,   11,   13,   17,   19,   23,
2232        29,   31,   37,   41,   43,   47,   53,   59,
2233        61,   67,   71,   73,   79,   83,   89,   97,
2234       101,  103,  107,  109,  113,  127,  131,  137,
2235       139,  149,  151,  157,  163,  167,  173,  179,
2236       181,  191,  193,  197,  199,  211,  223,  227,
2237       229,  233,  239,  241,  251,  257,  263,  269,
2238       271,  277,  281,  283,  293,  307,  311,  313,
2239       317,  331,  337,  347,  349,  353,  359,  367,
2240       373,  379,  383,  389,  397,  401,  409,  419,
2241       421,  431,  433,  439,  443,  449,  457,  461,
2242       463,  467,  479,  487,  491,  499,  503,  509,
2243       521,  523,  541,  547,  557,  563,  569,  571,
2244       577,  587,  593,  599,  601,  607,  613,  617,
2245       619,  631,  641,  643,  647,  653,  659,  661,
2246       673,  677,  683,  691,  701,  709,  719,  727,
2247       733,  739,  743,  751,  757,  761,  769,  773,
2248       787,  797,  809,  811,  821,  823,  827,  829,
2249       839,  853,  857,  859,  863,  877,  881,  883,
2250       887,  907,  911,  919,  929,  937,  941,  947,
2251       953,  967,  971,  977,  983,  991,  997, -103
2252 };
2253 
2254 /*
2255  * Small divisors test (X must be positive)
2256  *
2257  * Return values:
2258  * 0: no small factor (possible prime, more tests needed)
2259  * 1: certain prime
2260  * MBEDTLS_ERR_MPI_NOT_ACCEPTABLE: certain non-prime
2261  * other negative: error
2262  */
2263 static int mpi_check_small_factors( const mbedtls_mpi *X )
2264 {
2265     int ret = 0;
2266     size_t i;
2267     mbedtls_mpi_uint r;
2268 
2269     if( ( X->p[0] & 1 ) == 0 )
2270         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2271 
2272     for( i = 0; small_prime[i] > 0; i++ )
2273     {
2274         if( mbedtls_mpi_cmp_int( X, small_prime[i] ) <= 0 )
2275             return( 1 );
2276 
2277         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, small_prime[i] ) );
2278 
2279         if( r == 0 )
2280             return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2281     }
2282 
2283 cleanup:
2284     return( ret );
2285 }
2286 
2287 /*
2288  * Miller-Rabin pseudo-primality test  (HAC 4.24)
2289  */
2290 static int mpi_miller_rabin( const mbedtls_mpi *X, size_t rounds,
2291                              int (*f_rng)(void *, unsigned char *, size_t),
2292                              void *p_rng )
2293 {
2294     int ret, count;
2295     size_t i, j, k, s;
2296     mbedtls_mpi W, R, T, A, RR;
2297 
2298     mbedtls_mpi_init( &W ); mbedtls_mpi_init( &R ); mbedtls_mpi_init( &T ); mbedtls_mpi_init( &A );
2299     mbedtls_mpi_init( &RR );
2300 
2301     /*
2302      * W = |X| - 1
2303      * R = W >> lsb( W )
2304      */
2305     MBEDTLS_MPI_CHK( mbedtls_mpi_sub_int( &W, X, 1 ) );
2306     s = mbedtls_mpi_lsb( &W );
2307     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &R, &W ) );
2308     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &R, s ) );
2309 
2310     i = mbedtls_mpi_bitlen( X );
2311 
2312     for( i = 0; i < rounds; i++ )
2313     {
2314         /*
2315          * pick a random A, 1 < A < |X| - 1
2316          */
2317         count = 0;
2318         do {
2319             MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
2320 
2321             j = mbedtls_mpi_bitlen( &A );
2322             k = mbedtls_mpi_bitlen( &W );
2323             if (j > k) {
2324                 A.p[A.n - 1] &= ( (mbedtls_mpi_uint) 1 << ( k - ( A.n - 1 ) * biL - 1 ) ) - 1;
2325             }
2326 
2327             if (count++ > 30) {
2328                 ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2329                 goto cleanup;
2330             }
2331 
2332         } while ( mbedtls_mpi_cmp_mpi( &A, &W ) >= 0 ||
2333                   mbedtls_mpi_cmp_int( &A, 1 )  <= 0    );
2334 
2335         /*
2336          * A = A^R mod |X|
2337          */
2338         MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &A, &A, &R, X, &RR ) );
2339 
2340         if( mbedtls_mpi_cmp_mpi( &A, &W ) == 0 ||
2341             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2342             continue;
2343 
2344         j = 1;
2345         while( j < s && mbedtls_mpi_cmp_mpi( &A, &W ) != 0 )
2346         {
2347             /*
2348              * A = A * A mod |X|
2349              */
2350             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &T, &A, &A ) );
2351             MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &A, &T, X  ) );
2352 
2353             if( mbedtls_mpi_cmp_int( &A, 1 ) == 0 )
2354                 break;
2355 
2356             j++;
2357         }
2358 
2359         /*
2360          * not prime if A != |X| - 1 or A == 1
2361          */
2362         if( mbedtls_mpi_cmp_mpi( &A, &W ) != 0 ||
2363             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2364         {
2365             ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2366             break;
2367         }
2368     }
2369 
2370 cleanup:
2371     mbedtls_mpi_free( &W ); mbedtls_mpi_free( &R ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &A );
2372     mbedtls_mpi_free( &RR );
2373 
2374     return( ret );
2375 }
2376 
2377 /*
2378  * Pseudo-primality test: small factors, then Miller-Rabin
2379  */
2380 static int mpi_is_prime_internal( const mbedtls_mpi *X, int rounds,
2381                   int (*f_rng)(void *, unsigned char *, size_t),
2382                   void *p_rng )
2383 {
2384     int ret;
2385     mbedtls_mpi XX;
2386 
2387     XX.s = 1;
2388     XX.n = X->n;
2389     XX.p = X->p;
2390 
2391     if( mbedtls_mpi_cmp_int( &XX, 0 ) == 0 ||
2392         mbedtls_mpi_cmp_int( &XX, 1 ) == 0 )
2393         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2394 
2395     if( mbedtls_mpi_cmp_int( &XX, 2 ) == 0 )
2396         return( 0 );
2397 
2398     if( ( ret = mpi_check_small_factors( &XX ) ) != 0 )
2399     {
2400         if( ret == 1 )
2401             return( 0 );
2402 
2403         return( ret );
2404     }
2405 
2406     return( mpi_miller_rabin( &XX, rounds, f_rng, p_rng ) );
2407 }
2408 
2409 /*
2410  * Pseudo-primality test, error probability 2^-80
2411  */
2412 int mbedtls_mpi_is_prime( const mbedtls_mpi *X,
2413                   int (*f_rng)(void *, unsigned char *, size_t),
2414                   void *p_rng )
2415 {
2416     return mpi_is_prime_internal( X, 40, f_rng, p_rng );
2417 }
2418 
2419 /*
2420  * Prime number generation
2421  */
2422 int mbedtls_mpi_gen_prime( mbedtls_mpi *X, size_t nbits, int dh_flag,
2423                    int (*f_rng)(void *, unsigned char *, size_t),
2424                    void *p_rng )
2425 {
2426     int ret;
2427     size_t k, n;
2428     int rounds;
2429     mbedtls_mpi_uint r;
2430     mbedtls_mpi Y;
2431 
2432     if( nbits < 3 || nbits > MBEDTLS_MPI_MAX_BITS )
2433         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
2434 
2435     mbedtls_mpi_init( &Y );
2436 
2437     n = BITS_TO_LIMBS( nbits );
2438 
2439     /*
2440      * 2^-80 error probability, number of rounds chosen per HAC, table 4.4
2441      */
2442     rounds = ( ( nbits >= 1300 ) ?  2 : ( nbits >=  850 ) ?  3 :
2443                ( nbits >=  650 ) ?  4 : ( nbits >=  350 ) ?  8 :
2444                ( nbits >=  250 ) ? 12 : ( nbits >=  150 ) ? 18 : 27 );
2445 
2446     MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
2447 
2448     k = mbedtls_mpi_bitlen( X );
2449     if( k > nbits ) MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( X, k - nbits + 1 ) );
2450 
2451     mbedtls_mpi_set_bit( X, nbits-1, 1 );
2452 
2453     X->p[0] |= 1;
2454 
2455     if( dh_flag == 0 )
2456     {
2457         while( ( ret = mpi_is_prime_internal( X, rounds, f_rng, p_rng ) ) != 0 )
2458         {
2459             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2460                 goto cleanup;
2461 
2462             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 2 ) );
2463         }
2464     }
2465     else
2466     {
2467         /*
2468          * An necessary condition for Y and X = 2Y + 1 to be prime
2469          * is X = 2 mod 3 (which is equivalent to Y = 2 mod 3).
2470          * Make sure it is satisfied, while keeping X = 3 mod 4
2471          */
2472 
2473         X->p[0] |= 2;
2474 
2475         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, 3 ) );
2476         if( r == 0 )
2477             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 8 ) );
2478         else if( r == 1 )
2479             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 4 ) );
2480 
2481         /* Set Y = (X-1) / 2, which is X / 2 because X is odd */
2482         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, X ) );
2483         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, 1 ) );
2484 
2485         while( 1 )
2486         {
2487             /*
2488              * First, check small factors for X and Y
2489              * before doing Miller-Rabin on any of them
2490              */
2491             if( ( ret = mpi_check_small_factors(  X         ) ) == 0 &&
2492                 ( ret = mpi_check_small_factors( &Y         ) ) == 0 &&
2493                 ( ret = mpi_miller_rabin(  X, rounds, f_rng, p_rng  ) )
2494                                                                 == 0 &&
2495                 ( ret = mpi_miller_rabin( &Y, rounds, f_rng, p_rng  ) )
2496                                                                 == 0 )
2497             {
2498                 break;
2499             }
2500 
2501             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2502                 goto cleanup;
2503 
2504             /*
2505              * Next candidates. We want to preserve Y = (X-1) / 2 and
2506              * Y = 1 mod 2 and Y = 2 mod 3 (eq X = 3 mod 4 and X = 2 mod 3)
2507              * so up Y by 6 and X by 12.
2508              */
2509             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int(  X,  X, 12 ) );
2510             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( &Y, &Y, 6  ) );
2511         }
2512     }
2513 
2514 cleanup:
2515 
2516     mbedtls_mpi_free( &Y );
2517 
2518     return( ret );
2519 }
2520 
2521 #endif /* MBEDTLS_GENPRIME */
2522 
2523 #if defined(MBEDTLS_SELF_TEST)
2524 
2525 #define GCD_PAIR_COUNT  3
2526 
2527 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
2528 {
2529     { 693, 609, 21 },
2530     { 1764, 868, 28 },
2531     { 768454923, 542167814, 1 }
2532 };
2533 
2534 /*
2535  * Checkup routine
2536  */
2537 int mbedtls_mpi_self_test( int verbose )
2538 {
2539     int ret, i;
2540     mbedtls_mpi A, E, N, X, Y, U, V;
2541 
2542     mbedtls_mpi_init( &A ); mbedtls_mpi_init( &E ); mbedtls_mpi_init( &N ); mbedtls_mpi_init( &X );
2543     mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &U ); mbedtls_mpi_init( &V );
2544 
2545     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &A, 16,
2546         "EFE021C2645FD1DC586E69184AF4A31E" \
2547         "D5F53E93B5F123FA41680867BA110131" \
2548         "944FE7952E2517337780CB0DB80E61AA" \
2549         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
2550 
2551     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &E, 16,
2552         "B2E7EFD37075B9F03FF989C7C5051C20" \
2553         "34D2A323810251127E7BF8625A4F49A5" \
2554         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2555         "5B5C25763222FEFCCFC38B832366C29E" ) );
2556 
2557     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &N, 16,
2558         "0066A198186C18C10B2F5ED9B522752A" \
2559         "9830B69916E535C8F047518A889A43A5" \
2560         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2561 
2562     MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &X, &A, &N ) );
2563 
2564     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2565         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2566         "9E857EA95A03512E2BAE7391688D264A" \
2567         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2568         "8001B72E848A38CAE1C65F78E56ABDEF" \
2569         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2570         "ECF677152EF804370C1A305CAF3B5BF1" \
2571         "30879B56C61DE584A0F53A2447A51E" ) );
2572 
2573     if( verbose != 0 )
2574         mbedtls_printf( "  MPI test #1 (mul_mpi): " );
2575 
2576     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2577     {
2578         if( verbose != 0 )
2579             mbedtls_printf( "failed\n" );
2580 
2581         ret = 1;
2582         goto cleanup;
2583     }
2584 
2585     if( verbose != 0 )
2586         mbedtls_printf( "passed\n" );
2587 
2588     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( &X, &Y, &A, &N ) );
2589 
2590     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2591         "256567336059E52CAE22925474705F39A94" ) );
2592 
2593     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &V, 16,
2594         "6613F26162223DF488E9CD48CC132C7A" \
2595         "0AC93C701B001B092E4E5B9F73BCD27B" \
2596         "9EE50D0657C77F374E903CDFA4C642" ) );
2597 
2598     if( verbose != 0 )
2599         mbedtls_printf( "  MPI test #2 (div_mpi): " );
2600 
2601     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 ||
2602         mbedtls_mpi_cmp_mpi( &Y, &V ) != 0 )
2603     {
2604         if( verbose != 0 )
2605             mbedtls_printf( "failed\n" );
2606 
2607         ret = 1;
2608         goto cleanup;
2609     }
2610 
2611     if( verbose != 0 )
2612         mbedtls_printf( "passed\n" );
2613 
2614     MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2615 
2616     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2617         "36E139AEA55215609D2816998ED020BB" \
2618         "BD96C37890F65171D948E9BC7CBAA4D9" \
2619         "325D24D6A3C12710F10A09FA08AB87" ) );
2620 
2621     if( verbose != 0 )
2622         mbedtls_printf( "  MPI test #3 (exp_mod): " );
2623 
2624     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2625     {
2626         if( verbose != 0 )
2627             mbedtls_printf( "failed\n" );
2628 
2629         ret = 1;
2630         goto cleanup;
2631     }
2632 
2633     if( verbose != 0 )
2634         mbedtls_printf( "passed\n" );
2635 
2636     MBEDTLS_MPI_CHK( mbedtls_mpi_inv_mod( &X, &A, &N ) );
2637 
2638     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2639         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2640         "C3DBA76456363A10869622EAC2DD84EC" \
2641         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2642 
2643     if( verbose != 0 )
2644         mbedtls_printf( "  MPI test #4 (inv_mod): " );
2645 
2646     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2647     {
2648         if( verbose != 0 )
2649             mbedtls_printf( "failed\n" );
2650 
2651         ret = 1;
2652         goto cleanup;
2653     }
2654 
2655     if( verbose != 0 )
2656         mbedtls_printf( "passed\n" );
2657 
2658     if( verbose != 0 )
2659         mbedtls_printf( "  MPI test #5 (simple gcd): " );
2660 
2661     for( i = 0; i < GCD_PAIR_COUNT; i++ )
2662     {
2663         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &X, gcd_pairs[i][0] ) );
2664         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Y, gcd_pairs[i][1] ) );
2665 
2666         MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &A, &X, &Y ) );
2667 
2668         if( mbedtls_mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2669         {
2670             if( verbose != 0 )
2671                 mbedtls_printf( "failed at %d\n", i );
2672 
2673             ret = 1;
2674             goto cleanup;
2675         }
2676     }
2677 
2678     if( verbose != 0 )
2679         mbedtls_printf( "passed\n" );
2680 
2681 cleanup:
2682 
2683     if( ret != 0 && verbose != 0 )
2684         mbedtls_printf( "Unexpected error, return code = %08X\n", ret );
2685 
2686     mbedtls_mpi_free( &A ); mbedtls_mpi_free( &E ); mbedtls_mpi_free( &N ); mbedtls_mpi_free( &X );
2687     mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &U ); mbedtls_mpi_free( &V );
2688 
2689     if( verbose != 0 )
2690         mbedtls_printf( "\n" );
2691 
2692     return( ret );
2693 }
2694 
2695 #endif /* MBEDTLS_SELF_TEST */
2696 
2697 #endif /* MBEDTLS_BIGNUM_C */
2698