xref: /reactos/dll/3rdparty/mbedtls/bignum.c (revision 279107d5)
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, char **p )
508 {
509     int ret;
510     mbedtls_mpi_uint r;
511 
512     if( radix < 2 || radix > 16 )
513         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
514 
515     MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, radix ) );
516     MBEDTLS_MPI_CHK( mbedtls_mpi_div_int( X, NULL, X, radix ) );
517 
518     if( mbedtls_mpi_cmp_int( X, 0 ) != 0 )
519         MBEDTLS_MPI_CHK( mpi_write_hlp( X, radix, p ) );
520 
521     if( r < 10 )
522         *(*p)++ = (char)( r + 0x30 );
523     else
524         *(*p)++ = (char)( r + 0x37 );
525 
526 cleanup:
527 
528     return( ret );
529 }
530 
531 /*
532  * Export into an ASCII string
533  */
534 int mbedtls_mpi_write_string( const mbedtls_mpi *X, int radix,
535                               char *buf, size_t buflen, size_t *olen )
536 {
537     int ret = 0;
538     size_t n;
539     char *p;
540     mbedtls_mpi T;
541 
542     if( radix < 2 || radix > 16 )
543         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
544 
545     n = mbedtls_mpi_bitlen( X );
546     if( radix >=  4 ) n >>= 1;
547     if( radix >= 16 ) n >>= 1;
548     /*
549      * Round up the buffer length to an even value to ensure that there is
550      * enough room for hexadecimal values that can be represented in an odd
551      * number of digits.
552      */
553     n += 3 + ( ( n + 1 ) & 1 );
554 
555     if( buflen < n )
556     {
557         *olen = n;
558         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
559     }
560 
561     p = buf;
562     mbedtls_mpi_init( &T );
563 
564     if( X->s == -1 )
565         *p++ = '-';
566 
567     if( radix == 16 )
568     {
569         int c;
570         size_t i, j, k;
571 
572         for( i = X->n, k = 0; i > 0; i-- )
573         {
574             for( j = ciL; j > 0; j-- )
575             {
576                 c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
577 
578                 if( c == 0 && k == 0 && ( i + j ) != 2 )
579                     continue;
580 
581                 *(p++) = "0123456789ABCDEF" [c / 16];
582                 *(p++) = "0123456789ABCDEF" [c % 16];
583                 k = 1;
584             }
585         }
586     }
587     else
588     {
589         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T, X ) );
590 
591         if( T.s == -1 )
592             T.s = 1;
593 
594         MBEDTLS_MPI_CHK( mpi_write_hlp( &T, radix, &p ) );
595     }
596 
597     *p++ = '\0';
598     *olen = p - buf;
599 
600 cleanup:
601 
602     mbedtls_mpi_free( &T );
603 
604     return( ret );
605 }
606 
607 #if defined(MBEDTLS_FS_IO)
608 /*
609  * Read X from an opened file
610  */
611 int mbedtls_mpi_read_file( mbedtls_mpi *X, int radix, FILE *fin )
612 {
613     mbedtls_mpi_uint d;
614     size_t slen;
615     char *p;
616     /*
617      * Buffer should have space for (short) label and decimal formatted MPI,
618      * newline characters and '\0'
619      */
620     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
621 
622     memset( s, 0, sizeof( s ) );
623     if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
624         return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
625 
626     slen = strlen( s );
627     if( slen == sizeof( s ) - 2 )
628         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
629 
630     if( slen > 0 && s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
631     if( slen > 0 && s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }
632 
633     p = s + slen;
634     while( p-- > s )
635         if( mpi_get_digit( &d, radix, *p ) != 0 )
636             break;
637 
638     return( mbedtls_mpi_read_string( X, radix, p + 1 ) );
639 }
640 
641 /*
642  * Write X into an opened file (or stdout if fout == NULL)
643  */
644 int mbedtls_mpi_write_file( const char *p, const mbedtls_mpi *X, int radix, FILE *fout )
645 {
646     int ret;
647     size_t n, slen, plen;
648     /*
649      * Buffer should have space for (short) label and decimal formatted MPI,
650      * newline characters and '\0'
651      */
652     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
653 
654     memset( s, 0, sizeof( s ) );
655 
656     MBEDTLS_MPI_CHK( mbedtls_mpi_write_string( X, radix, s, sizeof( s ) - 2, &n ) );
657 
658     if( p == NULL ) p = "";
659 
660     plen = strlen( p );
661     slen = strlen( s );
662     s[slen++] = '\r';
663     s[slen++] = '\n';
664 
665     if( fout != NULL )
666     {
667         if( fwrite( p, 1, plen, fout ) != plen ||
668             fwrite( s, 1, slen, fout ) != slen )
669             return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
670     }
671     else
672         mbedtls_printf( "%s%s", p, s );
673 
674 cleanup:
675 
676     return( ret );
677 }
678 #endif /* MBEDTLS_FS_IO */
679 
680 /*
681  * Import X from unsigned binary data, big endian
682  */
683 int mbedtls_mpi_read_binary( mbedtls_mpi *X, const unsigned char *buf, size_t buflen )
684 {
685     int ret;
686     size_t i, j;
687     size_t const limbs = CHARS_TO_LIMBS( buflen );
688 
689     /* Ensure that target MPI has exactly the necessary number of limbs */
690     if( X->n != limbs )
691     {
692         mbedtls_mpi_free( X );
693         mbedtls_mpi_init( X );
694         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, limbs ) );
695     }
696 
697     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
698 
699     for( i = buflen, j = 0; i > 0; i--, j++ )
700         X->p[j / ciL] |= ((mbedtls_mpi_uint) buf[i - 1]) << ((j % ciL) << 3);
701 
702 cleanup:
703 
704     return( ret );
705 }
706 
707 /*
708  * Export X into unsigned binary data, big endian
709  */
710 int mbedtls_mpi_write_binary( const mbedtls_mpi *X,
711                               unsigned char *buf, size_t buflen )
712 {
713     size_t stored_bytes = X->n * ciL;
714     size_t bytes_to_copy;
715     unsigned char *p;
716     size_t i;
717 
718     if( stored_bytes < buflen )
719     {
720         /* There is enough space in the output buffer. Write initial
721          * null bytes and record the position at which to start
722          * writing the significant bytes. In this case, the execution
723          * trace of this function does not depend on the value of the
724          * number. */
725         bytes_to_copy = stored_bytes;
726         p = buf + buflen - stored_bytes;
727         memset( buf, 0, buflen - stored_bytes );
728     }
729     else
730     {
731         /* The output buffer is smaller than the allocated size of X.
732          * However X may fit if its leading bytes are zero. */
733         bytes_to_copy = buflen;
734         p = buf;
735         for( i = bytes_to_copy; i < stored_bytes; i++ )
736         {
737             if( GET_BYTE( X, i ) != 0 )
738                 return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
739         }
740     }
741 
742     for( i = 0; i < bytes_to_copy; i++ )
743         p[bytes_to_copy - i - 1] = GET_BYTE( X, i );
744 
745     return( 0 );
746 }
747 
748 /*
749  * Left-shift: X <<= count
750  */
751 int mbedtls_mpi_shift_l( mbedtls_mpi *X, size_t count )
752 {
753     int ret;
754     size_t i, v0, t1;
755     mbedtls_mpi_uint r0 = 0, r1;
756 
757     v0 = count / (biL    );
758     t1 = count & (biL - 1);
759 
760     i = mbedtls_mpi_bitlen( X ) + count;
761 
762     if( X->n * biL < i )
763         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, BITS_TO_LIMBS( i ) ) );
764 
765     ret = 0;
766 
767     /*
768      * shift by count / limb_size
769      */
770     if( v0 > 0 )
771     {
772         for( i = X->n; i > v0; i-- )
773             X->p[i - 1] = X->p[i - v0 - 1];
774 
775         for( ; i > 0; i-- )
776             X->p[i - 1] = 0;
777     }
778 
779     /*
780      * shift by count % limb_size
781      */
782     if( t1 > 0 )
783     {
784         for( i = v0; i < X->n; i++ )
785         {
786             r1 = X->p[i] >> (biL - t1);
787             X->p[i] <<= t1;
788             X->p[i] |= r0;
789             r0 = r1;
790         }
791     }
792 
793 cleanup:
794 
795     return( ret );
796 }
797 
798 /*
799  * Right-shift: X >>= count
800  */
801 int mbedtls_mpi_shift_r( mbedtls_mpi *X, size_t count )
802 {
803     size_t i, v0, v1;
804     mbedtls_mpi_uint r0 = 0, r1;
805 
806     v0 = count /  biL;
807     v1 = count & (biL - 1);
808 
809     if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
810         return mbedtls_mpi_lset( X, 0 );
811 
812     /*
813      * shift by count / limb_size
814      */
815     if( v0 > 0 )
816     {
817         for( i = 0; i < X->n - v0; i++ )
818             X->p[i] = X->p[i + v0];
819 
820         for( ; i < X->n; i++ )
821             X->p[i] = 0;
822     }
823 
824     /*
825      * shift by count % limb_size
826      */
827     if( v1 > 0 )
828     {
829         for( i = X->n; i > 0; i-- )
830         {
831             r1 = X->p[i - 1] << (biL - v1);
832             X->p[i - 1] >>= v1;
833             X->p[i - 1] |= r0;
834             r0 = r1;
835         }
836     }
837 
838     return( 0 );
839 }
840 
841 /*
842  * Compare unsigned values
843  */
844 int mbedtls_mpi_cmp_abs( const mbedtls_mpi *X, const mbedtls_mpi *Y )
845 {
846     size_t i, j;
847 
848     for( i = X->n; i > 0; i-- )
849         if( X->p[i - 1] != 0 )
850             break;
851 
852     for( j = Y->n; j > 0; j-- )
853         if( Y->p[j - 1] != 0 )
854             break;
855 
856     if( i == 0 && j == 0 )
857         return( 0 );
858 
859     if( i > j ) return(  1 );
860     if( j > i ) return( -1 );
861 
862     for( ; i > 0; i-- )
863     {
864         if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
865         if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
866     }
867 
868     return( 0 );
869 }
870 
871 /*
872  * Compare signed values
873  */
874 int mbedtls_mpi_cmp_mpi( const mbedtls_mpi *X, const mbedtls_mpi *Y )
875 {
876     size_t i, j;
877 
878     for( i = X->n; i > 0; i-- )
879         if( X->p[i - 1] != 0 )
880             break;
881 
882     for( j = Y->n; j > 0; j-- )
883         if( Y->p[j - 1] != 0 )
884             break;
885 
886     if( i == 0 && j == 0 )
887         return( 0 );
888 
889     if( i > j ) return(  X->s );
890     if( j > i ) return( -Y->s );
891 
892     if( X->s > 0 && Y->s < 0 ) return(  1 );
893     if( Y->s > 0 && X->s < 0 ) return( -1 );
894 
895     for( ; i > 0; i-- )
896     {
897         if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
898         if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
899     }
900 
901     return( 0 );
902 }
903 
904 /*
905  * Compare signed values
906  */
907 int mbedtls_mpi_cmp_int( const mbedtls_mpi *X, mbedtls_mpi_sint z )
908 {
909     mbedtls_mpi Y;
910     mbedtls_mpi_uint p[1];
911 
912     *p  = ( z < 0 ) ? -z : z;
913     Y.s = ( z < 0 ) ? -1 : 1;
914     Y.n = 1;
915     Y.p = p;
916 
917     return( mbedtls_mpi_cmp_mpi( X, &Y ) );
918 }
919 
920 /*
921  * Unsigned addition: X = |A| + |B|  (HAC 14.7)
922  */
923 int mbedtls_mpi_add_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
924 {
925     int ret;
926     size_t i, j;
927     mbedtls_mpi_uint *o, *p, c, tmp;
928 
929     if( X == B )
930     {
931         const mbedtls_mpi *T = A; A = X; B = T;
932     }
933 
934     if( X != A )
935         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
936 
937     /*
938      * X should always be positive as a result of unsigned additions.
939      */
940     X->s = 1;
941 
942     for( j = B->n; j > 0; j-- )
943         if( B->p[j - 1] != 0 )
944             break;
945 
946     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
947 
948     o = B->p; p = X->p; c = 0;
949 
950     /*
951      * tmp is used because it might happen that p == o
952      */
953     for( i = 0; i < j; i++, o++, p++ )
954     {
955         tmp= *o;
956         *p +=  c; c  = ( *p <  c );
957         *p += tmp; c += ( *p < tmp );
958     }
959 
960     while( c != 0 )
961     {
962         if( i >= X->n )
963         {
964             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + 1 ) );
965             p = X->p + i;
966         }
967 
968         *p += c; c = ( *p < c ); i++; p++;
969     }
970 
971 cleanup:
972 
973     return( ret );
974 }
975 
976 /*
977  * Helper for mbedtls_mpi subtraction
978  */
979 static void mpi_sub_hlp( size_t n, mbedtls_mpi_uint *s, mbedtls_mpi_uint *d )
980 {
981     size_t i;
982     mbedtls_mpi_uint c, z;
983 
984     for( i = c = 0; i < n; i++, s++, d++ )
985     {
986         z = ( *d <  c );     *d -=  c;
987         c = ( *d < *s ) + z; *d -= *s;
988     }
989 
990     while( c != 0 )
991     {
992         z = ( *d < c ); *d -= c;
993         c = z; i++; d++;
994     }
995 }
996 
997 /*
998  * Unsigned subtraction: X = |A| - |B|  (HAC 14.9)
999  */
1000 int mbedtls_mpi_sub_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1001 {
1002     mbedtls_mpi TB;
1003     int ret;
1004     size_t n;
1005 
1006     if( mbedtls_mpi_cmp_abs( A, B ) < 0 )
1007         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1008 
1009     mbedtls_mpi_init( &TB );
1010 
1011     if( X == B )
1012     {
1013         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
1014         B = &TB;
1015     }
1016 
1017     if( X != A )
1018         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
1019 
1020     /*
1021      * X should always be positive as a result of unsigned subtractions.
1022      */
1023     X->s = 1;
1024 
1025     ret = 0;
1026 
1027     for( n = B->n; n > 0; n-- )
1028         if( B->p[n - 1] != 0 )
1029             break;
1030 
1031     mpi_sub_hlp( n, B->p, X->p );
1032 
1033 cleanup:
1034 
1035     mbedtls_mpi_free( &TB );
1036 
1037     return( ret );
1038 }
1039 
1040 /*
1041  * Signed addition: X = A + B
1042  */
1043 int mbedtls_mpi_add_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1044 {
1045     int ret, s = A->s;
1046 
1047     if( A->s * B->s < 0 )
1048     {
1049         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1050         {
1051             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1052             X->s =  s;
1053         }
1054         else
1055         {
1056             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1057             X->s = -s;
1058         }
1059     }
1060     else
1061     {
1062         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1063         X->s = s;
1064     }
1065 
1066 cleanup:
1067 
1068     return( ret );
1069 }
1070 
1071 /*
1072  * Signed subtraction: X = A - B
1073  */
1074 int mbedtls_mpi_sub_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1075 {
1076     int ret, s = A->s;
1077 
1078     if( A->s * B->s > 0 )
1079     {
1080         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1081         {
1082             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1083             X->s =  s;
1084         }
1085         else
1086         {
1087             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1088             X->s = -s;
1089         }
1090     }
1091     else
1092     {
1093         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1094         X->s = s;
1095     }
1096 
1097 cleanup:
1098 
1099     return( ret );
1100 }
1101 
1102 /*
1103  * Signed addition: X = A + b
1104  */
1105 int mbedtls_mpi_add_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1106 {
1107     mbedtls_mpi _B;
1108     mbedtls_mpi_uint p[1];
1109 
1110     p[0] = ( b < 0 ) ? -b : b;
1111     _B.s = ( b < 0 ) ? -1 : 1;
1112     _B.n = 1;
1113     _B.p = p;
1114 
1115     return( mbedtls_mpi_add_mpi( X, A, &_B ) );
1116 }
1117 
1118 /*
1119  * Signed subtraction: X = A - b
1120  */
1121 int mbedtls_mpi_sub_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1122 {
1123     mbedtls_mpi _B;
1124     mbedtls_mpi_uint p[1];
1125 
1126     p[0] = ( b < 0 ) ? -b : b;
1127     _B.s = ( b < 0 ) ? -1 : 1;
1128     _B.n = 1;
1129     _B.p = p;
1130 
1131     return( mbedtls_mpi_sub_mpi( X, A, &_B ) );
1132 }
1133 
1134 /*
1135  * Helper for mbedtls_mpi multiplication
1136  */
1137 static
1138 #if defined(__APPLE__) && defined(__arm__)
1139 /*
1140  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
1141  * appears to need this to prevent bad ARM code generation at -O3.
1142  */
1143 __attribute__ ((noinline))
1144 #endif
1145 void mpi_mul_hlp( size_t i, mbedtls_mpi_uint *s, mbedtls_mpi_uint *d, mbedtls_mpi_uint b )
1146 {
1147     mbedtls_mpi_uint c = 0, t = 0;
1148 
1149 #if defined(MULADDC_HUIT)
1150     for( ; i >= 8; i -= 8 )
1151     {
1152         MULADDC_INIT
1153         MULADDC_HUIT
1154         MULADDC_STOP
1155     }
1156 
1157     for( ; i > 0; i-- )
1158     {
1159         MULADDC_INIT
1160         MULADDC_CORE
1161         MULADDC_STOP
1162     }
1163 #else /* MULADDC_HUIT */
1164     for( ; i >= 16; i -= 16 )
1165     {
1166         MULADDC_INIT
1167         MULADDC_CORE   MULADDC_CORE
1168         MULADDC_CORE   MULADDC_CORE
1169         MULADDC_CORE   MULADDC_CORE
1170         MULADDC_CORE   MULADDC_CORE
1171 
1172         MULADDC_CORE   MULADDC_CORE
1173         MULADDC_CORE   MULADDC_CORE
1174         MULADDC_CORE   MULADDC_CORE
1175         MULADDC_CORE   MULADDC_CORE
1176         MULADDC_STOP
1177     }
1178 
1179     for( ; i >= 8; i -= 8 )
1180     {
1181         MULADDC_INIT
1182         MULADDC_CORE   MULADDC_CORE
1183         MULADDC_CORE   MULADDC_CORE
1184 
1185         MULADDC_CORE   MULADDC_CORE
1186         MULADDC_CORE   MULADDC_CORE
1187         MULADDC_STOP
1188     }
1189 
1190     for( ; i > 0; i-- )
1191     {
1192         MULADDC_INIT
1193         MULADDC_CORE
1194         MULADDC_STOP
1195     }
1196 #endif /* MULADDC_HUIT */
1197 
1198     t++;
1199 
1200     do {
1201         *d += c; c = ( *d < c ); d++;
1202     }
1203     while( c != 0 );
1204 }
1205 
1206 /*
1207  * Baseline multiplication: X = A * B  (HAC 14.12)
1208  */
1209 int mbedtls_mpi_mul_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1210 {
1211     int ret;
1212     size_t i, j;
1213     mbedtls_mpi TA, TB;
1214 
1215     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
1216 
1217     if( X == A ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) ); A = &TA; }
1218     if( X == B ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) ); B = &TB; }
1219 
1220     for( i = A->n; i > 0; i-- )
1221         if( A->p[i - 1] != 0 )
1222             break;
1223 
1224     for( j = B->n; j > 0; j-- )
1225         if( B->p[j - 1] != 0 )
1226             break;
1227 
1228     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + j ) );
1229     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
1230 
1231     for( i++; j > 0; j-- )
1232         mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1233 
1234     X->s = A->s * B->s;
1235 
1236 cleanup:
1237 
1238     mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TA );
1239 
1240     return( ret );
1241 }
1242 
1243 /*
1244  * Baseline multiplication: X = A * b
1245  */
1246 int mbedtls_mpi_mul_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_uint b )
1247 {
1248     mbedtls_mpi _B;
1249     mbedtls_mpi_uint p[1];
1250 
1251     _B.s = 1;
1252     _B.n = 1;
1253     _B.p = p;
1254     p[0] = b;
1255 
1256     return( mbedtls_mpi_mul_mpi( X, A, &_B ) );
1257 }
1258 
1259 /*
1260  * Unsigned integer divide - double mbedtls_mpi_uint dividend, u1/u0, and
1261  * mbedtls_mpi_uint divisor, d
1262  */
1263 static mbedtls_mpi_uint mbedtls_int_div_int( mbedtls_mpi_uint u1,
1264             mbedtls_mpi_uint u0, mbedtls_mpi_uint d, mbedtls_mpi_uint *r )
1265 {
1266 #if defined(MBEDTLS_HAVE_UDBL)
1267     mbedtls_t_udbl dividend, quotient;
1268 #else
1269     const mbedtls_mpi_uint radix = (mbedtls_mpi_uint) 1 << biH;
1270     const mbedtls_mpi_uint uint_halfword_mask = ( (mbedtls_mpi_uint) 1 << biH ) - 1;
1271     mbedtls_mpi_uint d0, d1, q0, q1, rAX, r0, quotient;
1272     mbedtls_mpi_uint u0_msw, u0_lsw;
1273     size_t s;
1274 #endif
1275 
1276     /*
1277      * Check for overflow
1278      */
1279     if( 0 == d || u1 >= d )
1280     {
1281         if (r != NULL) *r = ~0;
1282 
1283         return ( ~0 );
1284     }
1285 
1286 #if defined(MBEDTLS_HAVE_UDBL)
1287     dividend  = (mbedtls_t_udbl) u1 << biL;
1288     dividend |= (mbedtls_t_udbl) u0;
1289     quotient = dividend / d;
1290     if( quotient > ( (mbedtls_t_udbl) 1 << biL ) - 1 )
1291         quotient = ( (mbedtls_t_udbl) 1 << biL ) - 1;
1292 
1293     if( r != NULL )
1294         *r = (mbedtls_mpi_uint)( dividend - (quotient * d ) );
1295 
1296     return (mbedtls_mpi_uint) quotient;
1297 #else
1298 
1299     /*
1300      * Algorithm D, Section 4.3.1 - The Art of Computer Programming
1301      *   Vol. 2 - Seminumerical Algorithms, Knuth
1302      */
1303 
1304     /*
1305      * Normalize the divisor, d, and dividend, u0, u1
1306      */
1307     s = mbedtls_clz( d );
1308     d = d << s;
1309 
1310     u1 = u1 << s;
1311     u1 |= ( u0 >> ( biL - s ) ) & ( -(mbedtls_mpi_sint)s >> ( biL - 1 ) );
1312     u0 =  u0 << s;
1313 
1314     d1 = d >> biH;
1315     d0 = d & uint_halfword_mask;
1316 
1317     u0_msw = u0 >> biH;
1318     u0_lsw = u0 & uint_halfword_mask;
1319 
1320     /*
1321      * Find the first quotient and remainder
1322      */
1323     q1 = u1 / d1;
1324     r0 = u1 - d1 * q1;
1325 
1326     while( q1 >= radix || ( q1 * d0 > radix * r0 + u0_msw ) )
1327     {
1328         q1 -= 1;
1329         r0 += d1;
1330 
1331         if ( r0 >= radix ) break;
1332     }
1333 
1334     rAX = ( u1 * radix ) + ( u0_msw - q1 * d );
1335     q0 = rAX / d1;
1336     r0 = rAX - q0 * d1;
1337 
1338     while( q0 >= radix || ( q0 * d0 > radix * r0 + u0_lsw ) )
1339     {
1340         q0 -= 1;
1341         r0 += d1;
1342 
1343         if ( r0 >= radix ) break;
1344     }
1345 
1346     if (r != NULL)
1347         *r = ( rAX * radix + u0_lsw - q0 * d ) >> s;
1348 
1349     quotient = q1 * radix + q0;
1350 
1351     return quotient;
1352 #endif
1353 }
1354 
1355 /*
1356  * Division by mbedtls_mpi: A = Q * B + R  (HAC 14.20)
1357  */
1358 int mbedtls_mpi_div_mpi( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1359 {
1360     int ret;
1361     size_t i, n, t, k;
1362     mbedtls_mpi X, Y, Z, T1, T2;
1363 
1364     if( mbedtls_mpi_cmp_int( B, 0 ) == 0 )
1365         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1366 
1367     mbedtls_mpi_init( &X ); mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &Z );
1368     mbedtls_mpi_init( &T1 ); mbedtls_mpi_init( &T2 );
1369 
1370     if( mbedtls_mpi_cmp_abs( A, B ) < 0 )
1371     {
1372         if( Q != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_lset( Q, 0 ) );
1373         if( R != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, A ) );
1374         return( 0 );
1375     }
1376 
1377     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &X, A ) );
1378     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, B ) );
1379     X.s = Y.s = 1;
1380 
1381     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &Z, A->n + 2 ) );
1382     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Z,  0 ) );
1383     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T1, 2 ) );
1384     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T2, 3 ) );
1385 
1386     k = mbedtls_mpi_bitlen( &Y ) % biL;
1387     if( k < biL - 1 )
1388     {
1389         k = biL - 1 - k;
1390         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &X, k ) );
1391         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, k ) );
1392     }
1393     else k = 0;
1394 
1395     n = X.n - 1;
1396     t = Y.n - 1;
1397     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, biL * ( n - t ) ) );
1398 
1399     while( mbedtls_mpi_cmp_mpi( &X, &Y ) >= 0 )
1400     {
1401         Z.p[n - t]++;
1402         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &Y ) );
1403     }
1404     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, biL * ( n - t ) ) );
1405 
1406     for( i = n; i > t ; i-- )
1407     {
1408         if( X.p[i] >= Y.p[t] )
1409             Z.p[i - t - 1] = ~0;
1410         else
1411         {
1412             Z.p[i - t - 1] = mbedtls_int_div_int( X.p[i], X.p[i - 1],
1413                                                             Y.p[t], NULL);
1414         }
1415 
1416         Z.p[i - t - 1]++;
1417         do
1418         {
1419             Z.p[i - t - 1]--;
1420 
1421             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T1, 0 ) );
1422             T1.p[0] = ( t < 1 ) ? 0 : Y.p[t - 1];
1423             T1.p[1] = Y.p[t];
1424             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1425 
1426             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T2, 0 ) );
1427             T2.p[0] = ( i < 2 ) ? 0 : X.p[i - 2];
1428             T2.p[1] = ( i < 1 ) ? 0 : X.p[i - 1];
1429             T2.p[2] = X.p[i];
1430         }
1431         while( mbedtls_mpi_cmp_mpi( &T1, &T2 ) > 0 );
1432 
1433         MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1434         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1,  biL * ( i - t - 1 ) ) );
1435         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &T1 ) );
1436 
1437         if( mbedtls_mpi_cmp_int( &X, 0 ) < 0 )
1438         {
1439             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T1, &Y ) );
1440             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1, biL * ( i - t - 1 ) ) );
1441             MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &X, &X, &T1 ) );
1442             Z.p[i - t - 1]--;
1443         }
1444     }
1445 
1446     if( Q != NULL )
1447     {
1448         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( Q, &Z ) );
1449         Q->s = A->s * B->s;
1450     }
1451 
1452     if( R != NULL )
1453     {
1454         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &X, k ) );
1455         X.s = A->s;
1456         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, &X ) );
1457 
1458         if( mbedtls_mpi_cmp_int( R, 0 ) == 0 )
1459             R->s = 1;
1460     }
1461 
1462 cleanup:
1463 
1464     mbedtls_mpi_free( &X ); mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &Z );
1465     mbedtls_mpi_free( &T1 ); mbedtls_mpi_free( &T2 );
1466 
1467     return( ret );
1468 }
1469 
1470 /*
1471  * Division by int: A = Q * b + R
1472  */
1473 int mbedtls_mpi_div_int( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1474 {
1475     mbedtls_mpi _B;
1476     mbedtls_mpi_uint p[1];
1477 
1478     p[0] = ( b < 0 ) ? -b : b;
1479     _B.s = ( b < 0 ) ? -1 : 1;
1480     _B.n = 1;
1481     _B.p = p;
1482 
1483     return( mbedtls_mpi_div_mpi( Q, R, A, &_B ) );
1484 }
1485 
1486 /*
1487  * Modulo: R = A mod B
1488  */
1489 int mbedtls_mpi_mod_mpi( mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1490 {
1491     int ret;
1492 
1493     if( mbedtls_mpi_cmp_int( B, 0 ) < 0 )
1494         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1495 
1496     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( NULL, R, A, B ) );
1497 
1498     while( mbedtls_mpi_cmp_int( R, 0 ) < 0 )
1499       MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( R, R, B ) );
1500 
1501     while( mbedtls_mpi_cmp_mpi( R, B ) >= 0 )
1502       MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( R, R, B ) );
1503 
1504 cleanup:
1505 
1506     return( ret );
1507 }
1508 
1509 /*
1510  * Modulo: r = A mod b
1511  */
1512 int mbedtls_mpi_mod_int( mbedtls_mpi_uint *r, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1513 {
1514     size_t i;
1515     mbedtls_mpi_uint x, y, z;
1516 
1517     if( b == 0 )
1518         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1519 
1520     if( b < 0 )
1521         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1522 
1523     /*
1524      * handle trivial cases
1525      */
1526     if( b == 1 )
1527     {
1528         *r = 0;
1529         return( 0 );
1530     }
1531 
1532     if( b == 2 )
1533     {
1534         *r = A->p[0] & 1;
1535         return( 0 );
1536     }
1537 
1538     /*
1539      * general case
1540      */
1541     for( i = A->n, y = 0; i > 0; i-- )
1542     {
1543         x  = A->p[i - 1];
1544         y  = ( y << biH ) | ( x >> biH );
1545         z  = y / b;
1546         y -= z * b;
1547 
1548         x <<= biH;
1549         y  = ( y << biH ) | ( x >> biH );
1550         z  = y / b;
1551         y -= z * b;
1552     }
1553 
1554     /*
1555      * If A is negative, then the current y represents a negative value.
1556      * Flipping it to the positive side.
1557      */
1558     if( A->s < 0 && y != 0 )
1559         y = b - y;
1560 
1561     *r = y;
1562 
1563     return( 0 );
1564 }
1565 
1566 /*
1567  * Fast Montgomery initialization (thanks to Tom St Denis)
1568  */
1569 static void mpi_montg_init( mbedtls_mpi_uint *mm, const mbedtls_mpi *N )
1570 {
1571     mbedtls_mpi_uint x, m0 = N->p[0];
1572     unsigned int i;
1573 
1574     x  = m0;
1575     x += ( ( m0 + 2 ) & 4 ) << 1;
1576 
1577     for( i = biL; i >= 8; i /= 2 )
1578         x *= ( 2 - ( m0 * x ) );
1579 
1580     *mm = ~x + 1;
1581 }
1582 
1583 /*
1584  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1585  */
1586 static int mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,
1587                          const mbedtls_mpi *T )
1588 {
1589     size_t i, n, m;
1590     mbedtls_mpi_uint u0, u1, *d;
1591 
1592     if( T->n < N->n + 1 || T->p == NULL )
1593         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1594 
1595     memset( T->p, 0, T->n * ciL );
1596 
1597     d = T->p;
1598     n = N->n;
1599     m = ( B->n < n ) ? B->n : n;
1600 
1601     for( i = 0; i < n; i++ )
1602     {
1603         /*
1604          * T = (T + u0*B + u1*N) / 2^biL
1605          */
1606         u0 = A->p[i];
1607         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1608 
1609         mpi_mul_hlp( m, B->p, d, u0 );
1610         mpi_mul_hlp( n, N->p, d, u1 );
1611 
1612         *d++ = u0; d[n + 1] = 0;
1613     }
1614 
1615     memcpy( A->p, d, ( n + 1 ) * ciL );
1616 
1617     if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )
1618         mpi_sub_hlp( n, N->p, A->p );
1619     else
1620         /* prevent timing attacks */
1621         mpi_sub_hlp( n, A->p, T->p );
1622 
1623     return( 0 );
1624 }
1625 
1626 /*
1627  * Montgomery reduction: A = A * R^-1 mod N
1628  */
1629 static int mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N, mbedtls_mpi_uint mm, const mbedtls_mpi *T )
1630 {
1631     mbedtls_mpi_uint z = 1;
1632     mbedtls_mpi U;
1633 
1634     U.n = U.s = (int) z;
1635     U.p = &z;
1636 
1637     return( mpi_montmul( A, &U, N, mm, T ) );
1638 }
1639 
1640 /*
1641  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1642  */
1643 int mbedtls_mpi_exp_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *E, const mbedtls_mpi *N, mbedtls_mpi *_RR )
1644 {
1645     int ret;
1646     size_t wbits, wsize, one = 1;
1647     size_t i, j, nblimbs;
1648     size_t bufsize, nbits;
1649     mbedtls_mpi_uint ei, mm, state;
1650     mbedtls_mpi RR, T, W[ 2 << MBEDTLS_MPI_WINDOW_SIZE ], Apos;
1651     int neg;
1652 
1653     if( mbedtls_mpi_cmp_int( N, 0 ) <= 0 || ( N->p[0] & 1 ) == 0 )
1654         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1655 
1656     if( mbedtls_mpi_cmp_int( E, 0 ) < 0 )
1657         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1658 
1659     /*
1660      * Init temps and window size
1661      */
1662     mpi_montg_init( &mm, N );
1663     mbedtls_mpi_init( &RR ); mbedtls_mpi_init( &T );
1664     mbedtls_mpi_init( &Apos );
1665     memset( W, 0, sizeof( W ) );
1666 
1667     i = mbedtls_mpi_bitlen( E );
1668 
1669     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1670             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1671 
1672     if( wsize > MBEDTLS_MPI_WINDOW_SIZE )
1673         wsize = MBEDTLS_MPI_WINDOW_SIZE;
1674 
1675     j = N->n + 1;
1676     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
1677     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[1],  j ) );
1678     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T, j * 2 ) );
1679 
1680     /*
1681      * Compensate for negative A (and correct at the end)
1682      */
1683     neg = ( A->s == -1 );
1684     if( neg )
1685     {
1686         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Apos, A ) );
1687         Apos.s = 1;
1688         A = &Apos;
1689     }
1690 
1691     /*
1692      * If 1st call, pre-compute R^2 mod N
1693      */
1694     if( _RR == NULL || _RR->p == NULL )
1695     {
1696         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &RR, 1 ) );
1697         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &RR, N->n * 2 * biL ) );
1698         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &RR, &RR, N ) );
1699 
1700         if( _RR != NULL )
1701             memcpy( _RR, &RR, sizeof( mbedtls_mpi ) );
1702     }
1703     else
1704         memcpy( &RR, _RR, sizeof( mbedtls_mpi ) );
1705 
1706     /*
1707      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1708      */
1709     if( mbedtls_mpi_cmp_mpi( A, N ) >= 0 )
1710         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &W[1], A, N ) );
1711     else
1712         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[1], A ) );
1713 
1714     MBEDTLS_MPI_CHK( mpi_montmul( &W[1], &RR, N, mm, &T ) );
1715 
1716     /*
1717      * X = R^2 * R^-1 mod N = R mod N
1718      */
1719     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &RR ) );
1720     MBEDTLS_MPI_CHK( mpi_montred( X, N, mm, &T ) );
1721 
1722     if( wsize > 1 )
1723     {
1724         /*
1725          * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1726          */
1727         j =  one << ( wsize - 1 );
1728 
1729         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[j], N->n + 1 ) );
1730         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[j], &W[1]    ) );
1731 
1732         for( i = 0; i < wsize - 1; i++ )
1733             MBEDTLS_MPI_CHK( mpi_montmul( &W[j], &W[j], N, mm, &T ) );
1734 
1735         /*
1736          * W[i] = W[i - 1] * W[1]
1737          */
1738         for( i = j + 1; i < ( one << wsize ); i++ )
1739         {
1740             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[i], N->n + 1 ) );
1741             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[i], &W[i - 1] ) );
1742 
1743             MBEDTLS_MPI_CHK( mpi_montmul( &W[i], &W[1], N, mm, &T ) );
1744         }
1745     }
1746 
1747     nblimbs = E->n;
1748     bufsize = 0;
1749     nbits   = 0;
1750     wbits   = 0;
1751     state   = 0;
1752 
1753     while( 1 )
1754     {
1755         if( bufsize == 0 )
1756         {
1757             if( nblimbs == 0 )
1758                 break;
1759 
1760             nblimbs--;
1761 
1762             bufsize = sizeof( mbedtls_mpi_uint ) << 3;
1763         }
1764 
1765         bufsize--;
1766 
1767         ei = (E->p[nblimbs] >> bufsize) & 1;
1768 
1769         /*
1770          * skip leading 0s
1771          */
1772         if( ei == 0 && state == 0 )
1773             continue;
1774 
1775         if( ei == 0 && state == 1 )
1776         {
1777             /*
1778              * out of window, square X
1779              */
1780             MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1781             continue;
1782         }
1783 
1784         /*
1785          * add ei to current window
1786          */
1787         state = 2;
1788 
1789         nbits++;
1790         wbits |= ( ei << ( wsize - nbits ) );
1791 
1792         if( nbits == wsize )
1793         {
1794             /*
1795              * X = X^wsize R^-1 mod N
1796              */
1797             for( i = 0; i < wsize; i++ )
1798                 MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1799 
1800             /*
1801              * X = X * W[wbits] R^-1 mod N
1802              */
1803             MBEDTLS_MPI_CHK( mpi_montmul( X, &W[wbits], N, mm, &T ) );
1804 
1805             state--;
1806             nbits = 0;
1807             wbits = 0;
1808         }
1809     }
1810 
1811     /*
1812      * process the remaining bits
1813      */
1814     for( i = 0; i < nbits; i++ )
1815     {
1816         MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1817 
1818         wbits <<= 1;
1819 
1820         if( ( wbits & ( one << wsize ) ) != 0 )
1821             MBEDTLS_MPI_CHK( mpi_montmul( X, &W[1], N, mm, &T ) );
1822     }
1823 
1824     /*
1825      * X = A^E * R * R^-1 mod N = A^E mod N
1826      */
1827     MBEDTLS_MPI_CHK( mpi_montred( X, N, mm, &T ) );
1828 
1829     if( neg && E->n != 0 && ( E->p[0] & 1 ) != 0 )
1830     {
1831         X->s = -1;
1832         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( X, N, X ) );
1833     }
1834 
1835 cleanup:
1836 
1837     for( i = ( one << ( wsize - 1 ) ); i < ( one << wsize ); i++ )
1838         mbedtls_mpi_free( &W[i] );
1839 
1840     mbedtls_mpi_free( &W[1] ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &Apos );
1841 
1842     if( _RR == NULL || _RR->p == NULL )
1843         mbedtls_mpi_free( &RR );
1844 
1845     return( ret );
1846 }
1847 
1848 /*
1849  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1850  */
1851 int mbedtls_mpi_gcd( mbedtls_mpi *G, const mbedtls_mpi *A, const mbedtls_mpi *B )
1852 {
1853     int ret;
1854     size_t lz, lzt;
1855     mbedtls_mpi TG, TA, TB;
1856 
1857     mbedtls_mpi_init( &TG ); mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
1858 
1859     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) );
1860     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
1861 
1862     lz = mbedtls_mpi_lsb( &TA );
1863     lzt = mbedtls_mpi_lsb( &TB );
1864 
1865     if( lzt < lz )
1866         lz = lzt;
1867 
1868     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, lz ) );
1869     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, lz ) );
1870 
1871     TA.s = TB.s = 1;
1872 
1873     while( mbedtls_mpi_cmp_int( &TA, 0 ) != 0 )
1874     {
1875         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, mbedtls_mpi_lsb( &TA ) ) );
1876         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, mbedtls_mpi_lsb( &TB ) ) );
1877 
1878         if( mbedtls_mpi_cmp_mpi( &TA, &TB ) >= 0 )
1879         {
1880             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TA, &TA, &TB ) );
1881             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, 1 ) );
1882         }
1883         else
1884         {
1885             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TB, &TB, &TA ) );
1886             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, 1 ) );
1887         }
1888     }
1889 
1890     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &TB, lz ) );
1891     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( G, &TB ) );
1892 
1893 cleanup:
1894 
1895     mbedtls_mpi_free( &TG ); mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TB );
1896 
1897     return( ret );
1898 }
1899 
1900 /*
1901  * Fill X with size bytes of random.
1902  *
1903  * Use a temporary bytes representation to make sure the result is the same
1904  * regardless of the platform endianness (useful when f_rng is actually
1905  * deterministic, eg for tests).
1906  */
1907 int mbedtls_mpi_fill_random( mbedtls_mpi *X, size_t size,
1908                      int (*f_rng)(void *, unsigned char *, size_t),
1909                      void *p_rng )
1910 {
1911     int ret;
1912     unsigned char buf[MBEDTLS_MPI_MAX_SIZE];
1913 
1914     if( size > MBEDTLS_MPI_MAX_SIZE )
1915         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1916 
1917     MBEDTLS_MPI_CHK( f_rng( p_rng, buf, size ) );
1918     MBEDTLS_MPI_CHK( mbedtls_mpi_read_binary( X, buf, size ) );
1919 
1920 cleanup:
1921     mbedtls_zeroize( buf, sizeof( buf ) );
1922     return( ret );
1923 }
1924 
1925 /*
1926  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1927  */
1928 int mbedtls_mpi_inv_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *N )
1929 {
1930     int ret;
1931     mbedtls_mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1932 
1933     if( mbedtls_mpi_cmp_int( N, 1 ) <= 0 )
1934         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1935 
1936     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TU ); mbedtls_mpi_init( &U1 ); mbedtls_mpi_init( &U2 );
1937     mbedtls_mpi_init( &G ); mbedtls_mpi_init( &TB ); mbedtls_mpi_init( &TV );
1938     mbedtls_mpi_init( &V1 ); mbedtls_mpi_init( &V2 );
1939 
1940     MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &G, A, N ) );
1941 
1942     if( mbedtls_mpi_cmp_int( &G, 1 ) != 0 )
1943     {
1944         ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
1945         goto cleanup;
1946     }
1947 
1948     MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &TA, A, N ) );
1949     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TU, &TA ) );
1950     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, N ) );
1951     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TV, N ) );
1952 
1953     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U1, 1 ) );
1954     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U2, 0 ) );
1955     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V1, 0 ) );
1956     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V2, 1 ) );
1957 
1958     do
1959     {
1960         while( ( TU.p[0] & 1 ) == 0 )
1961         {
1962             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TU, 1 ) );
1963 
1964             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1965             {
1966                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &U1, &U1, &TB ) );
1967                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &TA ) );
1968             }
1969 
1970             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U1, 1 ) );
1971             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U2, 1 ) );
1972         }
1973 
1974         while( ( TV.p[0] & 1 ) == 0 )
1975         {
1976             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TV, 1 ) );
1977 
1978             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1979             {
1980                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, &TB ) );
1981                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &TA ) );
1982             }
1983 
1984             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V1, 1 ) );
1985             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V2, 1 ) );
1986         }
1987 
1988         if( mbedtls_mpi_cmp_mpi( &TU, &TV ) >= 0 )
1989         {
1990             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TU, &TU, &TV ) );
1991             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U1, &U1, &V1 ) );
1992             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &V2 ) );
1993         }
1994         else
1995         {
1996             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TV, &TV, &TU ) );
1997             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, &U1 ) );
1998             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &U2 ) );
1999         }
2000     }
2001     while( mbedtls_mpi_cmp_int( &TU, 0 ) != 0 );
2002 
2003     while( mbedtls_mpi_cmp_int( &V1, 0 ) < 0 )
2004         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, N ) );
2005 
2006     while( mbedtls_mpi_cmp_mpi( &V1, N ) >= 0 )
2007         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, N ) );
2008 
2009     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &V1 ) );
2010 
2011 cleanup:
2012 
2013     mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TU ); mbedtls_mpi_free( &U1 ); mbedtls_mpi_free( &U2 );
2014     mbedtls_mpi_free( &G ); mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TV );
2015     mbedtls_mpi_free( &V1 ); mbedtls_mpi_free( &V2 );
2016 
2017     return( ret );
2018 }
2019 
2020 #if defined(MBEDTLS_GENPRIME)
2021 
2022 static const int small_prime[] =
2023 {
2024         3,    5,    7,   11,   13,   17,   19,   23,
2025        29,   31,   37,   41,   43,   47,   53,   59,
2026        61,   67,   71,   73,   79,   83,   89,   97,
2027       101,  103,  107,  109,  113,  127,  131,  137,
2028       139,  149,  151,  157,  163,  167,  173,  179,
2029       181,  191,  193,  197,  199,  211,  223,  227,
2030       229,  233,  239,  241,  251,  257,  263,  269,
2031       271,  277,  281,  283,  293,  307,  311,  313,
2032       317,  331,  337,  347,  349,  353,  359,  367,
2033       373,  379,  383,  389,  397,  401,  409,  419,
2034       421,  431,  433,  439,  443,  449,  457,  461,
2035       463,  467,  479,  487,  491,  499,  503,  509,
2036       521,  523,  541,  547,  557,  563,  569,  571,
2037       577,  587,  593,  599,  601,  607,  613,  617,
2038       619,  631,  641,  643,  647,  653,  659,  661,
2039       673,  677,  683,  691,  701,  709,  719,  727,
2040       733,  739,  743,  751,  757,  761,  769,  773,
2041       787,  797,  809,  811,  821,  823,  827,  829,
2042       839,  853,  857,  859,  863,  877,  881,  883,
2043       887,  907,  911,  919,  929,  937,  941,  947,
2044       953,  967,  971,  977,  983,  991,  997, -103
2045 };
2046 
2047 /*
2048  * Small divisors test (X must be positive)
2049  *
2050  * Return values:
2051  * 0: no small factor (possible prime, more tests needed)
2052  * 1: certain prime
2053  * MBEDTLS_ERR_MPI_NOT_ACCEPTABLE: certain non-prime
2054  * other negative: error
2055  */
2056 static int mpi_check_small_factors( const mbedtls_mpi *X )
2057 {
2058     int ret = 0;
2059     size_t i;
2060     mbedtls_mpi_uint r;
2061 
2062     if( ( X->p[0] & 1 ) == 0 )
2063         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2064 
2065     for( i = 0; small_prime[i] > 0; i++ )
2066     {
2067         if( mbedtls_mpi_cmp_int( X, small_prime[i] ) <= 0 )
2068             return( 1 );
2069 
2070         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, small_prime[i] ) );
2071 
2072         if( r == 0 )
2073             return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2074     }
2075 
2076 cleanup:
2077     return( ret );
2078 }
2079 
2080 /*
2081  * Miller-Rabin pseudo-primality test  (HAC 4.24)
2082  */
2083 static int mpi_miller_rabin( const mbedtls_mpi *X, size_t rounds,
2084                              int (*f_rng)(void *, unsigned char *, size_t),
2085                              void *p_rng )
2086 {
2087     int ret, count;
2088     size_t i, j, k, s;
2089     mbedtls_mpi W, R, T, A, RR;
2090 
2091     mbedtls_mpi_init( &W ); mbedtls_mpi_init( &R ); mbedtls_mpi_init( &T ); mbedtls_mpi_init( &A );
2092     mbedtls_mpi_init( &RR );
2093 
2094     /*
2095      * W = |X| - 1
2096      * R = W >> lsb( W )
2097      */
2098     MBEDTLS_MPI_CHK( mbedtls_mpi_sub_int( &W, X, 1 ) );
2099     s = mbedtls_mpi_lsb( &W );
2100     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &R, &W ) );
2101     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &R, s ) );
2102 
2103     i = mbedtls_mpi_bitlen( X );
2104 
2105     for( i = 0; i < rounds; i++ )
2106     {
2107         /*
2108          * pick a random A, 1 < A < |X| - 1
2109          */
2110         count = 0;
2111         do {
2112             MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
2113 
2114             j = mbedtls_mpi_bitlen( &A );
2115             k = mbedtls_mpi_bitlen( &W );
2116             if (j > k) {
2117                 A.p[A.n - 1] &= ( (mbedtls_mpi_uint) 1 << ( k - ( A.n - 1 ) * biL - 1 ) ) - 1;
2118             }
2119 
2120             if (count++ > 30) {
2121                 return MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2122             }
2123 
2124         } while ( mbedtls_mpi_cmp_mpi( &A, &W ) >= 0 ||
2125                   mbedtls_mpi_cmp_int( &A, 1 )  <= 0    );
2126 
2127         /*
2128          * A = A^R mod |X|
2129          */
2130         MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &A, &A, &R, X, &RR ) );
2131 
2132         if( mbedtls_mpi_cmp_mpi( &A, &W ) == 0 ||
2133             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2134             continue;
2135 
2136         j = 1;
2137         while( j < s && mbedtls_mpi_cmp_mpi( &A, &W ) != 0 )
2138         {
2139             /*
2140              * A = A * A mod |X|
2141              */
2142             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &T, &A, &A ) );
2143             MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &A, &T, X  ) );
2144 
2145             if( mbedtls_mpi_cmp_int( &A, 1 ) == 0 )
2146                 break;
2147 
2148             j++;
2149         }
2150 
2151         /*
2152          * not prime if A != |X| - 1 or A == 1
2153          */
2154         if( mbedtls_mpi_cmp_mpi( &A, &W ) != 0 ||
2155             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2156         {
2157             ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2158             break;
2159         }
2160     }
2161 
2162 cleanup:
2163     mbedtls_mpi_free( &W ); mbedtls_mpi_free( &R ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &A );
2164     mbedtls_mpi_free( &RR );
2165 
2166     return( ret );
2167 }
2168 
2169 /*
2170  * Pseudo-primality test: small factors, then Miller-Rabin
2171  */
2172 static int mpi_is_prime_internal( const mbedtls_mpi *X, int rounds,
2173                   int (*f_rng)(void *, unsigned char *, size_t),
2174                   void *p_rng )
2175 {
2176     int ret;
2177     mbedtls_mpi XX;
2178 
2179     XX.s = 1;
2180     XX.n = X->n;
2181     XX.p = X->p;
2182 
2183     if( mbedtls_mpi_cmp_int( &XX, 0 ) == 0 ||
2184         mbedtls_mpi_cmp_int( &XX, 1 ) == 0 )
2185         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2186 
2187     if( mbedtls_mpi_cmp_int( &XX, 2 ) == 0 )
2188         return( 0 );
2189 
2190     if( ( ret = mpi_check_small_factors( &XX ) ) != 0 )
2191     {
2192         if( ret == 1 )
2193             return( 0 );
2194 
2195         return( ret );
2196     }
2197 
2198     return( mpi_miller_rabin( &XX, rounds, f_rng, p_rng ) );
2199 }
2200 
2201 /*
2202  * Pseudo-primality test, error probability 2^-80
2203  */
2204 int mbedtls_mpi_is_prime( const mbedtls_mpi *X,
2205                   int (*f_rng)(void *, unsigned char *, size_t),
2206                   void *p_rng )
2207 {
2208     return mpi_is_prime_internal( X, 40, f_rng, p_rng );
2209 }
2210 
2211 /*
2212  * Prime number generation
2213  */
2214 int mbedtls_mpi_gen_prime( mbedtls_mpi *X, size_t nbits, int dh_flag,
2215                    int (*f_rng)(void *, unsigned char *, size_t),
2216                    void *p_rng )
2217 {
2218     int ret;
2219     size_t k, n;
2220     int rounds;
2221     mbedtls_mpi_uint r;
2222     mbedtls_mpi Y;
2223 
2224     if( nbits < 3 || nbits > MBEDTLS_MPI_MAX_BITS )
2225         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
2226 
2227     mbedtls_mpi_init( &Y );
2228 
2229     n = BITS_TO_LIMBS( nbits );
2230 
2231     /*
2232      * 2^-80 error probability, number of rounds chosen per HAC, table 4.4
2233      */
2234     rounds = ( ( nbits >= 1300 ) ?  2 : ( nbits >=  850 ) ?  3 :
2235                ( nbits >=  650 ) ?  4 : ( nbits >=  350 ) ?  8 :
2236                ( nbits >=  250 ) ? 12 : ( nbits >=  150 ) ? 18 : 27 );
2237 
2238     MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
2239 
2240     k = mbedtls_mpi_bitlen( X );
2241     if( k > nbits ) MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( X, k - nbits + 1 ) );
2242 
2243     mbedtls_mpi_set_bit( X, nbits-1, 1 );
2244 
2245     X->p[0] |= 1;
2246 
2247     if( dh_flag == 0 )
2248     {
2249         while( ( ret = mpi_is_prime_internal( X, rounds, f_rng, p_rng ) ) != 0 )
2250         {
2251             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2252                 goto cleanup;
2253 
2254             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 2 ) );
2255         }
2256     }
2257     else
2258     {
2259         /*
2260          * An necessary condition for Y and X = 2Y + 1 to be prime
2261          * is X = 2 mod 3 (which is equivalent to Y = 2 mod 3).
2262          * Make sure it is satisfied, while keeping X = 3 mod 4
2263          */
2264 
2265         X->p[0] |= 2;
2266 
2267         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, 3 ) );
2268         if( r == 0 )
2269             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 8 ) );
2270         else if( r == 1 )
2271             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 4 ) );
2272 
2273         /* Set Y = (X-1) / 2, which is X / 2 because X is odd */
2274         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, X ) );
2275         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, 1 ) );
2276 
2277         while( 1 )
2278         {
2279             /*
2280              * First, check small factors for X and Y
2281              * before doing Miller-Rabin on any of them
2282              */
2283             if( ( ret = mpi_check_small_factors(  X         ) ) == 0 &&
2284                 ( ret = mpi_check_small_factors( &Y         ) ) == 0 &&
2285                 ( ret = mpi_miller_rabin(  X, rounds, f_rng, p_rng  ) )
2286                                                                 == 0 &&
2287                 ( ret = mpi_miller_rabin( &Y, rounds, f_rng, p_rng  ) )
2288                                                                 == 0 )
2289             {
2290                 break;
2291             }
2292 
2293             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2294                 goto cleanup;
2295 
2296             /*
2297              * Next candidates. We want to preserve Y = (X-1) / 2 and
2298              * Y = 1 mod 2 and Y = 2 mod 3 (eq X = 3 mod 4 and X = 2 mod 3)
2299              * so up Y by 6 and X by 12.
2300              */
2301             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int(  X,  X, 12 ) );
2302             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( &Y, &Y, 6  ) );
2303         }
2304     }
2305 
2306 cleanup:
2307 
2308     mbedtls_mpi_free( &Y );
2309 
2310     return( ret );
2311 }
2312 
2313 #endif /* MBEDTLS_GENPRIME */
2314 
2315 #if defined(MBEDTLS_SELF_TEST)
2316 
2317 #define GCD_PAIR_COUNT  3
2318 
2319 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
2320 {
2321     { 693, 609, 21 },
2322     { 1764, 868, 28 },
2323     { 768454923, 542167814, 1 }
2324 };
2325 
2326 /*
2327  * Checkup routine
2328  */
2329 int mbedtls_mpi_self_test( int verbose )
2330 {
2331     int ret, i;
2332     mbedtls_mpi A, E, N, X, Y, U, V;
2333 
2334     mbedtls_mpi_init( &A ); mbedtls_mpi_init( &E ); mbedtls_mpi_init( &N ); mbedtls_mpi_init( &X );
2335     mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &U ); mbedtls_mpi_init( &V );
2336 
2337     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &A, 16,
2338         "EFE021C2645FD1DC586E69184AF4A31E" \
2339         "D5F53E93B5F123FA41680867BA110131" \
2340         "944FE7952E2517337780CB0DB80E61AA" \
2341         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
2342 
2343     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &E, 16,
2344         "B2E7EFD37075B9F03FF989C7C5051C20" \
2345         "34D2A323810251127E7BF8625A4F49A5" \
2346         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2347         "5B5C25763222FEFCCFC38B832366C29E" ) );
2348 
2349     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &N, 16,
2350         "0066A198186C18C10B2F5ED9B522752A" \
2351         "9830B69916E535C8F047518A889A43A5" \
2352         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2353 
2354     MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &X, &A, &N ) );
2355 
2356     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2357         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2358         "9E857EA95A03512E2BAE7391688D264A" \
2359         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2360         "8001B72E848A38CAE1C65F78E56ABDEF" \
2361         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2362         "ECF677152EF804370C1A305CAF3B5BF1" \
2363         "30879B56C61DE584A0F53A2447A51E" ) );
2364 
2365     if( verbose != 0 )
2366         mbedtls_printf( "  MPI test #1 (mul_mpi): " );
2367 
2368     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2369     {
2370         if( verbose != 0 )
2371             mbedtls_printf( "failed\n" );
2372 
2373         ret = 1;
2374         goto cleanup;
2375     }
2376 
2377     if( verbose != 0 )
2378         mbedtls_printf( "passed\n" );
2379 
2380     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( &X, &Y, &A, &N ) );
2381 
2382     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2383         "256567336059E52CAE22925474705F39A94" ) );
2384 
2385     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &V, 16,
2386         "6613F26162223DF488E9CD48CC132C7A" \
2387         "0AC93C701B001B092E4E5B9F73BCD27B" \
2388         "9EE50D0657C77F374E903CDFA4C642" ) );
2389 
2390     if( verbose != 0 )
2391         mbedtls_printf( "  MPI test #2 (div_mpi): " );
2392 
2393     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 ||
2394         mbedtls_mpi_cmp_mpi( &Y, &V ) != 0 )
2395     {
2396         if( verbose != 0 )
2397             mbedtls_printf( "failed\n" );
2398 
2399         ret = 1;
2400         goto cleanup;
2401     }
2402 
2403     if( verbose != 0 )
2404         mbedtls_printf( "passed\n" );
2405 
2406     MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2407 
2408     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2409         "36E139AEA55215609D2816998ED020BB" \
2410         "BD96C37890F65171D948E9BC7CBAA4D9" \
2411         "325D24D6A3C12710F10A09FA08AB87" ) );
2412 
2413     if( verbose != 0 )
2414         mbedtls_printf( "  MPI test #3 (exp_mod): " );
2415 
2416     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 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_inv_mod( &X, &A, &N ) );
2429 
2430     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2431         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2432         "C3DBA76456363A10869622EAC2DD84EC" \
2433         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2434 
2435     if( verbose != 0 )
2436         mbedtls_printf( "  MPI test #4 (inv_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     if( verbose != 0 )
2451         mbedtls_printf( "  MPI test #5 (simple gcd): " );
2452 
2453     for( i = 0; i < GCD_PAIR_COUNT; i++ )
2454     {
2455         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &X, gcd_pairs[i][0] ) );
2456         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Y, gcd_pairs[i][1] ) );
2457 
2458         MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &A, &X, &Y ) );
2459 
2460         if( mbedtls_mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2461         {
2462             if( verbose != 0 )
2463                 mbedtls_printf( "failed at %d\n", i );
2464 
2465             ret = 1;
2466             goto cleanup;
2467         }
2468     }
2469 
2470     if( verbose != 0 )
2471         mbedtls_printf( "passed\n" );
2472 
2473 cleanup:
2474 
2475     if( ret != 0 && verbose != 0 )
2476         mbedtls_printf( "Unexpected error, return code = %08X\n", ret );
2477 
2478     mbedtls_mpi_free( &A ); mbedtls_mpi_free( &E ); mbedtls_mpi_free( &N ); mbedtls_mpi_free( &X );
2479     mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &U ); mbedtls_mpi_free( &V );
2480 
2481     if( verbose != 0 )
2482         mbedtls_printf( "\n" );
2483 
2484     return( ret );
2485 }
2486 
2487 #endif /* MBEDTLS_SELF_TEST */
2488 
2489 #endif /* MBEDTLS_BIGNUM_C */
2490