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