xref: /optee_os/lib/libmbedtls/mbedtls/library/bignum.c (revision 817466cb476de705a8e3dabe1ef165fe27a18c2f)
1 /*
2  *  Multi-precision integer library
3  *
4  *  Copyright (C) 2006-2015, ARM Limited, All Rights Reserved
5  *  SPDX-License-Identifier: Apache-2.0
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 static void 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 static int mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,
1551                          const mbedtls_mpi *T )
1552 {
1553     size_t i, n, m;
1554     mbedtls_mpi_uint u0, u1, *d;
1555 
1556     if( T->n < N->n + 1 || T->p == NULL )
1557         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1558 
1559     memset( T->p, 0, T->n * ciL );
1560 
1561     d = T->p;
1562     n = N->n;
1563     m = ( B->n < n ) ? B->n : n;
1564 
1565     for( i = 0; i < n; i++ )
1566     {
1567         /*
1568          * T = (T + u0*B + u1*N) / 2^biL
1569          */
1570         u0 = A->p[i];
1571         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1572 
1573         mpi_mul_hlp( m, B->p, d, u0 );
1574         mpi_mul_hlp( n, N->p, d, u1 );
1575 
1576         *d++ = u0; d[n + 1] = 0;
1577     }
1578 
1579     memcpy( A->p, d, ( n + 1 ) * ciL );
1580 
1581     if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )
1582         mpi_sub_hlp( n, N->p, A->p );
1583     else
1584         /* prevent timing attacks */
1585         mpi_sub_hlp( n, A->p, T->p );
1586 
1587     return( 0 );
1588 }
1589 
1590 /*
1591  * Montgomery reduction: A = A * R^-1 mod N
1592  */
1593 static int mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N, mbedtls_mpi_uint mm, const mbedtls_mpi *T )
1594 {
1595     mbedtls_mpi_uint z = 1;
1596     mbedtls_mpi U;
1597 
1598     U.n = U.s = (int) z;
1599     U.p = &z;
1600 
1601     return( mpi_montmul( A, &U, N, mm, T ) );
1602 }
1603 
1604 /*
1605  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1606  */
1607 int mbedtls_mpi_exp_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *E, const mbedtls_mpi *N, mbedtls_mpi *_RR )
1608 {
1609     int ret;
1610     size_t wbits, wsize, one = 1;
1611     size_t i, j, nblimbs;
1612     size_t bufsize, nbits;
1613     mbedtls_mpi_uint ei, mm, state;
1614     mbedtls_mpi RR, T, W[ 2 << MBEDTLS_MPI_WINDOW_SIZE ], Apos;
1615     int neg;
1616 
1617     if( mbedtls_mpi_cmp_int( N, 0 ) < 0 || ( N->p[0] & 1 ) == 0 )
1618         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1619 
1620     if( mbedtls_mpi_cmp_int( E, 0 ) < 0 )
1621         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1622 
1623     /*
1624      * Init temps and window size
1625      */
1626     mpi_montg_init( &mm, N );
1627     mbedtls_mpi_init( &RR ); mbedtls_mpi_init( &T );
1628     mbedtls_mpi_init( &Apos );
1629     memset( W, 0, sizeof( W ) );
1630 
1631     i = mbedtls_mpi_bitlen( E );
1632 
1633     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1634             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1635 
1636     if( wsize > MBEDTLS_MPI_WINDOW_SIZE )
1637         wsize = MBEDTLS_MPI_WINDOW_SIZE;
1638 
1639     j = N->n + 1;
1640     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( X, j ) );
1641     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[1],  j ) );
1642     MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &T, j * 2 ) );
1643 
1644     /*
1645      * Compensate for negative A (and correct at the end)
1646      */
1647     neg = ( A->s == -1 );
1648     if( neg )
1649     {
1650         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Apos, A ) );
1651         Apos.s = 1;
1652         A = &Apos;
1653     }
1654 
1655     /*
1656      * If 1st call, pre-compute R^2 mod N
1657      */
1658     if( _RR == NULL || _RR->p == NULL )
1659     {
1660         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &RR, 1 ) );
1661         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &RR, N->n * 2 * biL ) );
1662         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &RR, &RR, N ) );
1663 
1664         if( _RR != NULL )
1665             memcpy( _RR, &RR, sizeof( mbedtls_mpi ) );
1666     }
1667     else
1668         memcpy( &RR, _RR, sizeof( mbedtls_mpi ) );
1669 
1670     /*
1671      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1672      */
1673     if( mbedtls_mpi_cmp_mpi( A, N ) >= 0 )
1674         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &W[1], A, N ) );
1675     else
1676         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[1], A ) );
1677 
1678     MBEDTLS_MPI_CHK( mpi_montmul( &W[1], &RR, N, mm, &T ) );
1679 
1680     /*
1681      * X = R^2 * R^-1 mod N = R mod N
1682      */
1683     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &RR ) );
1684     MBEDTLS_MPI_CHK( mpi_montred( X, N, mm, &T ) );
1685 
1686     if( wsize > 1 )
1687     {
1688         /*
1689          * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1690          */
1691         j =  one << ( wsize - 1 );
1692 
1693         MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[j], N->n + 1 ) );
1694         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[j], &W[1]    ) );
1695 
1696         for( i = 0; i < wsize - 1; i++ )
1697             MBEDTLS_MPI_CHK( mpi_montmul( &W[j], &W[j], N, mm, &T ) );
1698 
1699         /*
1700          * W[i] = W[i - 1] * W[1]
1701          */
1702         for( i = j + 1; i < ( one << wsize ); i++ )
1703         {
1704             MBEDTLS_MPI_CHK( mbedtls_mpi_grow( &W[i], N->n + 1 ) );
1705             MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &W[i], &W[i - 1] ) );
1706 
1707             MBEDTLS_MPI_CHK( mpi_montmul( &W[i], &W[1], N, mm, &T ) );
1708         }
1709     }
1710 
1711     nblimbs = E->n;
1712     bufsize = 0;
1713     nbits   = 0;
1714     wbits   = 0;
1715     state   = 0;
1716 
1717     while( 1 )
1718     {
1719         if( bufsize == 0 )
1720         {
1721             if( nblimbs == 0 )
1722                 break;
1723 
1724             nblimbs--;
1725 
1726             bufsize = sizeof( mbedtls_mpi_uint ) << 3;
1727         }
1728 
1729         bufsize--;
1730 
1731         ei = (E->p[nblimbs] >> bufsize) & 1;
1732 
1733         /*
1734          * skip leading 0s
1735          */
1736         if( ei == 0 && state == 0 )
1737             continue;
1738 
1739         if( ei == 0 && state == 1 )
1740         {
1741             /*
1742              * out of window, square X
1743              */
1744             MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1745             continue;
1746         }
1747 
1748         /*
1749          * add ei to current window
1750          */
1751         state = 2;
1752 
1753         nbits++;
1754         wbits |= ( ei << ( wsize - nbits ) );
1755 
1756         if( nbits == wsize )
1757         {
1758             /*
1759              * X = X^wsize R^-1 mod N
1760              */
1761             for( i = 0; i < wsize; i++ )
1762                 MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1763 
1764             /*
1765              * X = X * W[wbits] R^-1 mod N
1766              */
1767             MBEDTLS_MPI_CHK( mpi_montmul( X, &W[wbits], N, mm, &T ) );
1768 
1769             state--;
1770             nbits = 0;
1771             wbits = 0;
1772         }
1773     }
1774 
1775     /*
1776      * process the remaining bits
1777      */
1778     for( i = 0; i < nbits; i++ )
1779     {
1780         MBEDTLS_MPI_CHK( mpi_montmul( X, X, N, mm, &T ) );
1781 
1782         wbits <<= 1;
1783 
1784         if( ( wbits & ( one << wsize ) ) != 0 )
1785             MBEDTLS_MPI_CHK( mpi_montmul( X, &W[1], N, mm, &T ) );
1786     }
1787 
1788     /*
1789      * X = A^E * R * R^-1 mod N = A^E mod N
1790      */
1791     MBEDTLS_MPI_CHK( mpi_montred( X, N, mm, &T ) );
1792 
1793     if( neg && E->n != 0 && ( E->p[0] & 1 ) != 0 )
1794     {
1795         X->s = -1;
1796         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( X, N, X ) );
1797     }
1798 
1799 cleanup:
1800 
1801     for( i = ( one << ( wsize - 1 ) ); i < ( one << wsize ); i++ )
1802         mbedtls_mpi_free( &W[i] );
1803 
1804     mbedtls_mpi_free( &W[1] ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &Apos );
1805 
1806     if( _RR == NULL || _RR->p == NULL )
1807         mbedtls_mpi_free( &RR );
1808 
1809     return( ret );
1810 }
1811 
1812 /*
1813  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1814  */
1815 int mbedtls_mpi_gcd( mbedtls_mpi *G, const mbedtls_mpi *A, const mbedtls_mpi *B )
1816 {
1817     int ret;
1818     size_t lz, lzt;
1819     mbedtls_mpi TG, TA, TB;
1820 
1821     mbedtls_mpi_init( &TG ); mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TB );
1822 
1823     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TA, A ) );
1824     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, B ) );
1825 
1826     lz = mbedtls_mpi_lsb( &TA );
1827     lzt = mbedtls_mpi_lsb( &TB );
1828 
1829     if( lzt < lz )
1830         lz = lzt;
1831 
1832     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, lz ) );
1833     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, lz ) );
1834 
1835     TA.s = TB.s = 1;
1836 
1837     while( mbedtls_mpi_cmp_int( &TA, 0 ) != 0 )
1838     {
1839         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, mbedtls_mpi_lsb( &TA ) ) );
1840         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, mbedtls_mpi_lsb( &TB ) ) );
1841 
1842         if( mbedtls_mpi_cmp_mpi( &TA, &TB ) >= 0 )
1843         {
1844             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TA, &TA, &TB ) );
1845             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TA, 1 ) );
1846         }
1847         else
1848         {
1849             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_abs( &TB, &TB, &TA ) );
1850             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TB, 1 ) );
1851         }
1852     }
1853 
1854     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_l( &TB, lz ) );
1855     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( G, &TB ) );
1856 
1857 cleanup:
1858 
1859     mbedtls_mpi_free( &TG ); mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TB );
1860 
1861     return( ret );
1862 }
1863 
1864 /*
1865  * Fill X with size bytes of random.
1866  *
1867  * Use a temporary bytes representation to make sure the result is the same
1868  * regardless of the platform endianness (useful when f_rng is actually
1869  * deterministic, eg for tests).
1870  */
1871 int mbedtls_mpi_fill_random( mbedtls_mpi *X, size_t size,
1872                      int (*f_rng)(void *, unsigned char *, size_t),
1873                      void *p_rng )
1874 {
1875     int ret;
1876     unsigned char buf[MBEDTLS_MPI_MAX_SIZE];
1877 
1878     if( size > MBEDTLS_MPI_MAX_SIZE )
1879         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1880 
1881     MBEDTLS_MPI_CHK( f_rng( p_rng, buf, size ) );
1882     MBEDTLS_MPI_CHK( mbedtls_mpi_read_binary( X, buf, size ) );
1883 
1884 cleanup:
1885     return( ret );
1886 }
1887 
1888 /*
1889  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1890  */
1891 int mbedtls_mpi_inv_mod( mbedtls_mpi *X, const mbedtls_mpi *A, const mbedtls_mpi *N )
1892 {
1893     int ret;
1894     mbedtls_mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1895 
1896     if( mbedtls_mpi_cmp_int( N, 1 ) <= 0 )
1897         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
1898 
1899     mbedtls_mpi_init( &TA ); mbedtls_mpi_init( &TU ); mbedtls_mpi_init( &U1 ); mbedtls_mpi_init( &U2 );
1900     mbedtls_mpi_init( &G ); mbedtls_mpi_init( &TB ); mbedtls_mpi_init( &TV );
1901     mbedtls_mpi_init( &V1 ); mbedtls_mpi_init( &V2 );
1902 
1903     MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &G, A, N ) );
1904 
1905     if( mbedtls_mpi_cmp_int( &G, 1 ) != 0 )
1906     {
1907         ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
1908         goto cleanup;
1909     }
1910 
1911     MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &TA, A, N ) );
1912     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TU, &TA ) );
1913     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TB, N ) );
1914     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &TV, N ) );
1915 
1916     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U1, 1 ) );
1917     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &U2, 0 ) );
1918     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V1, 0 ) );
1919     MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &V2, 1 ) );
1920 
1921     do
1922     {
1923         while( ( TU.p[0] & 1 ) == 0 )
1924         {
1925             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TU, 1 ) );
1926 
1927             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1928             {
1929                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &U1, &U1, &TB ) );
1930                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &TA ) );
1931             }
1932 
1933             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U1, 1 ) );
1934             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &U2, 1 ) );
1935         }
1936 
1937         while( ( TV.p[0] & 1 ) == 0 )
1938         {
1939             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &TV, 1 ) );
1940 
1941             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1942             {
1943                 MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, &TB ) );
1944                 MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &TA ) );
1945             }
1946 
1947             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V1, 1 ) );
1948             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &V2, 1 ) );
1949         }
1950 
1951         if( mbedtls_mpi_cmp_mpi( &TU, &TV ) >= 0 )
1952         {
1953             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TU, &TU, &TV ) );
1954             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U1, &U1, &V1 ) );
1955             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &U2, &U2, &V2 ) );
1956         }
1957         else
1958         {
1959             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &TV, &TV, &TU ) );
1960             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, &U1 ) );
1961             MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V2, &V2, &U2 ) );
1962         }
1963     }
1964     while( mbedtls_mpi_cmp_int( &TU, 0 ) != 0 );
1965 
1966     while( mbedtls_mpi_cmp_int( &V1, 0 ) < 0 )
1967         MBEDTLS_MPI_CHK( mbedtls_mpi_add_mpi( &V1, &V1, N ) );
1968 
1969     while( mbedtls_mpi_cmp_mpi( &V1, N ) >= 0 )
1970         MBEDTLS_MPI_CHK( mbedtls_mpi_sub_mpi( &V1, &V1, N ) );
1971 
1972     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( X, &V1 ) );
1973 
1974 cleanup:
1975 
1976     mbedtls_mpi_free( &TA ); mbedtls_mpi_free( &TU ); mbedtls_mpi_free( &U1 ); mbedtls_mpi_free( &U2 );
1977     mbedtls_mpi_free( &G ); mbedtls_mpi_free( &TB ); mbedtls_mpi_free( &TV );
1978     mbedtls_mpi_free( &V1 ); mbedtls_mpi_free( &V2 );
1979 
1980     return( ret );
1981 }
1982 
1983 #if defined(MBEDTLS_GENPRIME)
1984 
1985 static const int small_prime[] =
1986 {
1987         3,    5,    7,   11,   13,   17,   19,   23,
1988        29,   31,   37,   41,   43,   47,   53,   59,
1989        61,   67,   71,   73,   79,   83,   89,   97,
1990       101,  103,  107,  109,  113,  127,  131,  137,
1991       139,  149,  151,  157,  163,  167,  173,  179,
1992       181,  191,  193,  197,  199,  211,  223,  227,
1993       229,  233,  239,  241,  251,  257,  263,  269,
1994       271,  277,  281,  283,  293,  307,  311,  313,
1995       317,  331,  337,  347,  349,  353,  359,  367,
1996       373,  379,  383,  389,  397,  401,  409,  419,
1997       421,  431,  433,  439,  443,  449,  457,  461,
1998       463,  467,  479,  487,  491,  499,  503,  509,
1999       521,  523,  541,  547,  557,  563,  569,  571,
2000       577,  587,  593,  599,  601,  607,  613,  617,
2001       619,  631,  641,  643,  647,  653,  659,  661,
2002       673,  677,  683,  691,  701,  709,  719,  727,
2003       733,  739,  743,  751,  757,  761,  769,  773,
2004       787,  797,  809,  811,  821,  823,  827,  829,
2005       839,  853,  857,  859,  863,  877,  881,  883,
2006       887,  907,  911,  919,  929,  937,  941,  947,
2007       953,  967,  971,  977,  983,  991,  997, -103
2008 };
2009 
2010 /*
2011  * Small divisors test (X must be positive)
2012  *
2013  * Return values:
2014  * 0: no small factor (possible prime, more tests needed)
2015  * 1: certain prime
2016  * MBEDTLS_ERR_MPI_NOT_ACCEPTABLE: certain non-prime
2017  * other negative: error
2018  */
2019 static int mpi_check_small_factors( const mbedtls_mpi *X )
2020 {
2021     int ret = 0;
2022     size_t i;
2023     mbedtls_mpi_uint r;
2024 
2025     if( ( X->p[0] & 1 ) == 0 )
2026         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2027 
2028     for( i = 0; small_prime[i] > 0; i++ )
2029     {
2030         if( mbedtls_mpi_cmp_int( X, small_prime[i] ) <= 0 )
2031             return( 1 );
2032 
2033         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, small_prime[i] ) );
2034 
2035         if( r == 0 )
2036             return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2037     }
2038 
2039 cleanup:
2040     return( ret );
2041 }
2042 
2043 /*
2044  * Miller-Rabin pseudo-primality test  (HAC 4.24)
2045  */
2046 static int mpi_miller_rabin( const mbedtls_mpi *X,
2047                              int (*f_rng)(void *, unsigned char *, size_t),
2048                              void *p_rng )
2049 {
2050     int ret, count;
2051     size_t i, j, k, n, s;
2052     mbedtls_mpi W, R, T, A, RR;
2053 
2054     mbedtls_mpi_init( &W ); mbedtls_mpi_init( &R ); mbedtls_mpi_init( &T ); mbedtls_mpi_init( &A );
2055     mbedtls_mpi_init( &RR );
2056 
2057     /*
2058      * W = |X| - 1
2059      * R = W >> lsb( W )
2060      */
2061     MBEDTLS_MPI_CHK( mbedtls_mpi_sub_int( &W, X, 1 ) );
2062     s = mbedtls_mpi_lsb( &W );
2063     MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &R, &W ) );
2064     MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &R, s ) );
2065 
2066     i = mbedtls_mpi_bitlen( X );
2067     /*
2068      * HAC, table 4.4
2069      */
2070     n = ( ( i >= 1300 ) ?  2 : ( i >=  850 ) ?  3 :
2071           ( i >=  650 ) ?  4 : ( i >=  350 ) ?  8 :
2072           ( i >=  250 ) ? 12 : ( i >=  150 ) ? 18 : 27 );
2073 
2074     for( i = 0; i < n; i++ )
2075     {
2076         /*
2077          * pick a random A, 1 < A < |X| - 1
2078          */
2079         MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
2080 
2081         if( mbedtls_mpi_cmp_mpi( &A, &W ) >= 0 )
2082         {
2083             j = mbedtls_mpi_bitlen( &A ) - mbedtls_mpi_bitlen( &W );
2084             MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &A, j + 1 ) );
2085         }
2086         A.p[0] |= 3;
2087 
2088         count = 0;
2089         do {
2090             MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
2091 
2092             j = mbedtls_mpi_bitlen( &A );
2093             k = mbedtls_mpi_bitlen( &W );
2094             if (j > k) {
2095                 MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &A, j - k ) );
2096             }
2097 
2098             if (count++ > 30) {
2099                 return MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2100             }
2101 
2102         } while ( mbedtls_mpi_cmp_mpi( &A, &W ) >= 0 ||
2103                   mbedtls_mpi_cmp_int( &A, 1 )  <= 0    );
2104 
2105         /*
2106          * A = A^R mod |X|
2107          */
2108         MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &A, &A, &R, X, &RR ) );
2109 
2110         if( mbedtls_mpi_cmp_mpi( &A, &W ) == 0 ||
2111             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2112             continue;
2113 
2114         j = 1;
2115         while( j < s && mbedtls_mpi_cmp_mpi( &A, &W ) != 0 )
2116         {
2117             /*
2118              * A = A * A mod |X|
2119              */
2120             MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &T, &A, &A ) );
2121             MBEDTLS_MPI_CHK( mbedtls_mpi_mod_mpi( &A, &T, X  ) );
2122 
2123             if( mbedtls_mpi_cmp_int( &A, 1 ) == 0 )
2124                 break;
2125 
2126             j++;
2127         }
2128 
2129         /*
2130          * not prime if A != |X| - 1 or A == 1
2131          */
2132         if( mbedtls_mpi_cmp_mpi( &A, &W ) != 0 ||
2133             mbedtls_mpi_cmp_int( &A,  1 ) == 0 )
2134         {
2135             ret = MBEDTLS_ERR_MPI_NOT_ACCEPTABLE;
2136             break;
2137         }
2138     }
2139 
2140 cleanup:
2141     mbedtls_mpi_free( &W ); mbedtls_mpi_free( &R ); mbedtls_mpi_free( &T ); mbedtls_mpi_free( &A );
2142     mbedtls_mpi_free( &RR );
2143 
2144     return( ret );
2145 }
2146 
2147 /*
2148  * Pseudo-primality test: small factors, then Miller-Rabin
2149  */
2150 int mbedtls_mpi_is_prime( const mbedtls_mpi *X,
2151                   int (*f_rng)(void *, unsigned char *, size_t),
2152                   void *p_rng )
2153 {
2154     int ret;
2155     mbedtls_mpi XX;
2156 
2157     XX.s = 1;
2158     XX.n = X->n;
2159     XX.p = X->p;
2160 
2161     if( mbedtls_mpi_cmp_int( &XX, 0 ) == 0 ||
2162         mbedtls_mpi_cmp_int( &XX, 1 ) == 0 )
2163         return( MBEDTLS_ERR_MPI_NOT_ACCEPTABLE );
2164 
2165     if( mbedtls_mpi_cmp_int( &XX, 2 ) == 0 )
2166         return( 0 );
2167 
2168     if( ( ret = mpi_check_small_factors( &XX ) ) != 0 )
2169     {
2170         if( ret == 1 )
2171             return( 0 );
2172 
2173         return( ret );
2174     }
2175 
2176     return( mpi_miller_rabin( &XX, f_rng, p_rng ) );
2177 }
2178 
2179 /*
2180  * Prime number generation
2181  */
2182 int mbedtls_mpi_gen_prime( mbedtls_mpi *X, size_t nbits, int dh_flag,
2183                    int (*f_rng)(void *, unsigned char *, size_t),
2184                    void *p_rng )
2185 {
2186     int ret;
2187     size_t k, n;
2188     mbedtls_mpi_uint r;
2189     mbedtls_mpi Y;
2190 
2191     if( nbits < 3 || nbits > MBEDTLS_MPI_MAX_BITS )
2192         return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
2193 
2194     mbedtls_mpi_init( &Y );
2195 
2196     n = BITS_TO_LIMBS( nbits );
2197 
2198     MBEDTLS_MPI_CHK( mbedtls_mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
2199 
2200     k = mbedtls_mpi_bitlen( X );
2201     if( k > nbits ) MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( X, k - nbits + 1 ) );
2202 
2203     mbedtls_mpi_set_bit( X, nbits-1, 1 );
2204 
2205     X->p[0] |= 1;
2206 
2207     if( dh_flag == 0 )
2208     {
2209         while( ( ret = mbedtls_mpi_is_prime( X, f_rng, p_rng ) ) != 0 )
2210         {
2211             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2212                 goto cleanup;
2213 
2214             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 2 ) );
2215         }
2216     }
2217     else
2218     {
2219         /*
2220          * An necessary condition for Y and X = 2Y + 1 to be prime
2221          * is X = 2 mod 3 (which is equivalent to Y = 2 mod 3).
2222          * Make sure it is satisfied, while keeping X = 3 mod 4
2223          */
2224 
2225         X->p[0] |= 2;
2226 
2227         MBEDTLS_MPI_CHK( mbedtls_mpi_mod_int( &r, X, 3 ) );
2228         if( r == 0 )
2229             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 8 ) );
2230         else if( r == 1 )
2231             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( X, X, 4 ) );
2232 
2233         /* Set Y = (X-1) / 2, which is X / 2 because X is odd */
2234         MBEDTLS_MPI_CHK( mbedtls_mpi_copy( &Y, X ) );
2235         MBEDTLS_MPI_CHK( mbedtls_mpi_shift_r( &Y, 1 ) );
2236 
2237         while( 1 )
2238         {
2239             /*
2240              * First, check small factors for X and Y
2241              * before doing Miller-Rabin on any of them
2242              */
2243             if( ( ret = mpi_check_small_factors(  X         ) ) == 0 &&
2244                 ( ret = mpi_check_small_factors( &Y         ) ) == 0 &&
2245                 ( ret = mpi_miller_rabin(  X, f_rng, p_rng  ) ) == 0 &&
2246                 ( ret = mpi_miller_rabin( &Y, f_rng, p_rng  ) ) == 0 )
2247             {
2248                 break;
2249             }
2250 
2251             if( ret != MBEDTLS_ERR_MPI_NOT_ACCEPTABLE )
2252                 goto cleanup;
2253 
2254             /*
2255              * Next candidates. We want to preserve Y = (X-1) / 2 and
2256              * Y = 1 mod 2 and Y = 2 mod 3 (eq X = 3 mod 4 and X = 2 mod 3)
2257              * so up Y by 6 and X by 12.
2258              */
2259             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int(  X,  X, 12 ) );
2260             MBEDTLS_MPI_CHK( mbedtls_mpi_add_int( &Y, &Y, 6  ) );
2261         }
2262     }
2263 
2264 cleanup:
2265 
2266     mbedtls_mpi_free( &Y );
2267 
2268     return( ret );
2269 }
2270 
2271 #endif /* MBEDTLS_GENPRIME */
2272 
2273 #if defined(MBEDTLS_SELF_TEST)
2274 
2275 #define GCD_PAIR_COUNT  3
2276 
2277 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
2278 {
2279     { 693, 609, 21 },
2280     { 1764, 868, 28 },
2281     { 768454923, 542167814, 1 }
2282 };
2283 
2284 /*
2285  * Checkup routine
2286  */
2287 int mbedtls_mpi_self_test( int verbose )
2288 {
2289     int ret, i;
2290     mbedtls_mpi A, E, N, X, Y, U, V;
2291 
2292     mbedtls_mpi_init( &A ); mbedtls_mpi_init( &E ); mbedtls_mpi_init( &N ); mbedtls_mpi_init( &X );
2293     mbedtls_mpi_init( &Y ); mbedtls_mpi_init( &U ); mbedtls_mpi_init( &V );
2294 
2295     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &A, 16,
2296         "EFE021C2645FD1DC586E69184AF4A31E" \
2297         "D5F53E93B5F123FA41680867BA110131" \
2298         "944FE7952E2517337780CB0DB80E61AA" \
2299         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
2300 
2301     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &E, 16,
2302         "B2E7EFD37075B9F03FF989C7C5051C20" \
2303         "34D2A323810251127E7BF8625A4F49A5" \
2304         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2305         "5B5C25763222FEFCCFC38B832366C29E" ) );
2306 
2307     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &N, 16,
2308         "0066A198186C18C10B2F5ED9B522752A" \
2309         "9830B69916E535C8F047518A889A43A5" \
2310         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2311 
2312     MBEDTLS_MPI_CHK( mbedtls_mpi_mul_mpi( &X, &A, &N ) );
2313 
2314     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2315         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2316         "9E857EA95A03512E2BAE7391688D264A" \
2317         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2318         "8001B72E848A38CAE1C65F78E56ABDEF" \
2319         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2320         "ECF677152EF804370C1A305CAF3B5BF1" \
2321         "30879B56C61DE584A0F53A2447A51E" ) );
2322 
2323     if( verbose != 0 )
2324         mbedtls_printf( "  MPI test #1 (mul_mpi): " );
2325 
2326     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2327     {
2328         if( verbose != 0 )
2329             mbedtls_printf( "failed\n" );
2330 
2331         ret = 1;
2332         goto cleanup;
2333     }
2334 
2335     if( verbose != 0 )
2336         mbedtls_printf( "passed\n" );
2337 
2338     MBEDTLS_MPI_CHK( mbedtls_mpi_div_mpi( &X, &Y, &A, &N ) );
2339 
2340     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2341         "256567336059E52CAE22925474705F39A94" ) );
2342 
2343     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &V, 16,
2344         "6613F26162223DF488E9CD48CC132C7A" \
2345         "0AC93C701B001B092E4E5B9F73BCD27B" \
2346         "9EE50D0657C77F374E903CDFA4C642" ) );
2347 
2348     if( verbose != 0 )
2349         mbedtls_printf( "  MPI test #2 (div_mpi): " );
2350 
2351     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 ||
2352         mbedtls_mpi_cmp_mpi( &Y, &V ) != 0 )
2353     {
2354         if( verbose != 0 )
2355             mbedtls_printf( "failed\n" );
2356 
2357         ret = 1;
2358         goto cleanup;
2359     }
2360 
2361     if( verbose != 0 )
2362         mbedtls_printf( "passed\n" );
2363 
2364     MBEDTLS_MPI_CHK( mbedtls_mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2365 
2366     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2367         "36E139AEA55215609D2816998ED020BB" \
2368         "BD96C37890F65171D948E9BC7CBAA4D9" \
2369         "325D24D6A3C12710F10A09FA08AB87" ) );
2370 
2371     if( verbose != 0 )
2372         mbedtls_printf( "  MPI test #3 (exp_mod): " );
2373 
2374     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2375     {
2376         if( verbose != 0 )
2377             mbedtls_printf( "failed\n" );
2378 
2379         ret = 1;
2380         goto cleanup;
2381     }
2382 
2383     if( verbose != 0 )
2384         mbedtls_printf( "passed\n" );
2385 
2386     MBEDTLS_MPI_CHK( mbedtls_mpi_inv_mod( &X, &A, &N ) );
2387 
2388     MBEDTLS_MPI_CHK( mbedtls_mpi_read_string( &U, 16,
2389         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2390         "C3DBA76456363A10869622EAC2DD84EC" \
2391         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2392 
2393     if( verbose != 0 )
2394         mbedtls_printf( "  MPI test #4 (inv_mod): " );
2395 
2396     if( mbedtls_mpi_cmp_mpi( &X, &U ) != 0 )
2397     {
2398         if( verbose != 0 )
2399             mbedtls_printf( "failed\n" );
2400 
2401         ret = 1;
2402         goto cleanup;
2403     }
2404 
2405     if( verbose != 0 )
2406         mbedtls_printf( "passed\n" );
2407 
2408     if( verbose != 0 )
2409         mbedtls_printf( "  MPI test #5 (simple gcd): " );
2410 
2411     for( i = 0; i < GCD_PAIR_COUNT; i++ )
2412     {
2413         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &X, gcd_pairs[i][0] ) );
2414         MBEDTLS_MPI_CHK( mbedtls_mpi_lset( &Y, gcd_pairs[i][1] ) );
2415 
2416         MBEDTLS_MPI_CHK( mbedtls_mpi_gcd( &A, &X, &Y ) );
2417 
2418         if( mbedtls_mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2419         {
2420             if( verbose != 0 )
2421                 mbedtls_printf( "failed at %d\n", i );
2422 
2423             ret = 1;
2424             goto cleanup;
2425         }
2426     }
2427 
2428     if( verbose != 0 )
2429         mbedtls_printf( "passed\n" );
2430 
2431 cleanup:
2432 
2433     if( ret != 0 && verbose != 0 )
2434         mbedtls_printf( "Unexpected error, return code = %08X\n", ret );
2435 
2436     mbedtls_mpi_free( &A ); mbedtls_mpi_free( &E ); mbedtls_mpi_free( &N ); mbedtls_mpi_free( &X );
2437     mbedtls_mpi_free( &Y ); mbedtls_mpi_free( &U ); mbedtls_mpi_free( &V );
2438 
2439     if( verbose != 0 )
2440         mbedtls_printf( "\n" );
2441 
2442     return( ret );
2443 }
2444 
2445 #endif /* MBEDTLS_SELF_TEST */
2446 
2447 #endif /* MBEDTLS_BIGNUM_C */
2448