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