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