xref: /optee_os/lib/libmbedtls/mbedtls/library/bignum.c (revision 62f21181c547da3bd098908300e5699e9ae5cca9)
1 // SPDX-License-Identifier: Apache-2.0
2 /*
3  *  Multi-precision integer library
4  *
5  *  Copyright (C) 2006-2015, ARM Limited, All Rights Reserved
6  *
7  *  Licensed under the Apache License, Version 2.0 (the "License"); you may
8  *  not use this file except in compliance with the License.
9  *  You may obtain a copy of the License at
10  *
11  *  http://www.apache.org/licenses/LICENSE-2.0
12  *
13  *  Unless required by applicable law or agreed to in writing, software
14  *  distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
15  *  WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16  *  See the License for the specific language governing permissions and
17  *  limitations under the License.
18  *
19  *  This file is part of mbed TLS (https://tls.mbed.org)
20  */
21 
22 /*
23  *  The following sources were referenced in the design of this Multi-precision
24  *  Integer library:
25  *
26  *  [1] Handbook of Applied Cryptography - 1997
27  *      Menezes, van Oorschot and Vanstone
28  *
29  *  [2] Multi-Precision Math
30  *      Tom St Denis
31  *      https://github.com/libtom/libtommath/blob/develop/tommath.pdf
32  *
33  *  [3] GNU Multi-Precision Arithmetic Library
34  *      https://gmplib.org/manual/index.html
35  *
36  */
37 
38 #if !defined(MBEDTLS_CONFIG_FILE)
39 #include "mbedtls/config.h"
40 #else
41 #include MBEDTLS_CONFIG_FILE
42 #endif
43 
44 #if defined(MBEDTLS_BIGNUM_C)
45 
46 #include "mbedtls/bignum.h"
47 #include "mbedtls/bn_mul.h"
48 
49 #include <string.h>
50 
51 #if defined(MBEDTLS_PLATFORM_C)
52 #include "mbedtls/platform.h"
53 #else
54 #include <stdio.h>
55 #include <stdlib.h>
56 #define mbedtls_printf     printf
57 #define mbedtls_calloc    calloc
58 #define mbedtls_free       free
59 #endif
60 
61 /* Implementation that should never be optimized out by the compiler */
62 static void mbedtls_mpi_zeroize( mbedtls_mpi_uint *v, size_t n ) {
63     volatile mbedtls_mpi_uint *p = v; while( n-- ) *p++ = 0;
64 }
65 
66 #define ciL    (sizeof(mbedtls_mpi_uint))         /* chars in limb  */
67 #define biL    (ciL << 3)               /* bits  in limb  */
68 #define biH    (ciL << 2)               /* half limb size */
69 
70 #define MPI_SIZE_T_MAX  ( (size_t) -1 ) /* SIZE_T_MAX is not standard */
71 
72 /*
73  * Convert between bits/chars and number of limbs
74  * Divide first in order to avoid potential overflows
75  */
76 #define BITS_TO_LIMBS(i)  ( (i) / biL + ( (i) % biL != 0 ) )
77 #define CHARS_TO_LIMBS(i) ( (i) / ciL + ( (i) % ciL != 0 ) )
78 
79 /*
80  * Initialize one MPI
81  */
82 void mbedtls_mpi_init( mbedtls_mpi *X )
83 {
84     if( X == NULL )
85         return;
86 
87     X->s = 1;
88     X->n = 0;
89     X->p = NULL;
90 }
91 
92 /*
93  * Unallocate one MPI
94  */
95 void mbedtls_mpi_free( mbedtls_mpi *X )
96 {
97     if( X == NULL )
98         return;
99 
100     if( X->p != NULL )
101     {
102         mbedtls_mpi_zeroize( X->p, X->n );
103         mbedtls_free( X->p );
104     }
105 
106     X->s = 1;
107     X->n = 0;
108     X->p = NULL;
109 }
110 
111 /*
112  * Enlarge to the specified number of limbs
113  */
114 int mbedtls_mpi_grow( mbedtls_mpi *X, size_t nblimbs )
115 {
116     mbedtls_mpi_uint *p;
117 
118     if( nblimbs > MBEDTLS_MPI_MAX_LIMBS )
119         return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
120 
121     if( X->n < nblimbs )
122     {
123         if( ( p = (mbedtls_mpi_uint*)mbedtls_calloc( nblimbs, ciL ) ) == NULL )
124             return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
125 
126         if( X->p != NULL )
127         {
128             memcpy( p, X->p, X->n * ciL );
129             mbedtls_mpi_zeroize( X->p, X->n );
130             mbedtls_free( X->p );
131         }
132 
133         X->n = nblimbs;
134         X->p = p;
135     }
136 
137     return( 0 );
138 }
139 
140 /*
141  * Resize down as much as possible,
142  * while keeping at least the specified number of limbs
143  */
144 int mbedtls_mpi_shrink( mbedtls_mpi *X, size_t nblimbs )
145 {
146     mbedtls_mpi_uint *p;
147     size_t i;
148 
149     /* Actually resize up in this case */
150     if( X->n <= nblimbs )
151         return( mbedtls_mpi_grow( X, nblimbs ) );
152 
153     for( i = X->n - 1; i > 0; i-- )
154         if( X->p[i] != 0 )
155             break;
156     i++;
157 
158     if( i < nblimbs )
159         i = nblimbs;
160 
161     if( ( p = (mbedtls_mpi_uint*)mbedtls_calloc( i, ciL ) ) == NULL )
162         return( MBEDTLS_ERR_MPI_ALLOC_FAILED );
163 
164     if( X->p != NULL )
165     {
166         memcpy( p, X->p, i * ciL );
167         mbedtls_mpi_zeroize( X->p, X->n );
168         mbedtls_free( X->p );
169     }
170 
171     X->n = i;
172     X->p = p;
173 
174     return( 0 );
175 }
176 
177 /*
178  * Copy the contents of Y into X
179  */
180 int mbedtls_mpi_copy( mbedtls_mpi *X, const mbedtls_mpi *Y )
181 {
182     int ret;
183     size_t i;
184 
185     if( X == Y )
186         return( 0 );
187 
188     if( Y->p == NULL )
189     {
190         mbedtls_mpi_free( X );
191         return( 0 );
192     }
193 
194     for( i = Y->n - 1; i > 0; i-- )
195         if( Y->p[i] != 0 )
196             break;
197     i++;
198 
199     X->s = Y->s;
200 
201     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i ) );
202 
203     memset( X->p, 0, X->n * ciL );
204     memcpy( X->p, Y->p, i * ciL );
205 
206 cleanup:
207 
208     return( ret );
209 }
210 
211 /*
212  * Swap the contents of X and Y
213  */
214 void mbedtls_mpi_swap( mbedtls_mpi *X, mbedtls_mpi *Y )
215 {
216     mbedtls_mpi T;
217 
218     memcpy( &T,  X, sizeof( mbedtls_mpi ) );
219     memcpy(  X,  Y, sizeof( mbedtls_mpi ) );
220     memcpy(  Y, &T, sizeof( mbedtls_mpi ) );
221 }
222 
223 /*
224  * Conditionally assign X = Y, without leaking information
225  * about whether the assignment was made or not.
226  * (Leaking information about the respective sizes of X and Y is ok however.)
227  */
228 int mbedtls_mpi_safe_cond_assign( mbedtls_mpi *X, const mbedtls_mpi *Y, unsigned char assign )
229 {
230     int ret = 0;
231     size_t i;
232 
233     /* make sure assign is 0 or 1 in a time-constant manner */
234     assign = (assign | (unsigned char)-assign) >> 7;
235 
236     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, Y->n ) );
237 
238     X->s = X->s * ( 1 - assign ) + Y->s * assign;
239 
240     for( i = 0; i < Y->n; i++ )
241         X->p[i] = X->p[i] * ( 1 - assign ) + Y->p[i] * assign;
242 
243     for( ; i < X->n; i++ )
244         X->p[i] *= ( 1 - assign );
245 
246 cleanup:
247     return( ret );
248 }
249 
250 /*
251  * Conditionally swap X and Y, without leaking information
252  * about whether the swap was made or not.
253  * Here it is not ok to simply swap the pointers, which whould lead to
254  * different memory access patterns when X and Y are used afterwards.
255  */
256 int mbedtls_mpi_safe_cond_swap( mbedtls_mpi *X, mbedtls_mpi *Y, unsigned char swap )
257 {
258     int ret, s;
259     size_t i;
260     mbedtls_mpi_uint tmp;
261 
262     if( X == Y )
263         return( 0 );
264 
265     /* make sure swap is 0 or 1 in a time-constant manner */
266     swap = (swap | (unsigned char)-swap) >> 7;
267 
268     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, Y->n ) );
269     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( Y, X->n ) );
270 
271     s = X->s;
272     X->s = X->s * ( 1 - swap ) + Y->s * swap;
273     Y->s = Y->s * ( 1 - swap ) +    s * swap;
274 
275 
276     for( i = 0; i < X->n; i++ )
277     {
278         tmp = X->p[i];
279         X->p[i] = X->p[i] * ( 1 - swap ) + Y->p[i] * swap;
280         Y->p[i] = Y->p[i] * ( 1 - swap ) +     tmp * swap;
281     }
282 
283 cleanup:
284     return( ret );
285 }
286 
287 /*
288  * Set value from integer
289  */
290 int mbedtls_mpi_lset( mbedtls_mpi *X, mbedtls_mpi_sint z )
291 {
292     int ret;
293 
294     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, 1 ) );
295     memset( X->p, 0, X->n * ciL );
296 
297     X->p[0] = ( z < 0 ) ? -z : z;
298     X->s    = ( z < 0 ) ? -1 : 1;
299 
300 cleanup:
301 
302     return( ret );
303 }
304 
305 /*
306  * Get a specific bit
307  */
308 int mbedtls_mpi_get_bit( const mbedtls_mpi *X, size_t pos )
309 {
310     if( X->n * biL <= pos )
311         return( 0 );
312 
313     return( ( X->p[pos / biL] >> ( pos % biL ) ) & 0x01 );
314 }
315 
316 /*
317  * Set a bit to a specific value of 0 or 1
318  */
319 int mbedtls_mpi_set_bit( mbedtls_mpi *X, size_t pos, unsigned char val )
320 {
321     int ret = 0;
322     size_t off = pos / biL;
323     size_t idx = pos % biL;
324 
325     if( val != 0 && val != 1 )
326         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
327 
328     if( X->n * biL <= pos )
329     {
330         if( val == 0 )
331             return( 0 );
332 
333         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, off + 1 ) );
334     }
335 
336     X->p[off] &= ~( (mbedtls_mpi_uint) 0x01 << idx );
337     X->p[off] |= (mbedtls_mpi_uint) val << idx;
338 
339 cleanup:
340 
341     return( ret );
342 }
343 
344 /*
345  * Return the number of less significant zero-bits
346  */
347 size_t mbedtls_mpi_lsb( const mbedtls_mpi *X )
348 {
349     size_t i, j, count = 0;
350 
351     for( i = 0; i < X->n; i++ )
352         for( j = 0; j < biL; j++, count++ )
353             if( ( ( X->p[i] >> j ) & 1 ) != 0 )
354                 return( count );
355 
356     return( 0 );
357 }
358 
359 /*
360  * Count leading zero bits in a given integer
361  */
362 static size_t mbedtls_clz( const mbedtls_mpi_uint x )
363 {
364     size_t j;
365     mbedtls_mpi_uint mask = (mbedtls_mpi_uint) 1 << (biL - 1);
366 
367     for( j = 0; j < biL; j++ )
368     {
369         if( x & mask ) break;
370 
371         mask >>= 1;
372     }
373 
374     return j;
375 }
376 
377 /*
378  * Return the number of bits
379  */
380 size_t mbedtls_mpi_bitlen( const mbedtls_mpi *X )
381 {
382     size_t i, j;
383 
384     if( X->n == 0 )
385         return( 0 );
386 
387     for( i = X->n - 1; i > 0; i-- )
388         if( X->p[i] != 0 )
389             break;
390 
391     j = biL - mbedtls_clz( X->p[i] );
392 
393     return( ( i * biL ) + j );
394 }
395 
396 /*
397  * Return the total size in bytes
398  */
399 size_t mbedtls_mpi_size( const mbedtls_mpi *X )
400 {
401     return( ( mbedtls_mpi_bitlen( X ) + 7 ) >> 3 );
402 }
403 
404 /*
405  * Convert an ASCII character to digit value
406  */
407 static int mpi_get_digit( mbedtls_mpi_uint *d, int radix, char c )
408 {
409     *d = 255;
410 
411     if( c >= 0x30 && c <= 0x39 ) *d = c - 0x30;
412     if( c >= 0x41 && c <= 0x46 ) *d = c - 0x37;
413     if( c >= 0x61 && c <= 0x66 ) *d = c - 0x57;
414 
415     if( *d >= (mbedtls_mpi_uint) radix )
416         return( MBEDTLS_ERR_MPI_INVALID_CHARACTER );
417 
418     return( 0 );
419 }
420 
421 /*
422  * Import from an ASCII string
423  */
424 int mbedtls_mpi_read_string( mbedtls_mpi *X, int radix, const char *s )
425 {
426     int ret;
427     size_t i, j, slen, n;
428     mbedtls_mpi_uint d;
429     mbedtls_mpi T;
430 
431     if( radix < 2 || radix > 16 )
432         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
433 
434     mbedtls_mpi_init( &T );
435 
436     slen = strlen( s );
437 
438     if( radix == 16 )
439     {
440         if( slen > MPI_SIZE_T_MAX >> 2 )
441             return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
442 
443         n = BITS_TO_LIMBS( slen << 2 );
444 
445         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, n ) );
446         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
447 
448         for( i = slen, j = 0; i > 0; i--, j++ )
449         {
450             if( i == 1 && s[i - 1] == '-' )
451             {
452                 X->s = -1;
453                 break;
454             }
455 
456             MBEDTLS_MPI_CHK( mpi_get_digit( &d, radix, s[i - 1] ) );
457             X->p[j / ( 2 * ciL )] |= d << ( ( j % ( 2 * ciL ) ) << 2 );
458         }
459     }
460     else
461     {
462         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
463 
464         for( i = 0; i < slen; i++ )
465         {
466             if( i == 0 && s[i] == '-' )
467             {
468                 X->s = -1;
469                 continue;
470             }
471 
472             MBEDTLS_MPI_CHK( mpi_get_digit( &d, radix, s[i] ) );
473             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T, X, radix ) );
474 
475             if( X->s == 1 )
476             {
477                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, &T, d ) );
478             }
479             else
480             {
481                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_int( X, &T, d ) );
482             }
483         }
484     }
485 
486 cleanup:
487 
488     mbedtls_mpi_free( &T );
489 
490     return( ret );
491 }
492 
493 /*
494  * Helper to write the digits high-order first
495  */
496 static int mpi_write_hlp( mbedtls_mpi *X, int radix, char **p )
497 {
498     int ret;
499     mbedtls_mpi_uint r;
500 
501     if( radix < 2 || radix > 16 )
502         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
503 
504     MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, radix ) );
505     MBEDTLS_MPI_CHK( mbedtls_mpi_div_int( X, NULL, X, radix ) );
506 
507     if( mbedtls_mpi_cmp_int( X, 0 ) != 0 )
508         MBEDTLS_MPI_CHK( mpi_write_hlp( X, radix, p ) );
509 
510     if( r < 10 )
511         *(*p)++ = (char)( r + 0x30 );
512     else
513         *(*p)++ = (char)( r + 0x37 );
514 
515 cleanup:
516 
517     return( ret );
518 }
519 
520 /*
521  * Export into an ASCII string
522  */
523 int mbedtls_mpi_write_string( const mbedtls_mpi *X, int radix,
524                               char *buf, size_t buflen, size_t *olen )
525 {
526     int ret = 0;
527     size_t n;
528     char *p;
529     mbedtls_mpi T;
530 
531     if( radix < 2 || radix > 16 )
532         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
533 
534     n = mbedtls_mpi_bitlen( X );
535     if( radix >=  4 ) n >>= 1;
536     if( radix >= 16 ) n >>= 1;
537     /*
538      * Round up the buffer length to an even value to ensure that there is
539      * enough room for hexadecimal values that can be represented in an odd
540      * number of digits.
541      */
542     n += 3 + ( ( n + 1 ) & 1 );
543 
544     if( buflen < n )
545     {
546         *olen = n;
547         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
548     }
549 
550     p = buf;
551     mbedtls_mpi_init( &T );
552 
553     if( X->s == -1 )
554         *p++ = '-';
555 
556     if( radix == 16 )
557     {
558         int c;
559         size_t i, j, k;
560 
561         for( i = X->n, k = 0; i > 0; i-- )
562         {
563             for( j = ciL; j > 0; j-- )
564             {
565                 c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
566 
567                 if( c == 0 && k == 0 && ( i + j ) != 2 )
568                     continue;
569 
570                 *(p++) = "0123456789ABCDEF" [c / 16];
571                 *(p++) = "0123456789ABCDEF" [c % 16];
572                 k = 1;
573             }
574         }
575     }
576     else
577     {
578         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T, X ) );
579 
580         if( T.s == -1 )
581             T.s = 1;
582 
583         MBEDTLS_MPI_CHK( mpi_write_hlp( &T, radix, &p ) );
584     }
585 
586     *p++ = '\0';
587     *olen = p - buf;
588 
589 cleanup:
590 
591     mbedtls_mpi_free( &T );
592 
593     return( ret );
594 }
595 
596 #if defined(MBEDTLS_FS_IO)
597 /*
598  * Read X from an opened file
599  */
600 int mbedtls_mpi_read_file( mbedtls_mpi *X, int radix, FILE *fin )
601 {
602     mbedtls_mpi_uint d;
603     size_t slen;
604     char *p;
605     /*
606      * Buffer should have space for (short) label and decimal formatted MPI,
607      * newline characters and '\0'
608      */
609     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
610 
611     memset( s, 0, sizeof( s ) );
612     if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
613         return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
614 
615     slen = strlen( s );
616     if( slen == sizeof( s ) - 2 )
617         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
618 
619     if( slen > 0 && s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
620     if( slen > 0 && s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }
621 
622     p = s + slen;
623     while( p-- > s )
624         if( mpi_get_digit( &d, radix, *p ) != 0 )
625             break;
626 
627     return( mbedtls_mpi_read_string( X, radix, p + 1 ) );
628 }
629 
630 /*
631  * Write X into an opened file (or stdout if fout == NULL)
632  */
633 int mbedtls_mpi_write_file( const char *p, const mbedtls_mpi *X, int radix, FILE *fout )
634 {
635     int ret;
636     size_t n, slen, plen;
637     /*
638      * Buffer should have space for (short) label and decimal formatted MPI,
639      * newline characters and '\0'
640      */
641     char s[ MBEDTLS_MPI_RW_BUFFER_SIZE ];
642 
643     memset( s, 0, sizeof( s ) );
644 
645     MBEDTLS_MPI_CHK( mbedtls_mpi_write_string( X, radix, s, sizeof( s ) - 2, &n ) );
646 
647     if( p == NULL ) p = "";
648 
649     plen = strlen( p );
650     slen = strlen( s );
651     s[slen++] = '\r';
652     s[slen++] = '\n';
653 
654     if( fout != NULL )
655     {
656         if( fwrite( p, 1, plen, fout ) != plen ||
657             fwrite( s, 1, slen, fout ) != slen )
658             return( MBEDTLS_ERR_MPI_FILE_IO_ERROR );
659     }
660     else
661         mbedtls_printf( "%s%s", p, s );
662 
663 cleanup:
664 
665     return( ret );
666 }
667 #endif /* MBEDTLS_FS_IO */
668 
669 /*
670  * Import X from unsigned binary data, big endian
671  */
672 int mbedtls_mpi_read_binary( mbedtls_mpi *X, const unsigned char *buf, size_t buflen )
673 {
674     int ret;
675     size_t i, j, n;
676 
677     for( n = 0; n < buflen; n++ )
678         if( buf[n] != 0 )
679             break;
680 
681     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, CHARS_TO_LIMBS( buflen - n ) ) );
682     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
683 
684     for( i = buflen, j = 0; i > n; i--, j++ )
685         X->p[j / ciL] |= ((mbedtls_mpi_uint) buf[i - 1]) << ((j % ciL) << 3);
686 
687 cleanup:
688 
689     return( ret );
690 }
691 
692 /*
693  * Export X into unsigned binary data, big endian
694  */
695 int mbedtls_mpi_write_binary( const mbedtls_mpi *X, unsigned char *buf, size_t buflen )
696 {
697     size_t i, j, n;
698 
699     n = mbedtls_mpi_size( X );
700 
701     if( buflen < n )
702         return( MBEDTLS_ERR_MPI_BUFFER_TOO_SMALL );
703 
704     memset( buf, 0, buflen );
705 
706     for( i = buflen - 1, j = 0; n > 0; i--, j++, n-- )
707         buf[i] = (unsigned char)( X->p[j / ciL] >> ((j % ciL) << 3) );
708 
709     return( 0 );
710 }
711 
712 /*
713  * Left-shift: X <<= count
714  */
715 int mbedtls_mpi_shift_l( mbedtls_mpi *X, size_t count )
716 {
717     int ret;
718     size_t i, v0, t1;
719     mbedtls_mpi_uint r0 = 0, r1;
720 
721     v0 = count / (biL    );
722     t1 = count & (biL - 1);
723 
724     i = mbedtls_mpi_bitlen( X ) + count;
725 
726     if( X->n * biL < i )
727         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, BITS_TO_LIMBS( i ) ) );
728 
729     ret = 0;
730 
731     /*
732      * shift by count / limb_size
733      */
734     if( v0 > 0 )
735     {
736         for( i = X->n; i > v0; i-- )
737             X->p[i - 1] = X->p[i - v0 - 1];
738 
739         for( ; i > 0; i-- )
740             X->p[i - 1] = 0;
741     }
742 
743     /*
744      * shift by count % limb_size
745      */
746     if( t1 > 0 )
747     {
748         for( i = v0; i < X->n; i++ )
749         {
750             r1 = X->p[i] >> (biL - t1);
751             X->p[i] <<= t1;
752             X->p[i] |= r0;
753             r0 = r1;
754         }
755     }
756 
757 cleanup:
758 
759     return( ret );
760 }
761 
762 /*
763  * Right-shift: X >>= count
764  */
765 int mbedtls_mpi_shift_r( mbedtls_mpi *X, size_t count )
766 {
767     size_t i, v0, v1;
768     mbedtls_mpi_uint r0 = 0, r1;
769 
770     v0 = count /  biL;
771     v1 = count & (biL - 1);
772 
773     if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
774         return mbedtls_mpi_lset( X, 0 );
775 
776     /*
777      * shift by count / limb_size
778      */
779     if( v0 > 0 )
780     {
781         for( i = 0; i < X->n - v0; i++ )
782             X->p[i] = X->p[i + v0];
783 
784         for( ; i < X->n; i++ )
785             X->p[i] = 0;
786     }
787 
788     /*
789      * shift by count % limb_size
790      */
791     if( v1 > 0 )
792     {
793         for( i = X->n; i > 0; i-- )
794         {
795             r1 = X->p[i - 1] << (biL - v1);
796             X->p[i - 1] >>= v1;
797             X->p[i - 1] |= r0;
798             r0 = r1;
799         }
800     }
801 
802     return( 0 );
803 }
804 
805 /*
806  * Compare unsigned values
807  */
808 int mbedtls_mpi_cmp_abs( const mbedtls_mpi *X, const mbedtls_mpi *Y )
809 {
810     size_t i, j;
811 
812     for( i = X->n; i > 0; i-- )
813         if( X->p[i - 1] != 0 )
814             break;
815 
816     for( j = Y->n; j > 0; j-- )
817         if( Y->p[j - 1] != 0 )
818             break;
819 
820     if( i == 0 && j == 0 )
821         return( 0 );
822 
823     if( i > j ) return(  1 );
824     if( j > i ) return( -1 );
825 
826     for( ; i > 0; i-- )
827     {
828         if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
829         if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
830     }
831 
832     return( 0 );
833 }
834 
835 /*
836  * Compare signed values
837  */
838 int mbedtls_mpi_cmp_mpi( const mbedtls_mpi *X, const mbedtls_mpi *Y )
839 {
840     size_t i, j;
841 
842     for( i = X->n; i > 0; i-- )
843         if( X->p[i - 1] != 0 )
844             break;
845 
846     for( j = Y->n; j > 0; j-- )
847         if( Y->p[j - 1] != 0 )
848             break;
849 
850     if( i == 0 && j == 0 )
851         return( 0 );
852 
853     if( i > j ) return(  X->s );
854     if( j > i ) return( -Y->s );
855 
856     if( X->s > 0 && Y->s < 0 ) return(  1 );
857     if( Y->s > 0 && X->s < 0 ) return( -1 );
858 
859     for( ; i > 0; i-- )
860     {
861         if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
862         if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
863     }
864 
865     return( 0 );
866 }
867 
868 /*
869  * Compare signed values
870  */
871 int mbedtls_mpi_cmp_int( const mbedtls_mpi *X, mbedtls_mpi_sint z )
872 {
873     mbedtls_mpi Y;
874     mbedtls_mpi_uint p[1];
875 
876     *p  = ( z < 0 ) ? -z : z;
877     Y.s = ( z < 0 ) ? -1 : 1;
878     Y.n = 1;
879     Y.p = p;
880 
881     return( mbedtls_mpi_cmp_mpi( X, &Y ) );
882 }
883 
884 /*
885  * Unsigned addition: X = |A| + |B|  (HAC 14.7)
886  */
887 int mbedtls_mpi_add_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
888 {
889     int ret;
890     size_t i, j;
891     mbedtls_mpi_uint *o, *p, c, tmp;
892 
893     if( X == B )
894     {
895         const mbedtls_mpi *T = A; A = X; B = T;
896     }
897 
898     if( X != A )
899         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
900 
901     /*
902      * X should always be positive as a result of unsigned additions.
903      */
904     X->s = 1;
905 
906     for( j = B->n; j > 0; j-- )
907         if( B->p[j - 1] != 0 )
908             break;
909 
910     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
911 
912     o = B->p; p = X->p; c = 0;
913 
914     /*
915      * tmp is used because it might happen that p == o
916      */
917     for( i = 0; i < j; i++, o++, p++ )
918     {
919         tmp= *o;
920         *p +=  c; c  = ( *p <  c );
921         *p += tmp; c += ( *p < tmp );
922     }
923 
924     while( c != 0 )
925     {
926         if( i >= X->n )
927         {
928             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + 1 ) );
929             p = X->p + i;
930         }
931 
932         *p += c; c = ( *p < c ); i++; p++;
933     }
934 
935 cleanup:
936 
937     return( ret );
938 }
939 
940 /*
941  * Helper for mbedtls_mpi subtraction
942  */
943 static void mpi_sub_hlp( size_t n, mbedtls_mpi_uint *s, mbedtls_mpi_uint *d )
944 {
945     size_t i;
946     mbedtls_mpi_uint c, z;
947 
948     for( i = c = 0; i < n; i++, s++, d++ )
949     {
950         z = ( *d <  c );     *d -=  c;
951         c = ( *d < *s ) + z; *d -= *s;
952     }
953 
954     while( c != 0 )
955     {
956         z = ( *d < c ); *d -= c;
957         c = z; i++; d++;
958     }
959 }
960 
961 /*
962  * Unsigned subtraction: X = |A| - |B|  (HAC 14.9)
963  */
964 int mbedtls_mpi_sub_abs( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
965 {
966     mbedtls_mpi TB;
967     int ret;
968     size_t n;
969 
970     if( mbedtls_mpi_cmp_abs( A, B ) < 0 )
971         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
972 
973     mbedtls_mpi_init( &TB );
974 
975     if( X == B )
976     {
977         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
978         B = &TB;
979     }
980 
981     if( X != A )
982         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, A ) );
983 
984     /*
985      * X should always be positive as a result of unsigned subtractions.
986      */
987     X->s = 1;
988 
989     ret = 0;
990 
991     for( n = B->n; n > 0; n-- )
992         if( B->p[n - 1] != 0 )
993             break;
994 
995     mpi_sub_hlp( n, B->p, X->p );
996 
997 cleanup:
998 
999     mbedtls_mpi_free( &TB );
1000 
1001     return( ret );
1002 }
1003 
1004 /*
1005  * Signed addition: X = A + B
1006  */
1007 int mbedtls_mpi_add_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1008 {
1009     int ret, s = A->s;
1010 
1011     if( A->s * B->s < 0 )
1012     {
1013         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1014         {
1015             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1016             X->s =  s;
1017         }
1018         else
1019         {
1020             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1021             X->s = -s;
1022         }
1023     }
1024     else
1025     {
1026         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1027         X->s = s;
1028     }
1029 
1030 cleanup:
1031 
1032     return( ret );
1033 }
1034 
1035 /*
1036  * Signed subtraction: X = A - B
1037  */
1038 int mbedtls_mpi_sub_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1039 {
1040     int ret, s = A->s;
1041 
1042     if( A->s * B->s > 0 )
1043     {
1044         if( mbedtls_mpi_cmp_abs( A, B ) >= 0 )
1045         {
1046             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, A, B ) );
1047             X->s =  s;
1048         }
1049         else
1050         {
1051             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( X, B, A ) );
1052             X->s = -s;
1053         }
1054     }
1055     else
1056     {
1057         MBEDTLS_MPI_CHK( mbedtls_mpi_add_abs( X, A, B ) );
1058         X->s = s;
1059     }
1060 
1061 cleanup:
1062 
1063     return( ret );
1064 }
1065 
1066 /*
1067  * Signed addition: X = A + b
1068  */
1069 int mbedtls_mpi_add_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1070 {
1071     mbedtls_mpi _B;
1072     mbedtls_mpi_uint p[1];
1073 
1074     p[0] = ( b < 0 ) ? -b : b;
1075     _B.s = ( b < 0 ) ? -1 : 1;
1076     _B.n = 1;
1077     _B.p = p;
1078 
1079     return( mbedtls_mpi_add_mpi( X, A, &_B ) );
1080 }
1081 
1082 /*
1083  * Signed subtraction: X = A - b
1084  */
1085 int mbedtls_mpi_sub_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1086 {
1087     mbedtls_mpi _B;
1088     mbedtls_mpi_uint p[1];
1089 
1090     p[0] = ( b < 0 ) ? -b : b;
1091     _B.s = ( b < 0 ) ? -1 : 1;
1092     _B.n = 1;
1093     _B.p = p;
1094 
1095     return( mbedtls_mpi_sub_mpi( X, A, &_B ) );
1096 }
1097 
1098 /*
1099  * Helper for mbedtls_mpi multiplication
1100  */
1101 static
1102 #if defined(__APPLE__) && defined(__arm__)
1103 /*
1104  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
1105  * appears to need this to prevent bad ARM code generation at -O3.
1106  */
1107 __attribute__ ((noinline))
1108 #endif
1109 void mpi_mul_hlp( size_t i, mbedtls_mpi_uint *s, mbedtls_mpi_uint *d, mbedtls_mpi_uint b )
1110 {
1111     mbedtls_mpi_uint c = 0, t = 0;
1112 
1113 #if defined(MULADDC_HUIT)
1114     for( ; i >= 8; i -= 8 )
1115     {
1116         MULADDC_INIT
1117         MULADDC_HUIT
1118         MULADDC_STOP
1119     }
1120 
1121     for( ; i > 0; i-- )
1122     {
1123         MULADDC_INIT
1124         MULADDC_CORE
1125         MULADDC_STOP
1126     }
1127 #else /* MULADDC_HUIT */
1128     for( ; i >= 16; i -= 16 )
1129     {
1130         MULADDC_INIT
1131         MULADDC_CORE   MULADDC_CORE
1132         MULADDC_CORE   MULADDC_CORE
1133         MULADDC_CORE   MULADDC_CORE
1134         MULADDC_CORE   MULADDC_CORE
1135 
1136         MULADDC_CORE   MULADDC_CORE
1137         MULADDC_CORE   MULADDC_CORE
1138         MULADDC_CORE   MULADDC_CORE
1139         MULADDC_CORE   MULADDC_CORE
1140         MULADDC_STOP
1141     }
1142 
1143     for( ; i >= 8; i -= 8 )
1144     {
1145         MULADDC_INIT
1146         MULADDC_CORE   MULADDC_CORE
1147         MULADDC_CORE   MULADDC_CORE
1148 
1149         MULADDC_CORE   MULADDC_CORE
1150         MULADDC_CORE   MULADDC_CORE
1151         MULADDC_STOP
1152     }
1153 
1154     for( ; i > 0; i-- )
1155     {
1156         MULADDC_INIT
1157         MULADDC_CORE
1158         MULADDC_STOP
1159     }
1160 #endif /* MULADDC_HUIT */
1161 
1162     t++;
1163 
1164     do {
1165         *d += c; c = ( *d < c ); d++;
1166     }
1167     while( c != 0 );
1168 }
1169 
1170 /*
1171  * Baseline multiplication: X = A * B  (HAC 14.12)
1172  */
1173 int mbedtls_mpi_mul_mpi( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *B )
1174 {
1175     int ret;
1176     size_t i, j;
1177     mbedtls_mpi TA, TB;
1178 
1179     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
1180 
1181     if( X == A ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) ); A = &TA; }
1182     if( X == B ) { MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) ); B = &TB; }
1183 
1184     for( i = A->n; i > 0; i-- )
1185         if( A->p[i - 1] != 0 )
1186             break;
1187 
1188     for( j = B->n; j > 0; j-- )
1189         if( B->p[j - 1] != 0 )
1190             break;
1191 
1192     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, i + j ) );
1193     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( X, 0 ) );
1194 
1195     for( i++; j > 0; j-- )
1196         mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1197 
1198     X->s = A->s * B->s;
1199 
1200 cleanup:
1201 
1202     mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TA );
1203 
1204     return( ret );
1205 }
1206 
1207 /*
1208  * Baseline multiplication: X = A * b
1209  */
1210 int mbedtls_mpi_mul_int( mbedtls_mpi *X, const mbedtls_mpi *A, mbedtls_mpi_uint b )
1211 {
1212     mbedtls_mpi _B;
1213     mbedtls_mpi_uint p[1];
1214 
1215     _B.s = 1;
1216     _B.n = 1;
1217     _B.p = p;
1218     p[0] = b;
1219 
1220     return( mbedtls_mpi_mul_mpi( X, A, &_B ) );
1221 }
1222 
1223 /*
1224  * Unsigned integer divide - double mbedtls_mpi_uint dividend, u1/u0, and
1225  * mbedtls_mpi_uint divisor, d
1226  */
1227 static mbedtls_mpi_uint mbedtls_int_div_int( mbedtls_mpi_uint u1,
1228             mbedtls_mpi_uint u0, mbedtls_mpi_uint d, mbedtls_mpi_uint *r )
1229 {
1230 #if defined(MBEDTLS_HAVE_UDBL)
1231     mbedtls_t_udbl dividend, quotient;
1232 #else
1233     const mbedtls_mpi_uint radix = (mbedtls_mpi_uint) 1 << biH;
1234     const mbedtls_mpi_uint uint_halfword_mask = ( (mbedtls_mpi_uint) 1 << biH ) - 1;
1235     mbedtls_mpi_uint d0, d1, q0, q1, rAX, r0, quotient;
1236     mbedtls_mpi_uint u0_msw, u0_lsw;
1237     size_t s;
1238 #endif
1239 
1240     /*
1241      * Check for overflow
1242      */
1243     if( 0 == d || u1 >= d )
1244     {
1245         if (r != NULL) *r = ~0;
1246 
1247         return ( ~0 );
1248     }
1249 
1250 #if defined(MBEDTLS_HAVE_UDBL)
1251     dividend  = (mbedtls_t_udbl) u1 << biL;
1252     dividend |= (mbedtls_t_udbl) u0;
1253     quotient = dividend / d;
1254     if( quotient > ( (mbedtls_t_udbl) 1 << biL ) - 1 )
1255         quotient = ( (mbedtls_t_udbl) 1 << biL ) - 1;
1256 
1257     if( r != NULL )
1258         *r = (mbedtls_mpi_uint)( dividend - (quotient * d ) );
1259 
1260     return (mbedtls_mpi_uint) quotient;
1261 #else
1262 
1263     /*
1264      * Algorithm D, Section 4.3.1 - The Art of Computer Programming
1265      *   Vol. 2 - Seminumerical Algorithms, Knuth
1266      */
1267 
1268     /*
1269      * Normalize the divisor, d, and dividend, u0, u1
1270      */
1271     s = mbedtls_clz( d );
1272     d = d << s;
1273 
1274     u1 = u1 << s;
1275     u1 |= ( u0 >> ( biL - s ) ) & ( -(mbedtls_mpi_sint)s >> ( biL - 1 ) );
1276     u0 =  u0 << s;
1277 
1278     d1 = d >> biH;
1279     d0 = d & uint_halfword_mask;
1280 
1281     u0_msw = u0 >> biH;
1282     u0_lsw = u0 & uint_halfword_mask;
1283 
1284     /*
1285      * Find the first quotient and remainder
1286      */
1287     q1 = u1 / d1;
1288     r0 = u1 - d1 * q1;
1289 
1290     while( q1 >= radix || ( q1 * d0 > radix * r0 + u0_msw ) )
1291     {
1292         q1 -= 1;
1293         r0 += d1;
1294 
1295         if ( r0 >= radix ) break;
1296     }
1297 
1298     rAX = ( u1 * radix ) + ( u0_msw - q1 * d );
1299     q0 = rAX / d1;
1300     r0 = rAX - q0 * d1;
1301 
1302     while( q0 >= radix || ( q0 * d0 > radix * r0 + u0_lsw ) )
1303     {
1304         q0 -= 1;
1305         r0 += d1;
1306 
1307         if ( r0 >= radix ) break;
1308     }
1309 
1310     if (r != NULL)
1311         *r = ( rAX * radix + u0_lsw - q0 * d ) >> s;
1312 
1313     quotient = q1 * radix + q0;
1314 
1315     return quotient;
1316 #endif
1317 }
1318 
1319 /*
1320  * Division by mbedtls_mpi: A = Q * B + R  (HAC 14.20)
1321  */
1322 int mbedtls_mpi_div_mpi( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1323 {
1324     int ret;
1325     size_t i, n, t, k;
1326     mbedtls_mpi X, Y, Z, T1, T2;
1327 
1328     if( mbedtls_mpi_cmp_int( B, 0 ) == 0 )
1329         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1330 
1331     mbedtls_mpi_init( &X ); mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &Z );
1332     mbedtls_mpi_init( &T1 ); mbedtls_mpi_init( &T2 );
1333 
1334     if( mbedtls_mpi_cmp_abs( A, B ) < 0 )
1335     {
1336         if( Q != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_lset( Q, 0 ) );
1337         if( R != NULL ) MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, A ) );
1338         return( 0 );
1339     }
1340 
1341     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &X, A ) );
1342     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, B ) );
1343     X.s = Y.s = 1;
1344 
1345     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &Z, A->n + 2 ) );
1346     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Z,  0 ) );
1347     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T1, 2 ) );
1348     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T2, 3 ) );
1349 
1350     k = mbedtls_mpi_bitlen( &Y ) % biL;
1351     if( k < biL - 1 )
1352     {
1353         k = biL - 1 - k;
1354         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &X, k ) );
1355         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, k ) );
1356     }
1357     else k = 0;
1358 
1359     n = X.n - 1;
1360     t = Y.n - 1;
1361     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &Y, biL * ( n - t ) ) );
1362 
1363     while( mbedtls_mpi_cmp_mpi( &X, &Y ) >= 0 )
1364     {
1365         Z.p[n - t]++;
1366         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &Y ) );
1367     }
1368     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, biL * ( n - t ) ) );
1369 
1370     for( i = n; i > t ; i-- )
1371     {
1372         if( X.p[i] >= Y.p[t] )
1373             Z.p[i - t - 1] = ~0;
1374         else
1375         {
1376             Z.p[i - t - 1] = mbedtls_int_div_int( X.p[i], X.p[i - 1],
1377                                                             Y.p[t], NULL);
1378         }
1379 
1380         Z.p[i - t - 1]++;
1381         do
1382         {
1383             Z.p[i - t - 1]--;
1384 
1385             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T1, 0 ) );
1386             T1.p[0] = ( t < 1 ) ? 0 : Y.p[t - 1];
1387             T1.p[1] = Y.p[t];
1388             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1389 
1390             MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &T2, 0 ) );
1391             T2.p[0] = ( i < 2 ) ? 0 : X.p[i - 2];
1392             T2.p[1] = ( i < 1 ) ? 0 : X.p[i - 1];
1393             T2.p[2] = X.p[i];
1394         }
1395         while( mbedtls_mpi_cmp_mpi( &T1, &T2 ) > 0 );
1396 
1397         MBEDTLS_MPI_CHK( mbedtls_mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1398         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1,  biL * ( i - t - 1 ) ) );
1399         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &X, &X, &T1 ) );
1400 
1401         if( mbedtls_mpi_cmp_int( &X, 0 ) < 0 )
1402         {
1403             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &T1, &Y ) );
1404             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &T1, biL * ( i - t - 1 ) ) );
1405             MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &X, &X, &T1 ) );
1406             Z.p[i - t - 1]--;
1407         }
1408     }
1409 
1410     if( Q != NULL )
1411     {
1412         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( Q, &Z ) );
1413         Q->s = A->s * B->s;
1414     }
1415 
1416     if( R != NULL )
1417     {
1418         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &X, k ) );
1419         X.s = A->s;
1420         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( R, &X ) );
1421 
1422         if( mbedtls_mpi_cmp_int( R, 0 ) == 0 )
1423             R->s = 1;
1424     }
1425 
1426 cleanup:
1427 
1428     mbedtls_mpi_free( &X ); mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &Z );
1429     mbedtls_mpi_free( &T1 ); mbedtls_mpi_free( &T2 );
1430 
1431     return( ret );
1432 }
1433 
1434 /*
1435  * Division by int: A = Q * b + R
1436  */
1437 int mbedtls_mpi_div_int( mbedtls_mpi *Q, mbedtls_mpi *R, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1438 {
1439     mbedtls_mpi _B;
1440     mbedtls_mpi_uint p[1];
1441 
1442     p[0] = ( b < 0 ) ? -b : b;
1443     _B.s = ( b < 0 ) ? -1 : 1;
1444     _B.n = 1;
1445     _B.p = p;
1446 
1447     return( mbedtls_mpi_div_mpi( Q, R, A, &_B ) );
1448 }
1449 
1450 /*
1451  * Modulo: R = A mod B
1452  */
1453 int mbedtls_mpi_mod_mpi( mbedtls_mpi *R, const mbedtls_mpi *A, const mbedtls_mpi *B )
1454 {
1455     int ret;
1456 
1457     if( mbedtls_mpi_cmp_int( B, 0 ) < 0 )
1458         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1459 
1460     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( NULL, R, A, B ) );
1461 
1462     while( mbedtls_mpi_cmp_int( R, 0 ) < 0 )
1463       MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( R, R, B ) );
1464 
1465     while( mbedtls_mpi_cmp_mpi( R, B ) >= 0 )
1466       MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( R, R, B ) );
1467 
1468 cleanup:
1469 
1470     return( ret );
1471 }
1472 
1473 /*
1474  * Modulo: r = A mod b
1475  */
1476 int mbedtls_mpi_mod_int( mbedtls_mpi_uint *r, const mbedtls_mpi *A, mbedtls_mpi_sint b )
1477 {
1478     size_t i;
1479     mbedtls_mpi_uint x, y, z;
1480 
1481     if( b == 0 )
1482         return( MBEDTLS_ERR_MPI_DIVISION_BY_ZERO );
1483 
1484     if( b < 0 )
1485         return( MBEDTLS_ERR_MPI_NEGATIVE_VALUE );
1486 
1487     /*
1488      * handle trivial cases
1489      */
1490     if( b == 1 )
1491     {
1492         *r = 0;
1493         return( 0 );
1494     }
1495 
1496     if( b == 2 )
1497     {
1498         *r = A->p[0] & 1;
1499         return( 0 );
1500     }
1501 
1502     /*
1503      * general case
1504      */
1505     for( i = A->n, y = 0; i > 0; i-- )
1506     {
1507         x  = A->p[i - 1];
1508         y  = ( y << biH ) | ( x >> biH );
1509         z  = y / b;
1510         y -= z * b;
1511 
1512         x <<= biH;
1513         y  = ( y << biH ) | ( x >> biH );
1514         z  = y / b;
1515         y -= z * b;
1516     }
1517 
1518     /*
1519      * If A is negative, then the current y represents a negative value.
1520      * Flipping it to the positive side.
1521      */
1522     if( A->s < 0 && y != 0 )
1523         y = b - y;
1524 
1525     *r = y;
1526 
1527     return( 0 );
1528 }
1529 
1530 /*
1531  * Fast Montgomery initialization (thanks to Tom St Denis)
1532  */
1533 void mbedtls_mpi_montg_init( mbedtls_mpi_uint *mm, const mbedtls_mpi *N )
1534 {
1535     mbedtls_mpi_uint x, m0 = N->p[0];
1536     unsigned int i;
1537 
1538     x  = m0;
1539     x += ( ( m0 + 2 ) & 4 ) << 1;
1540 
1541     for( i = biL; i >= 8; i /= 2 )
1542         x *= ( 2 - ( m0 * x ) );
1543 
1544     *mm = ~x + 1;
1545 }
1546 
1547 /*
1548  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1549  */
1550 int mbedtls_mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B,
1551 			 const mbedtls_mpi *N, mbedtls_mpi_uint mm,
1552                          const mbedtls_mpi *T )
1553 {
1554     size_t i, n, m;
1555     mbedtls_mpi_uint u0, u1, *d;
1556 
1557     if( T->n < N->n + 1 || T->p == NULL )
1558         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1559 
1560     memset( T->p, 0, T->n * ciL );
1561 
1562     d = T->p;
1563     n = N->n;
1564     m = ( B->n < n ) ? B->n : n;
1565 
1566     for( i = 0; i < n; i++ )
1567     {
1568         /*
1569          * T = (T + u0*B + u1*N) / 2^biL
1570          */
1571         u0 = A->p[i];
1572         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1573 
1574         mpi_mul_hlp( m, B->p, d, u0 );
1575         mpi_mul_hlp( n, N->p, d, u1 );
1576 
1577         *d++ = u0; d[n + 1] = 0;
1578     }
1579 
1580     memcpy( A->p, d, ( n + 1 ) * ciL );
1581 
1582     if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )
1583         mpi_sub_hlp( n, N->p, A->p );
1584     else
1585         /* prevent timing attacks */
1586         mpi_sub_hlp( n, A->p, T->p );
1587 
1588     return( 0 );
1589 }
1590 
1591 /*
1592  * Montgomery reduction: A = A * R^-1 mod N
1593  */
1594 int mbedtls_mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N,
1595 			 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( mbedtls_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     mbedtls_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( mbedtls_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( mbedtls_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( mbedtls_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( mbedtls_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( mbedtls_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( mbedtls_mpi_montmul( X, X, N, mm, &T ) );
1765 
1766             /*
1767              * X = X * W[wbits] R^-1 mod N
1768              */
1769             MBEDTLS_MPI_CHK( mbedtls_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( mbedtls_mpi_montmul( X, X, N, mm, &T ) );
1783 
1784         wbits <<= 1;
1785 
1786         if( ( wbits & ( one << wsize ) ) != 0 )
1787             MBEDTLS_MPI_CHK( mbedtls_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( mbedtls_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