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