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