import PolarSSL 1.2.6's bignum
[gnuk/gnuk.git] / polarssl / library / bignum.c
1 /*
2  *  Multi-precision integer library
3  *
4  *  Copyright (C) 2006-2010, Brainspark B.V.
5  *
6  *  This file is part of PolarSSL (http://www.polarssl.org)
7  *  Lead Maintainer: Paul Bakker <polarssl_maintainer at polarssl.org>
8  *
9  *  All rights reserved.
10  *
11  *  This program is free software; you can redistribute it and/or modify
12  *  it under the terms of the GNU General Public License as published by
13  *  the Free Software Foundation; either version 2 of the License, or
14  *  (at your option) any later version.
15  *
16  *  This program is distributed in the hope that it will be useful,
17  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19  *  GNU General Public License for more details.
20  *
21  *  You should have received a copy of the GNU General Public License along
22  *  with this program; if not, write to the Free Software Foundation, Inc.,
23  *  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24  */
25 /*
26  *  This MPI implementation is based on:
27  *
28  *  http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
29  *  http://www.stillhq.com/extracted/gnupg-api/mpi/
30  *  http://math.libtomcrypt.com/files/tommath.pdf
31  */
32
33 #include "polarssl/config.h"
34
35 #if defined(POLARSSL_BIGNUM_C)
36
37 #include "polarssl/bignum.h"
38 #include "polarssl/bn_mul.h"
39
40 #include <stdlib.h>
41
42 #define ciL    (sizeof(t_uint))         /* chars in limb  */
43 #define biL    (ciL << 3)               /* bits  in limb  */
44 #define biH    (ciL << 2)               /* half limb size */
45
46 /*
47  * Convert between bits/chars and number of limbs
48  */
49 #define BITS_TO_LIMBS(i)  (((i) + biL - 1) / biL)
50 #define CHARS_TO_LIMBS(i) (((i) + ciL - 1) / ciL)
51
52 /*
53  * Initialize one MPI
54  */
55 void mpi_init( mpi *X )
56 {
57     if( X == NULL )
58         return;
59
60     X->s = 1;
61     X->n = 0;
62     X->p = NULL;
63 }
64
65 /*
66  * Unallocate one MPI
67  */
68 void mpi_free( mpi *X )
69 {
70     if( X == NULL )
71         return;
72
73     if( X->p != NULL )
74     {
75         memset( X->p, 0, X->n * ciL );
76         free( X->p );
77     }
78
79     X->s = 1;
80     X->n = 0;
81     X->p = NULL;
82 }
83
84 /*
85  * Enlarge to the specified number of limbs
86  */
87 int mpi_grow( mpi *X, size_t nblimbs )
88 {
89     t_uint *p;
90
91     if( nblimbs > POLARSSL_MPI_MAX_LIMBS )
92         return( POLARSSL_ERR_MPI_MALLOC_FAILED );
93
94     if( X->n < nblimbs )
95     {
96         if( ( p = (t_uint *) malloc( nblimbs * ciL ) ) == NULL )
97             return( POLARSSL_ERR_MPI_MALLOC_FAILED );
98
99         memset( p, 0, nblimbs * ciL );
100
101         if( X->p != NULL )
102         {
103             memcpy( p, X->p, X->n * ciL );
104             memset( X->p, 0, X->n * ciL );
105             free( X->p );
106         }
107
108         X->n = nblimbs;
109         X->p = p;
110     }
111
112     return( 0 );
113 }
114
115 /*
116  * Copy the contents of Y into X
117  */
118 int mpi_copy( mpi *X, const mpi *Y )
119 {
120     int ret;
121     size_t i;
122
123     if( X == Y )
124         return( 0 );
125
126     for( i = Y->n - 1; i > 0; i-- )
127         if( Y->p[i] != 0 )
128             break;
129     i++;
130
131     X->s = Y->s;
132
133     MPI_CHK( mpi_grow( X, i ) );
134
135     memset( X->p, 0, X->n * ciL );
136     memcpy( X->p, Y->p, i * ciL );
137
138 cleanup:
139
140     return( ret );
141 }
142
143 /*
144  * Swap the contents of X and Y
145  */
146 void mpi_swap( mpi *X, mpi *Y )
147 {
148     mpi T;
149
150     memcpy( &T,  X, sizeof( mpi ) );
151     memcpy(  X,  Y, sizeof( mpi ) );
152     memcpy(  Y, &T, sizeof( mpi ) );
153 }
154
155 /*
156  * Set value from integer
157  */
158 int mpi_lset( mpi *X, t_sint z )
159 {
160     int ret;
161
162     MPI_CHK( mpi_grow( X, 1 ) );
163     memset( X->p, 0, X->n * ciL );
164
165     X->p[0] = ( z < 0 ) ? -z : z;
166     X->s    = ( z < 0 ) ? -1 : 1;
167
168 cleanup:
169
170     return( ret );
171 }
172
173 /*
174  * Get a specific bit
175  */
176 int mpi_get_bit( const mpi *X, size_t pos )
177 {
178     if( X->n * biL <= pos )
179         return( 0 );
180
181     return ( X->p[pos / biL] >> ( pos % biL ) ) & 0x01;
182 }
183
184 /*
185  * Set a bit to a specific value of 0 or 1
186  */
187 int mpi_set_bit( mpi *X, size_t pos, unsigned char val )
188 {
189     int ret = 0;
190     size_t off = pos / biL;
191     size_t idx = pos % biL;
192
193     if( val != 0 && val != 1 )
194         return POLARSSL_ERR_MPI_BAD_INPUT_DATA;
195         
196     if( X->n * biL <= pos )
197     {
198         if( val == 0 )
199             return ( 0 );
200
201         MPI_CHK( mpi_grow( X, off + 1 ) );
202     }
203
204     X->p[off] = ( X->p[off] & ~( 0x01 << idx ) ) | ( val << idx );
205
206 cleanup:
207     
208     return( ret );
209 }
210
211 /*
212  * Return the number of least significant bits
213  */
214 size_t mpi_lsb( const mpi *X )
215 {
216     size_t i, j, count = 0;
217
218     for( i = 0; i < X->n; i++ )
219         for( j = 0; j < biL; j++, count++ )
220             if( ( ( X->p[i] >> j ) & 1 ) != 0 )
221                 return( count );
222
223     return( 0 );
224 }
225
226 /*
227  * Return the number of most significant bits
228  */
229 size_t mpi_msb( const mpi *X )
230 {
231     size_t i, j;
232
233     for( i = X->n - 1; i > 0; i-- )
234         if( X->p[i] != 0 )
235             break;
236
237     for( j = biL; j > 0; j-- )
238         if( ( ( X->p[i] >> ( j - 1 ) ) & 1 ) != 0 )
239             break;
240
241     return( ( i * biL ) + j );
242 }
243
244 /*
245  * Return the total size in bytes
246  */
247 size_t mpi_size( const mpi *X )
248 {
249     return( ( mpi_msb( X ) + 7 ) >> 3 );
250 }
251
252 /*
253  * Convert an ASCII character to digit value
254  */
255 static int mpi_get_digit( t_uint *d, int radix, char c )
256 {
257     *d = 255;
258
259     if( c >= 0x30 && c <= 0x39 ) *d = c - 0x30;
260     if( c >= 0x41 && c <= 0x46 ) *d = c - 0x37;
261     if( c >= 0x61 && c <= 0x66 ) *d = c - 0x57;
262
263     if( *d >= (t_uint) radix )
264         return( POLARSSL_ERR_MPI_INVALID_CHARACTER );
265
266     return( 0 );
267 }
268
269 /*
270  * Import from an ASCII string
271  */
272 int mpi_read_string( mpi *X, int radix, const char *s )
273 {
274     int ret;
275     size_t i, j, slen, n;
276     t_uint d;
277     mpi T;
278
279     if( radix < 2 || radix > 16 )
280         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
281
282     mpi_init( &T );
283
284     slen = strlen( s );
285
286     if( radix == 16 )
287     {
288         n = BITS_TO_LIMBS( slen << 2 );
289
290         MPI_CHK( mpi_grow( X, n ) );
291         MPI_CHK( mpi_lset( X, 0 ) );
292
293         for( i = slen, j = 0; i > 0; i--, j++ )
294         {
295             if( i == 1 && s[i - 1] == '-' )
296             {
297                 X->s = -1;
298                 break;
299             }
300
301             MPI_CHK( mpi_get_digit( &d, radix, s[i - 1] ) );
302             X->p[j / (2 * ciL)] |= d << ( (j % (2 * ciL)) << 2 );
303         }
304     }
305     else
306     {
307         MPI_CHK( mpi_lset( X, 0 ) );
308
309         for( i = 0; i < slen; i++ )
310         {
311             if( i == 0 && s[i] == '-' )
312             {
313                 X->s = -1;
314                 continue;
315             }
316
317             MPI_CHK( mpi_get_digit( &d, radix, s[i] ) );
318             MPI_CHK( mpi_mul_int( &T, X, radix ) );
319
320             if( X->s == 1 )
321             {
322                 MPI_CHK( mpi_add_int( X, &T, d ) );
323             }
324             else
325             {
326                 MPI_CHK( mpi_sub_int( X, &T, d ) );
327             }
328         }
329     }
330
331 cleanup:
332
333     mpi_free( &T );
334
335     return( ret );
336 }
337
338 /*
339  * Helper to write the digits high-order first
340  */
341 static int mpi_write_hlp( mpi *X, int radix, char **p )
342 {
343     int ret;
344     t_uint r;
345
346     if( radix < 2 || radix > 16 )
347         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
348
349     MPI_CHK( mpi_mod_int( &r, X, radix ) );
350     MPI_CHK( mpi_div_int( X, NULL, X, radix ) );
351
352     if( mpi_cmp_int( X, 0 ) != 0 )
353         MPI_CHK( mpi_write_hlp( X, radix, p ) );
354
355     if( r < 10 )
356         *(*p)++ = (char)( r + 0x30 );
357     else
358         *(*p)++ = (char)( r + 0x37 );
359
360 cleanup:
361
362     return( ret );
363 }
364
365 /*
366  * Export into an ASCII string
367  */
368 int mpi_write_string( const mpi *X, int radix, char *s, size_t *slen )
369 {
370     int ret = 0;
371     size_t n;
372     char *p;
373     mpi T;
374
375     if( radix < 2 || radix > 16 )
376         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
377
378     n = mpi_msb( X );
379     if( radix >=  4 ) n >>= 1;
380     if( radix >= 16 ) n >>= 1;
381     n += 3;
382
383     if( *slen < n )
384     {
385         *slen = n;
386         return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );
387     }
388
389     p = s;
390     mpi_init( &T );
391
392     if( X->s == -1 )
393         *p++ = '-';
394
395     if( radix == 16 )
396     {
397         int c;
398         size_t i, j, k;
399
400         for( i = X->n, k = 0; i > 0; i-- )
401         {
402             for( j = ciL; j > 0; j-- )
403             {
404                 c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
405
406                 if( c == 0 && k == 0 && ( i + j + 3 ) != 0 )
407                     continue;
408
409                 *(p++) = "0123456789ABCDEF" [c / 16];
410                 *(p++) = "0123456789ABCDEF" [c % 16];
411                 k = 1;
412             }
413         }
414     }
415     else
416     {
417         MPI_CHK( mpi_copy( &T, X ) );
418
419         if( T.s == -1 )
420             T.s = 1;
421
422         MPI_CHK( mpi_write_hlp( &T, radix, &p ) );
423     }
424
425     *p++ = '\0';
426     *slen = p - s;
427
428 cleanup:
429
430     mpi_free( &T );
431
432     return( ret );
433 }
434
435 #if defined(POLARSSL_FS_IO)
436 /*
437  * Read X from an opened file
438  */
439 int mpi_read_file( mpi *X, int radix, FILE *fin )
440 {
441     t_uint d;
442     size_t slen;
443     char *p;
444     /*
445      * Buffer should have space for (short) label and decimal formatted MPI,
446      * newline characters and '\0'
447      */
448     char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
449
450     memset( s, 0, sizeof( s ) );
451     if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
452         return( POLARSSL_ERR_MPI_FILE_IO_ERROR );
453
454     slen = strlen( s );
455     if( slen == sizeof( s ) - 2 )
456         return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );
457
458     if( s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
459     if( s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }
460
461     p = s + slen;
462     while( --p >= s )
463         if( mpi_get_digit( &d, radix, *p ) != 0 )
464             break;
465
466     return( mpi_read_string( X, radix, p + 1 ) );
467 }
468
469 /*
470  * Write X into an opened file (or stdout if fout == NULL)
471  */
472 int mpi_write_file( const char *p, const mpi *X, int radix, FILE *fout )
473 {
474     int ret;
475     size_t n, slen, plen;
476     /*
477      * Buffer should have space for (short) label and decimal formatted MPI,
478      * newline characters and '\0'
479      */
480     char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
481
482     n = sizeof( s );
483     memset( s, 0, n );
484     n -= 2;
485
486     MPI_CHK( mpi_write_string( X, radix, s, (size_t *) &n ) );
487
488     if( p == NULL ) p = "";
489
490     plen = strlen( p );
491     slen = strlen( s );
492     s[slen++] = '\r';
493     s[slen++] = '\n';
494
495     if( fout != NULL )
496     {
497         if( fwrite( p, 1, plen, fout ) != plen ||
498             fwrite( s, 1, slen, fout ) != slen )
499             return( POLARSSL_ERR_MPI_FILE_IO_ERROR );
500     }
501     else
502         printf( "%s%s", p, s );
503
504 cleanup:
505
506     return( ret );
507 }
508 #endif /* POLARSSL_FS_IO */
509
510 /*
511  * Import X from unsigned binary data, big endian
512  */
513 int mpi_read_binary( mpi *X, const unsigned char *buf, size_t buflen )
514 {
515     int ret;
516     size_t i, j, n;
517
518     for( n = 0; n < buflen; n++ )
519         if( buf[n] != 0 )
520             break;
521
522     MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( buflen - n ) ) );
523     MPI_CHK( mpi_lset( X, 0 ) );
524
525     for( i = buflen, j = 0; i > n; i--, j++ )
526         X->p[j / ciL] |= ((t_uint) buf[i - 1]) << ((j % ciL) << 3);
527
528 cleanup:
529
530     return( ret );
531 }
532
533 /*
534  * Export X into unsigned binary data, big endian
535  */
536 int mpi_write_binary( const mpi *X, unsigned char *buf, size_t buflen )
537 {
538     size_t i, j, n;
539
540     n = mpi_size( X );
541
542     if( buflen < n )
543         return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );
544
545     memset( buf, 0, buflen );
546
547     for( i = buflen - 1, j = 0; n > 0; i--, j++, n-- )
548         buf[i] = (unsigned char)( X->p[j / ciL] >> ((j % ciL) << 3) );
549
550     return( 0 );
551 }
552
553 /*
554  * Left-shift: X <<= count
555  */
556 int mpi_shift_l( mpi *X, size_t count )
557 {
558     int ret;
559     size_t i, v0, t1;
560     t_uint r0 = 0, r1;
561
562     v0 = count / (biL    );
563     t1 = count & (biL - 1);
564
565     i = mpi_msb( X ) + count;
566
567     if( X->n * biL < i )
568         MPI_CHK( mpi_grow( X, BITS_TO_LIMBS( i ) ) );
569
570     ret = 0;
571
572     /*
573      * shift by count / limb_size
574      */
575     if( v0 > 0 )
576     {
577         for( i = X->n; i > v0; i-- )
578             X->p[i - 1] = X->p[i - v0 - 1];
579
580         for( ; i > 0; i-- )
581             X->p[i - 1] = 0;
582     }
583
584     /*
585      * shift by count % limb_size
586      */
587     if( t1 > 0 )
588     {
589         for( i = v0; i < X->n; i++ )
590         {
591             r1 = X->p[i] >> (biL - t1);
592             X->p[i] <<= t1;
593             X->p[i] |= r0;
594             r0 = r1;
595         }
596     }
597
598 cleanup:
599
600     return( ret );
601 }
602
603 /*
604  * Right-shift: X >>= count
605  */
606 int mpi_shift_r( mpi *X, size_t count )
607 {
608     size_t i, v0, v1;
609     t_uint r0 = 0, r1;
610
611     v0 = count /  biL;
612     v1 = count & (biL - 1);
613
614     if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
615         return mpi_lset( X, 0 );
616
617     /*
618      * shift by count / limb_size
619      */
620     if( v0 > 0 )
621     {
622         for( i = 0; i < X->n - v0; i++ )
623             X->p[i] = X->p[i + v0];
624
625         for( ; i < X->n; i++ )
626             X->p[i] = 0;
627     }
628
629     /*
630      * shift by count % limb_size
631      */
632     if( v1 > 0 )
633     {
634         for( i = X->n; i > 0; i-- )
635         {
636             r1 = X->p[i - 1] << (biL - v1);
637             X->p[i - 1] >>= v1;
638             X->p[i - 1] |= r0;
639             r0 = r1;
640         }
641     }
642
643     return( 0 );
644 }
645
646 /*
647  * Compare unsigned values
648  */
649 int mpi_cmp_abs( const mpi *X, const mpi *Y )
650 {
651     size_t i, j;
652
653     for( i = X->n; i > 0; i-- )
654         if( X->p[i - 1] != 0 )
655             break;
656
657     for( j = Y->n; j > 0; j-- )
658         if( Y->p[j - 1] != 0 )
659             break;
660
661     if( i == 0 && j == 0 )
662         return( 0 );
663
664     if( i > j ) return(  1 );
665     if( j > i ) return( -1 );
666
667     for( ; i > 0; i-- )
668     {
669         if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
670         if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
671     }
672
673     return( 0 );
674 }
675
676 /*
677  * Compare signed values
678  */
679 int mpi_cmp_mpi( const mpi *X, const mpi *Y )
680 {
681     size_t i, j;
682
683     for( i = X->n; i > 0; i-- )
684         if( X->p[i - 1] != 0 )
685             break;
686
687     for( j = Y->n; j > 0; j-- )
688         if( Y->p[j - 1] != 0 )
689             break;
690
691     if( i == 0 && j == 0 )
692         return( 0 );
693
694     if( i > j ) return(  X->s );
695     if( j > i ) return( -Y->s );
696
697     if( X->s > 0 && Y->s < 0 ) return(  1 );
698     if( Y->s > 0 && X->s < 0 ) return( -1 );
699
700     for( ; i > 0; i-- )
701     {
702         if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
703         if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
704     }
705
706     return( 0 );
707 }
708
709 /*
710  * Compare signed values
711  */
712 int mpi_cmp_int( const mpi *X, t_sint z )
713 {
714     mpi Y;
715     t_uint p[1];
716
717     *p  = ( z < 0 ) ? -z : z;
718     Y.s = ( z < 0 ) ? -1 : 1;
719     Y.n = 1;
720     Y.p = p;
721
722     return( mpi_cmp_mpi( X, &Y ) );
723 }
724
725 /*
726  * Unsigned addition: X = |A| + |B|  (HAC 14.7)
727  */
728 int mpi_add_abs( mpi *X, const mpi *A, const mpi *B )
729 {
730     int ret;
731     size_t i, j;
732     t_uint *o, *p, c;
733
734     if( X == B )
735     {
736         const mpi *T = A; A = X; B = T;
737     }
738
739     if( X != A )
740         MPI_CHK( mpi_copy( X, A ) );
741    
742     /*
743      * X should always be positive as a result of unsigned additions.
744      */
745     X->s = 1;
746
747     for( j = B->n; j > 0; j-- )
748         if( B->p[j - 1] != 0 )
749             break;
750
751     MPI_CHK( mpi_grow( X, j ) );
752
753     o = B->p; p = X->p; c = 0;
754
755     for( i = 0; i < j; i++, o++, p++ )
756     {
757         *p +=  c; c  = ( *p <  c );
758         *p += *o; c += ( *p < *o );
759     }
760
761     while( c != 0 )
762     {
763         if( i >= X->n )
764         {
765             MPI_CHK( mpi_grow( X, i + 1 ) );
766             p = X->p + i;
767         }
768
769         *p += c; c = ( *p < c ); i++; p++;
770     }
771
772 cleanup:
773
774     return( ret );
775 }
776
777 /*
778  * Helper for mpi substraction
779  */
780 static void mpi_sub_hlp( size_t n, t_uint *s, t_uint *d )
781 {
782     size_t i;
783     t_uint c, z;
784
785     for( i = c = 0; i < n; i++, s++, d++ )
786     {
787         z = ( *d <  c );     *d -=  c;
788         c = ( *d < *s ) + z; *d -= *s;
789     }
790
791     while( c != 0 )
792     {
793         z = ( *d < c ); *d -= c;
794         c = z; i++; d++;
795     }
796 }
797
798 /*
799  * Unsigned substraction: X = |A| - |B|  (HAC 14.9)
800  */
801 int mpi_sub_abs( mpi *X, const mpi *A, const mpi *B )
802 {
803     mpi TB;
804     int ret;
805     size_t n;
806
807     if( mpi_cmp_abs( A, B ) < 0 )
808         return( POLARSSL_ERR_MPI_NEGATIVE_VALUE );
809
810     mpi_init( &TB );
811
812     if( X == B )
813     {
814         MPI_CHK( mpi_copy( &TB, B ) );
815         B = &TB;
816     }
817
818     if( X != A )
819         MPI_CHK( mpi_copy( X, A ) );
820
821     /*
822      * X should always be positive as a result of unsigned substractions.
823      */
824     X->s = 1;
825
826     ret = 0;
827
828     for( n = B->n; n > 0; n-- )
829         if( B->p[n - 1] != 0 )
830             break;
831
832     mpi_sub_hlp( n, B->p, X->p );
833
834 cleanup:
835
836     mpi_free( &TB );
837
838     return( ret );
839 }
840
841 /*
842  * Signed addition: X = A + B
843  */
844 int mpi_add_mpi( mpi *X, const mpi *A, const mpi *B )
845 {
846     int ret, s = A->s;
847
848     if( A->s * B->s < 0 )
849     {
850         if( mpi_cmp_abs( A, B ) >= 0 )
851         {
852             MPI_CHK( mpi_sub_abs( X, A, B ) );
853             X->s =  s;
854         }
855         else
856         {
857             MPI_CHK( mpi_sub_abs( X, B, A ) );
858             X->s = -s;
859         }
860     }
861     else
862     {
863         MPI_CHK( mpi_add_abs( X, A, B ) );
864         X->s = s;
865     }
866
867 cleanup:
868
869     return( ret );
870 }
871
872 /*
873  * Signed substraction: X = A - B
874  */
875 int mpi_sub_mpi( mpi *X, const mpi *A, const mpi *B )
876 {
877     int ret, s = A->s;
878
879     if( A->s * B->s > 0 )
880     {
881         if( mpi_cmp_abs( A, B ) >= 0 )
882         {
883             MPI_CHK( mpi_sub_abs( X, A, B ) );
884             X->s =  s;
885         }
886         else
887         {
888             MPI_CHK( mpi_sub_abs( X, B, A ) );
889             X->s = -s;
890         }
891     }
892     else
893     {
894         MPI_CHK( mpi_add_abs( X, A, B ) );
895         X->s = s;
896     }
897
898 cleanup:
899
900     return( ret );
901 }
902
903 /*
904  * Signed addition: X = A + b
905  */
906 int mpi_add_int( mpi *X, const mpi *A, t_sint b )
907 {
908     mpi _B;
909     t_uint p[1];
910
911     p[0] = ( b < 0 ) ? -b : b;
912     _B.s = ( b < 0 ) ? -1 : 1;
913     _B.n = 1;
914     _B.p = p;
915
916     return( mpi_add_mpi( X, A, &_B ) );
917 }
918
919 /*
920  * Signed substraction: X = A - b
921  */
922 int mpi_sub_int( mpi *X, const mpi *A, t_sint b )
923 {
924     mpi _B;
925     t_uint p[1];
926
927     p[0] = ( b < 0 ) ? -b : b;
928     _B.s = ( b < 0 ) ? -1 : 1;
929     _B.n = 1;
930     _B.p = p;
931
932     return( mpi_sub_mpi( X, A, &_B ) );
933 }
934
935 /*
936  * Helper for mpi multiplication
937  */ 
938 static void mpi_mul_hlp( size_t i, t_uint *s, t_uint *d, t_uint b )
939 {
940     t_uint c = 0, t = 0;
941
942 #if defined(MULADDC_HUIT)
943     for( ; i >= 8; i -= 8 )
944     {
945         MULADDC_INIT
946         MULADDC_HUIT
947         MULADDC_STOP
948     }
949
950     for( ; i > 0; i-- )
951     {
952         MULADDC_INIT
953         MULADDC_CORE
954         MULADDC_STOP
955     }
956 #else
957     for( ; i >= 16; i -= 16 )
958     {
959         MULADDC_INIT
960         MULADDC_CORE   MULADDC_CORE
961         MULADDC_CORE   MULADDC_CORE
962         MULADDC_CORE   MULADDC_CORE
963         MULADDC_CORE   MULADDC_CORE
964
965         MULADDC_CORE   MULADDC_CORE
966         MULADDC_CORE   MULADDC_CORE
967         MULADDC_CORE   MULADDC_CORE
968         MULADDC_CORE   MULADDC_CORE
969         MULADDC_STOP
970     }
971
972     for( ; i >= 8; i -= 8 )
973     {
974         MULADDC_INIT
975         MULADDC_CORE   MULADDC_CORE
976         MULADDC_CORE   MULADDC_CORE
977
978         MULADDC_CORE   MULADDC_CORE
979         MULADDC_CORE   MULADDC_CORE
980         MULADDC_STOP
981     }
982
983     for( ; i > 0; i-- )
984     {
985         MULADDC_INIT
986         MULADDC_CORE
987         MULADDC_STOP
988     }
989 #endif
990
991     t++;
992
993     do {
994         *d += c; c = ( *d < c ); d++;
995     }
996     while( c != 0 );
997 }
998
999 /*
1000  * Baseline multiplication: X = A * B  (HAC 14.12)
1001  */
1002 int mpi_mul_mpi( mpi *X, const mpi *A, const mpi *B )
1003 {
1004     int ret;
1005     size_t i, j;
1006     mpi TA, TB;
1007
1008     mpi_init( &TA ); mpi_init( &TB );
1009
1010     if( X == A ) { MPI_CHK( mpi_copy( &TA, A ) ); A = &TA; }
1011     if( X == B ) { MPI_CHK( mpi_copy( &TB, B ) ); B = &TB; }
1012
1013     for( i = A->n; i > 0; i-- )
1014         if( A->p[i - 1] != 0 )
1015             break;
1016
1017     for( j = B->n; j > 0; j-- )
1018         if( B->p[j - 1] != 0 )
1019             break;
1020
1021     MPI_CHK( mpi_grow( X, i + j ) );
1022     MPI_CHK( mpi_lset( X, 0 ) );
1023
1024     for( i++; j > 0; j-- )
1025         mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1026
1027     X->s = A->s * B->s;
1028
1029 cleanup:
1030
1031     mpi_free( &TB ); mpi_free( &TA );
1032
1033     return( ret );
1034 }
1035
1036 /*
1037  * Baseline multiplication: X = A * b
1038  */
1039 int mpi_mul_int( mpi *X, const mpi *A, t_sint b )
1040 {
1041     mpi _B;
1042     t_uint p[1];
1043
1044     _B.s = 1;
1045     _B.n = 1;
1046     _B.p = p;
1047     p[0] = b;
1048
1049     return( mpi_mul_mpi( X, A, &_B ) );
1050 }
1051
1052 /*
1053  * Division by mpi: A = Q * B + R  (HAC 14.20)
1054  */
1055 int mpi_div_mpi( mpi *Q, mpi *R, const mpi *A, const mpi *B )
1056 {
1057     int ret;
1058     size_t i, n, t, k;
1059     mpi X, Y, Z, T1, T2;
1060
1061     if( mpi_cmp_int( B, 0 ) == 0 )
1062         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1063
1064     mpi_init( &X ); mpi_init( &Y ); mpi_init( &Z );
1065     mpi_init( &T1 ); mpi_init( &T2 );
1066
1067     if( mpi_cmp_abs( A, B ) < 0 )
1068     {
1069         if( Q != NULL ) MPI_CHK( mpi_lset( Q, 0 ) );
1070         if( R != NULL ) MPI_CHK( mpi_copy( R, A ) );
1071         return( 0 );
1072     }
1073
1074     MPI_CHK( mpi_copy( &X, A ) );
1075     MPI_CHK( mpi_copy( &Y, B ) );
1076     X.s = Y.s = 1;
1077
1078     MPI_CHK( mpi_grow( &Z, A->n + 2 ) );
1079     MPI_CHK( mpi_lset( &Z,  0 ) );
1080     MPI_CHK( mpi_grow( &T1, 2 ) );
1081     MPI_CHK( mpi_grow( &T2, 3 ) );
1082
1083     k = mpi_msb( &Y ) % biL;
1084     if( k < biL - 1 )
1085     {
1086         k = biL - 1 - k;
1087         MPI_CHK( mpi_shift_l( &X, k ) );
1088         MPI_CHK( mpi_shift_l( &Y, k ) );
1089     }
1090     else k = 0;
1091
1092     n = X.n - 1;
1093     t = Y.n - 1;
1094     MPI_CHK( mpi_shift_l( &Y, biL * (n - t) ) );
1095
1096     while( mpi_cmp_mpi( &X, &Y ) >= 0 )
1097     {
1098         Z.p[n - t]++;
1099         mpi_sub_mpi( &X, &X, &Y );
1100     }
1101     mpi_shift_r( &Y, biL * (n - t) );
1102
1103     for( i = n; i > t ; i-- )
1104     {
1105         if( X.p[i] >= Y.p[t] )
1106             Z.p[i - t - 1] = ~0;
1107         else
1108         {
1109 #if defined(POLARSSL_HAVE_UDBL)
1110             t_udbl r;
1111
1112             r  = (t_udbl) X.p[i] << biL;
1113             r |= (t_udbl) X.p[i - 1];
1114             r /= Y.p[t];
1115             if( r > ((t_udbl) 1 << biL) - 1)
1116                 r = ((t_udbl) 1 << biL) - 1;
1117
1118             Z.p[i - t - 1] = (t_uint) r;
1119 #else
1120             /*
1121              * __udiv_qrnnd_c, from gmp/longlong.h
1122              */
1123             t_uint q0, q1, r0, r1;
1124             t_uint d0, d1, d, m;
1125
1126             d  = Y.p[t];
1127             d0 = ( d << biH ) >> biH;
1128             d1 = ( d >> biH );
1129
1130             q1 = X.p[i] / d1;
1131             r1 = X.p[i] - d1 * q1;
1132             r1 <<= biH;
1133             r1 |= ( X.p[i - 1] >> biH );
1134
1135             m = q1 * d0;
1136             if( r1 < m )
1137             {
1138                 q1--, r1 += d;
1139                 while( r1 >= d && r1 < m )
1140                     q1--, r1 += d;
1141             }
1142             r1 -= m;
1143
1144             q0 = r1 / d1;
1145             r0 = r1 - d1 * q0;
1146             r0 <<= biH;
1147             r0 |= ( X.p[i - 1] << biH ) >> biH;
1148
1149             m = q0 * d0;
1150             if( r0 < m )
1151             {
1152                 q0--, r0 += d;
1153                 while( r0 >= d && r0 < m )
1154                     q0--, r0 += d;
1155             }
1156             r0 -= m;
1157
1158             Z.p[i - t - 1] = ( q1 << biH ) | q0;
1159 #endif
1160         }
1161
1162         Z.p[i - t - 1]++;
1163         do
1164         {
1165             Z.p[i - t - 1]--;
1166
1167             MPI_CHK( mpi_lset( &T1, 0 ) );
1168             T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1169             T1.p[1] = Y.p[t];
1170             MPI_CHK( mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1171
1172             MPI_CHK( mpi_lset( &T2, 0 ) );
1173             T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1174             T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1175             T2.p[2] = X.p[i];
1176         }
1177         while( mpi_cmp_mpi( &T1, &T2 ) > 0 );
1178
1179         MPI_CHK( mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1180         MPI_CHK( mpi_shift_l( &T1,  biL * (i - t - 1) ) );
1181         MPI_CHK( mpi_sub_mpi( &X, &X, &T1 ) );
1182
1183         if( mpi_cmp_int( &X, 0 ) < 0 )
1184         {
1185             MPI_CHK( mpi_copy( &T1, &Y ) );
1186             MPI_CHK( mpi_shift_l( &T1, biL * (i - t - 1) ) );
1187             MPI_CHK( mpi_add_mpi( &X, &X, &T1 ) );
1188             Z.p[i - t - 1]--;
1189         }
1190     }
1191
1192     if( Q != NULL )
1193     {
1194         mpi_copy( Q, &Z );
1195         Q->s = A->s * B->s;
1196     }
1197
1198     if( R != NULL )
1199     {
1200         mpi_shift_r( &X, k );
1201         X.s = A->s;
1202         mpi_copy( R, &X );
1203
1204         if( mpi_cmp_int( R, 0 ) == 0 )
1205             R->s = 1;
1206     }
1207
1208 cleanup:
1209
1210     mpi_free( &X ); mpi_free( &Y ); mpi_free( &Z );
1211     mpi_free( &T1 ); mpi_free( &T2 );
1212
1213     return( ret );
1214 }
1215
1216 /*
1217  * Division by int: A = Q * b + R
1218  */
1219 int mpi_div_int( mpi *Q, mpi *R, const mpi *A, t_sint b )
1220 {
1221     mpi _B;
1222     t_uint p[1];
1223
1224     p[0] = ( b < 0 ) ? -b : b;
1225     _B.s = ( b < 0 ) ? -1 : 1;
1226     _B.n = 1;
1227     _B.p = p;
1228
1229     return( mpi_div_mpi( Q, R, A, &_B ) );
1230 }
1231
1232 /*
1233  * Modulo: R = A mod B
1234  */
1235 int mpi_mod_mpi( mpi *R, const mpi *A, const mpi *B )
1236 {
1237     int ret;
1238
1239     if( mpi_cmp_int( B, 0 ) < 0 )
1240         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1241
1242     MPI_CHK( mpi_div_mpi( NULL, R, A, B ) );
1243
1244     while( mpi_cmp_int( R, 0 ) < 0 )
1245       MPI_CHK( mpi_add_mpi( R, R, B ) );
1246
1247     while( mpi_cmp_mpi( R, B ) >= 0 )
1248       MPI_CHK( mpi_sub_mpi( R, R, B ) );
1249
1250 cleanup:
1251
1252     return( ret );
1253 }
1254
1255 /*
1256  * Modulo: r = A mod b
1257  */
1258 int mpi_mod_int( t_uint *r, const mpi *A, t_sint b )
1259 {
1260     size_t i;
1261     t_uint x, y, z;
1262
1263     if( b == 0 )
1264         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1265
1266     if( b < 0 )
1267         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1268
1269     /*
1270      * handle trivial cases
1271      */
1272     if( b == 1 )
1273     {
1274         *r = 0;
1275         return( 0 );
1276     }
1277
1278     if( b == 2 )
1279     {
1280         *r = A->p[0] & 1;
1281         return( 0 );
1282     }
1283
1284     /*
1285      * general case
1286      */
1287     for( i = A->n, y = 0; i > 0; i-- )
1288     {
1289         x  = A->p[i - 1];
1290         y  = ( y << biH ) | ( x >> biH );
1291         z  = y / b;
1292         y -= z * b;
1293
1294         x <<= biH;
1295         y  = ( y << biH ) | ( x >> biH );
1296         z  = y / b;
1297         y -= z * b;
1298     }
1299
1300     /*
1301      * If A is negative, then the current y represents a negative value.
1302      * Flipping it to the positive side.
1303      */
1304     if( A->s < 0 && y != 0 )
1305         y = b - y;
1306
1307     *r = y;
1308
1309     return( 0 );
1310 }
1311
1312 /*
1313  * Fast Montgomery initialization (thanks to Tom St Denis)
1314  */
1315 static void mpi_montg_init( t_uint *mm, const mpi *N )
1316 {
1317     t_uint x, m0 = N->p[0];
1318
1319     x  = m0;
1320     x += ( ( m0 + 2 ) & 4 ) << 1;
1321     x *= ( 2 - ( m0 * x ) );
1322
1323     if( biL >= 16 ) x *= ( 2 - ( m0 * x ) );
1324     if( biL >= 32 ) x *= ( 2 - ( m0 * x ) );
1325     if( biL >= 64 ) x *= ( 2 - ( m0 * x ) );
1326
1327     *mm = ~x + 1;
1328 }
1329
1330 /*
1331  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1332  */
1333 static void mpi_montmul( mpi *A, const mpi *B, const mpi *N, t_uint mm, const mpi *T )
1334 {
1335     size_t i, n, m;
1336     t_uint u0, u1, *d;
1337
1338     memset( T->p, 0, T->n * ciL );
1339
1340     d = T->p;
1341     n = N->n;
1342     m = ( B->n < n ) ? B->n : n;
1343
1344     for( i = 0; i < n; i++ )
1345     {
1346         /*
1347          * T = (T + u0*B + u1*N) / 2^biL
1348          */
1349         u0 = A->p[i];
1350         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1351
1352         mpi_mul_hlp( m, B->p, d, u0 );
1353         mpi_mul_hlp( n, N->p, d, u1 );
1354
1355         *d++ = u0; d[n + 1] = 0;
1356     }
1357
1358     memcpy( A->p, d, (n + 1) * ciL );
1359
1360     if( mpi_cmp_abs( A, N ) >= 0 )
1361         mpi_sub_hlp( n, N->p, A->p );
1362     else
1363         /* prevent timing attacks */
1364         mpi_sub_hlp( n, A->p, T->p );
1365 }
1366
1367 /*
1368  * Montgomery reduction: A = A * R^-1 mod N
1369  */
1370 static void mpi_montred( mpi *A, const mpi *N, t_uint mm, const mpi *T )
1371 {
1372     t_uint z = 1;
1373     mpi U;
1374
1375     U.n = U.s = (int) z;
1376     U.p = &z;
1377
1378     mpi_montmul( A, &U, N, mm, T );
1379 }
1380
1381 /*
1382  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1383  */
1384 int mpi_exp_mod( mpi *X, const mpi *A, const mpi *E, const mpi *N, mpi *_RR )
1385 {
1386     int ret;
1387     size_t wbits, wsize, one = 1;
1388     size_t i, j, nblimbs;
1389     size_t bufsize, nbits;
1390     t_uint ei, mm, state;
1391     mpi RR, T, W[ 2 << POLARSSL_MPI_WINDOW_SIZE ], Apos;
1392     int neg;
1393
1394     if( mpi_cmp_int( N, 0 ) < 0 || ( N->p[0] & 1 ) == 0 )
1395         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1396
1397     if( mpi_cmp_int( E, 0 ) < 0 )
1398         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1399
1400     /*
1401      * Init temps and window size
1402      */
1403     mpi_montg_init( &mm, N );
1404     mpi_init( &RR ); mpi_init( &T );
1405     memset( W, 0, sizeof( W ) );
1406
1407     i = mpi_msb( E );
1408
1409     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1410             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1411
1412     if( wsize > POLARSSL_MPI_WINDOW_SIZE )
1413         wsize = POLARSSL_MPI_WINDOW_SIZE;
1414
1415     j = N->n + 1;
1416     MPI_CHK( mpi_grow( X, j ) );
1417     MPI_CHK( mpi_grow( &W[1],  j ) );
1418     MPI_CHK( mpi_grow( &T, j * 2 ) );
1419
1420     /*
1421      * Compensate for negative A (and correct at the end)
1422      */
1423     neg = ( A->s == -1 );
1424
1425     mpi_init( &Apos );
1426     if( neg )
1427     {
1428         MPI_CHK( mpi_copy( &Apos, A ) );
1429         Apos.s = 1;
1430         A = &Apos;
1431     }
1432
1433     /*
1434      * If 1st call, pre-compute R^2 mod N
1435      */
1436     if( _RR == NULL || _RR->p == NULL )
1437     {
1438         MPI_CHK( mpi_lset( &RR, 1 ) );
1439         MPI_CHK( mpi_shift_l( &RR, N->n * 2 * biL ) );
1440         MPI_CHK( mpi_mod_mpi( &RR, &RR, N ) );
1441
1442         if( _RR != NULL )
1443             memcpy( _RR, &RR, sizeof( mpi ) );
1444     }
1445     else
1446         memcpy( &RR, _RR, sizeof( mpi ) );
1447
1448     /*
1449      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1450      */
1451     if( mpi_cmp_mpi( A, N ) >= 0 )
1452         mpi_mod_mpi( &W[1], A, N );
1453     else   mpi_copy( &W[1], A );
1454
1455     mpi_montmul( &W[1], &RR, N, mm, &T );
1456
1457     /*
1458      * X = R^2 * R^-1 mod N = R mod N
1459      */
1460     MPI_CHK( mpi_copy( X, &RR ) );
1461     mpi_montred( X, N, mm, &T );
1462
1463     if( wsize > 1 )
1464     {
1465         /*
1466          * W[1 << (wsize - 1)] = W[1] ^ (wsize - 1)
1467          */
1468         j =  one << (wsize - 1);
1469
1470         MPI_CHK( mpi_grow( &W[j], N->n + 1 ) );
1471         MPI_CHK( mpi_copy( &W[j], &W[1]    ) );
1472
1473         for( i = 0; i < wsize - 1; i++ )
1474             mpi_montmul( &W[j], &W[j], N, mm, &T );
1475     
1476         /*
1477          * W[i] = W[i - 1] * W[1]
1478          */
1479         for( i = j + 1; i < (one << wsize); i++ )
1480         {
1481             MPI_CHK( mpi_grow( &W[i], N->n + 1 ) );
1482             MPI_CHK( mpi_copy( &W[i], &W[i - 1] ) );
1483
1484             mpi_montmul( &W[i], &W[1], N, mm, &T );
1485         }
1486     }
1487
1488     nblimbs = E->n;
1489     bufsize = 0;
1490     nbits   = 0;
1491     wbits   = 0;
1492     state   = 0;
1493
1494     while( 1 )
1495     {
1496         if( bufsize == 0 )
1497         {
1498             if( nblimbs-- == 0 )
1499                 break;
1500
1501             bufsize = sizeof( t_uint ) << 3;
1502         }
1503
1504         bufsize--;
1505
1506         ei = (E->p[nblimbs] >> bufsize) & 1;
1507
1508         /*
1509          * skip leading 0s
1510          */
1511         if( ei == 0 && state == 0 )
1512             continue;
1513
1514         if( ei == 0 && state == 1 )
1515         {
1516             /*
1517              * out of window, square X
1518              */
1519             mpi_montmul( X, X, N, mm, &T );
1520             continue;
1521         }
1522
1523         /*
1524          * add ei to current window
1525          */
1526         state = 2;
1527
1528         nbits++;
1529         wbits |= (ei << (wsize - nbits));
1530
1531         if( nbits == wsize )
1532         {
1533             /*
1534              * X = X^wsize R^-1 mod N
1535              */
1536             for( i = 0; i < wsize; i++ )
1537                 mpi_montmul( X, X, N, mm, &T );
1538
1539             /*
1540              * X = X * W[wbits] R^-1 mod N
1541              */
1542             mpi_montmul( X, &W[wbits], N, mm, &T );
1543
1544             state--;
1545             nbits = 0;
1546             wbits = 0;
1547         }
1548     }
1549
1550     /*
1551      * process the remaining bits
1552      */
1553     for( i = 0; i < nbits; i++ )
1554     {
1555         mpi_montmul( X, X, N, mm, &T );
1556
1557         wbits <<= 1;
1558
1559         if( (wbits & (one << wsize)) != 0 )
1560             mpi_montmul( X, &W[1], N, mm, &T );
1561     }
1562
1563     /*
1564      * X = A^E * R * R^-1 mod N = A^E mod N
1565      */
1566     mpi_montred( X, N, mm, &T );
1567
1568     if( neg )
1569     {
1570         X->s = -1;
1571         mpi_add_mpi( X, N, X );
1572     }
1573
1574 cleanup:
1575
1576     for( i = (one << (wsize - 1)); i < (one << wsize); i++ )
1577         mpi_free( &W[i] );
1578
1579     mpi_free( &W[1] ); mpi_free( &T ); mpi_free( &Apos );
1580
1581     if( _RR == NULL )
1582         mpi_free( &RR );
1583
1584     return( ret );
1585 }
1586
1587 /*
1588  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1589  */
1590 int mpi_gcd( mpi *G, const mpi *A, const mpi *B )
1591 {
1592     int ret;
1593     size_t lz, lzt;
1594     mpi TG, TA, TB;
1595
1596     mpi_init( &TG ); mpi_init( &TA ); mpi_init( &TB );
1597
1598     MPI_CHK( mpi_copy( &TA, A ) );
1599     MPI_CHK( mpi_copy( &TB, B ) );
1600
1601     lz = mpi_lsb( &TA );
1602     lzt = mpi_lsb( &TB );
1603
1604     if ( lzt < lz )
1605         lz = lzt;
1606
1607     MPI_CHK( mpi_shift_r( &TA, lz ) );
1608     MPI_CHK( mpi_shift_r( &TB, lz ) );
1609
1610     TA.s = TB.s = 1;
1611
1612     while( mpi_cmp_int( &TA, 0 ) != 0 )
1613     {
1614         MPI_CHK( mpi_shift_r( &TA, mpi_lsb( &TA ) ) );
1615         MPI_CHK( mpi_shift_r( &TB, mpi_lsb( &TB ) ) );
1616
1617         if( mpi_cmp_mpi( &TA, &TB ) >= 0 )
1618         {
1619             MPI_CHK( mpi_sub_abs( &TA, &TA, &TB ) );
1620             MPI_CHK( mpi_shift_r( &TA, 1 ) );
1621         }
1622         else
1623         {
1624             MPI_CHK( mpi_sub_abs( &TB, &TB, &TA ) );
1625             MPI_CHK( mpi_shift_r( &TB, 1 ) );
1626         }
1627     }
1628
1629     MPI_CHK( mpi_shift_l( &TB, lz ) );
1630     MPI_CHK( mpi_copy( G, &TB ) );
1631
1632 cleanup:
1633
1634     mpi_free( &TG ); mpi_free( &TA ); mpi_free( &TB );
1635
1636     return( ret );
1637 }
1638
1639 int mpi_fill_random( mpi *X, size_t size,
1640                      int (*f_rng)(void *, unsigned char *, size_t),
1641                      void *p_rng )
1642 {
1643     int ret;
1644
1645     MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( size ) ) );
1646     MPI_CHK( mpi_lset( X, 0 ) );
1647
1648     MPI_CHK( f_rng( p_rng, (unsigned char *) X->p, size ) );
1649
1650 cleanup:
1651     return( ret );
1652 }
1653
1654 /*
1655  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1656  */
1657 int mpi_inv_mod( mpi *X, const mpi *A, const mpi *N )
1658 {
1659     int ret;
1660     mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1661
1662     if( mpi_cmp_int( N, 0 ) <= 0 )
1663         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1664
1665     mpi_init( &TA ); mpi_init( &TU ); mpi_init( &U1 ); mpi_init( &U2 );
1666     mpi_init( &G ); mpi_init( &TB ); mpi_init( &TV );
1667     mpi_init( &V1 ); mpi_init( &V2 );
1668
1669     MPI_CHK( mpi_gcd( &G, A, N ) );
1670
1671     if( mpi_cmp_int( &G, 1 ) != 0 )
1672     {
1673         ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
1674         goto cleanup;
1675     }
1676
1677     MPI_CHK( mpi_mod_mpi( &TA, A, N ) );
1678     MPI_CHK( mpi_copy( &TU, &TA ) );
1679     MPI_CHK( mpi_copy( &TB, N ) );
1680     MPI_CHK( mpi_copy( &TV, N ) );
1681
1682     MPI_CHK( mpi_lset( &U1, 1 ) );
1683     MPI_CHK( mpi_lset( &U2, 0 ) );
1684     MPI_CHK( mpi_lset( &V1, 0 ) );
1685     MPI_CHK( mpi_lset( &V2, 1 ) );
1686
1687     do
1688     {
1689         while( ( TU.p[0] & 1 ) == 0 )
1690         {
1691             MPI_CHK( mpi_shift_r( &TU, 1 ) );
1692
1693             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1694             {
1695                 MPI_CHK( mpi_add_mpi( &U1, &U1, &TB ) );
1696                 MPI_CHK( mpi_sub_mpi( &U2, &U2, &TA ) );
1697             }
1698
1699             MPI_CHK( mpi_shift_r( &U1, 1 ) );
1700             MPI_CHK( mpi_shift_r( &U2, 1 ) );
1701         }
1702
1703         while( ( TV.p[0] & 1 ) == 0 )
1704         {
1705             MPI_CHK( mpi_shift_r( &TV, 1 ) );
1706
1707             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1708             {
1709                 MPI_CHK( mpi_add_mpi( &V1, &V1, &TB ) );
1710                 MPI_CHK( mpi_sub_mpi( &V2, &V2, &TA ) );
1711             }
1712
1713             MPI_CHK( mpi_shift_r( &V1, 1 ) );
1714             MPI_CHK( mpi_shift_r( &V2, 1 ) );
1715         }
1716
1717         if( mpi_cmp_mpi( &TU, &TV ) >= 0 )
1718         {
1719             MPI_CHK( mpi_sub_mpi( &TU, &TU, &TV ) );
1720             MPI_CHK( mpi_sub_mpi( &U1, &U1, &V1 ) );
1721             MPI_CHK( mpi_sub_mpi( &U2, &U2, &V2 ) );
1722         }
1723         else
1724         {
1725             MPI_CHK( mpi_sub_mpi( &TV, &TV, &TU ) );
1726             MPI_CHK( mpi_sub_mpi( &V1, &V1, &U1 ) );
1727             MPI_CHK( mpi_sub_mpi( &V2, &V2, &U2 ) );
1728         }
1729     }
1730     while( mpi_cmp_int( &TU, 0 ) != 0 );
1731
1732     while( mpi_cmp_int( &V1, 0 ) < 0 )
1733         MPI_CHK( mpi_add_mpi( &V1, &V1, N ) );
1734
1735     while( mpi_cmp_mpi( &V1, N ) >= 0 )
1736         MPI_CHK( mpi_sub_mpi( &V1, &V1, N ) );
1737
1738     MPI_CHK( mpi_copy( X, &V1 ) );
1739
1740 cleanup:
1741
1742     mpi_free( &TA ); mpi_free( &TU ); mpi_free( &U1 ); mpi_free( &U2 );
1743     mpi_free( &G ); mpi_free( &TB ); mpi_free( &TV );
1744     mpi_free( &V1 ); mpi_free( &V2 );
1745
1746     return( ret );
1747 }
1748
1749 #if defined(POLARSSL_GENPRIME)
1750
1751 static const int small_prime[] =
1752 {
1753         3,    5,    7,   11,   13,   17,   19,   23,
1754        29,   31,   37,   41,   43,   47,   53,   59,
1755        61,   67,   71,   73,   79,   83,   89,   97,
1756       101,  103,  107,  109,  113,  127,  131,  137,
1757       139,  149,  151,  157,  163,  167,  173,  179,
1758       181,  191,  193,  197,  199,  211,  223,  227,
1759       229,  233,  239,  241,  251,  257,  263,  269,
1760       271,  277,  281,  283,  293,  307,  311,  313,
1761       317,  331,  337,  347,  349,  353,  359,  367,
1762       373,  379,  383,  389,  397,  401,  409,  419,
1763       421,  431,  433,  439,  443,  449,  457,  461,
1764       463,  467,  479,  487,  491,  499,  503,  509,
1765       521,  523,  541,  547,  557,  563,  569,  571,
1766       577,  587,  593,  599,  601,  607,  613,  617,
1767       619,  631,  641,  643,  647,  653,  659,  661,
1768       673,  677,  683,  691,  701,  709,  719,  727,
1769       733,  739,  743,  751,  757,  761,  769,  773,
1770       787,  797,  809,  811,  821,  823,  827,  829,
1771       839,  853,  857,  859,  863,  877,  881,  883,
1772       887,  907,  911,  919,  929,  937,  941,  947,
1773       953,  967,  971,  977,  983,  991,  997, -103
1774 };
1775
1776 /*
1777  * Miller-Rabin primality test  (HAC 4.24)
1778  */
1779 int mpi_is_prime( mpi *X,
1780                   int (*f_rng)(void *, unsigned char *, size_t),
1781                   void *p_rng )
1782 {
1783     int ret, xs;
1784     size_t i, j, n, s;
1785     mpi W, R, T, A, RR;
1786
1787     if( mpi_cmp_int( X, 0 ) == 0 ||
1788         mpi_cmp_int( X, 1 ) == 0 )
1789         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1790
1791     if( mpi_cmp_int( X, 2 ) == 0 )
1792         return( 0 );
1793
1794     mpi_init( &W ); mpi_init( &R ); mpi_init( &T ); mpi_init( &A );
1795     mpi_init( &RR );
1796
1797     xs = X->s; X->s = 1;
1798
1799     /*
1800      * test trivial factors first
1801      */
1802     if( ( X->p[0] & 1 ) == 0 )
1803         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1804
1805     for( i = 0; small_prime[i] > 0; i++ )
1806     {
1807         t_uint r;
1808
1809         if( mpi_cmp_int( X, small_prime[i] ) <= 0 )
1810             return( 0 );
1811
1812         MPI_CHK( mpi_mod_int( &r, X, small_prime[i] ) );
1813
1814         if( r == 0 )
1815             return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
1816     }
1817
1818     /*
1819      * W = |X| - 1
1820      * R = W >> lsb( W )
1821      */
1822     MPI_CHK( mpi_sub_int( &W, X, 1 ) );
1823     s = mpi_lsb( &W );
1824     MPI_CHK( mpi_copy( &R, &W ) );
1825     MPI_CHK( mpi_shift_r( &R, s ) );
1826
1827     i = mpi_msb( X );
1828     /*
1829      * HAC, table 4.4
1830      */
1831     n = ( ( i >= 1300 ) ?  2 : ( i >=  850 ) ?  3 :
1832           ( i >=  650 ) ?  4 : ( i >=  350 ) ?  8 :
1833           ( i >=  250 ) ? 12 : ( i >=  150 ) ? 18 : 27 );
1834
1835     for( i = 0; i < n; i++ )
1836     {
1837         /*
1838          * pick a random A, 1 < A < |X| - 1
1839          */
1840         MPI_CHK( mpi_fill_random( &A, X->n * ciL, f_rng, p_rng ) );
1841
1842         if( mpi_cmp_mpi( &A, &W ) >= 0 )
1843         {
1844             j = mpi_msb( &A ) - mpi_msb( &W );
1845             MPI_CHK( mpi_shift_r( &A, j + 1 ) );
1846         }
1847         A.p[0] |= 3;
1848
1849         /*
1850          * A = A^R mod |X|
1851          */
1852         MPI_CHK( mpi_exp_mod( &A, &A, &R, X, &RR ) );
1853
1854         if( mpi_cmp_mpi( &A, &W ) == 0 ||
1855             mpi_cmp_int( &A,  1 ) == 0 )
1856             continue;
1857
1858         j = 1;
1859         while( j < s && mpi_cmp_mpi( &A, &W ) != 0 )
1860         {
1861             /*
1862              * A = A * A mod |X|
1863              */
1864             MPI_CHK( mpi_mul_mpi( &T, &A, &A ) );
1865             MPI_CHK( mpi_mod_mpi( &A, &T, X  ) );
1866
1867             if( mpi_cmp_int( &A, 1 ) == 0 )
1868                 break;
1869
1870             j++;
1871         }
1872
1873         /*
1874          * not prime if A != |X| - 1 or A == 1
1875          */
1876         if( mpi_cmp_mpi( &A, &W ) != 0 ||
1877             mpi_cmp_int( &A,  1 ) == 0 )
1878         {
1879             ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
1880             break;
1881         }
1882     }
1883
1884 cleanup:
1885
1886     X->s = xs;
1887
1888     mpi_free( &W ); mpi_free( &R ); mpi_free( &T ); mpi_free( &A );
1889     mpi_free( &RR );
1890
1891     return( ret );
1892 }
1893
1894 /*
1895  * Prime number generation
1896  */
1897 int mpi_gen_prime( mpi *X, size_t nbits, int dh_flag,
1898                    int (*f_rng)(void *, unsigned char *, size_t),
1899                    void *p_rng )
1900 {
1901     int ret;
1902     size_t k, n;
1903     mpi Y;
1904
1905     if( nbits < 3 || nbits > POLARSSL_MPI_MAX_BITS )
1906         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1907
1908     mpi_init( &Y );
1909
1910     n = BITS_TO_LIMBS( nbits );
1911
1912     MPI_CHK( mpi_fill_random( X, n * ciL, f_rng, p_rng ) );
1913
1914     k = mpi_msb( X );
1915     if( k < nbits ) MPI_CHK( mpi_shift_l( X, nbits - k ) );
1916     if( k > nbits ) MPI_CHK( mpi_shift_r( X, k - nbits ) );
1917
1918     X->p[0] |= 3;
1919
1920     if( dh_flag == 0 )
1921     {
1922         while( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) != 0 )
1923         {
1924             if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1925                 goto cleanup;
1926
1927             MPI_CHK( mpi_add_int( X, X, 2 ) );
1928         }
1929     }
1930     else
1931     {
1932         MPI_CHK( mpi_sub_int( &Y, X, 1 ) );
1933         MPI_CHK( mpi_shift_r( &Y, 1 ) );
1934
1935         while( 1 )
1936         {
1937             if( ( ret = mpi_is_prime( X, f_rng, p_rng ) ) == 0 )
1938             {
1939                 if( ( ret = mpi_is_prime( &Y, f_rng, p_rng ) ) == 0 )
1940                     break;
1941
1942                 if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1943                     goto cleanup;
1944             }
1945
1946             if( ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE )
1947                 goto cleanup;
1948
1949             MPI_CHK( mpi_add_int( &Y, X, 1 ) );
1950             MPI_CHK( mpi_add_int(  X, X, 2 ) );
1951             MPI_CHK( mpi_shift_r( &Y, 1 ) );
1952         }
1953     }
1954
1955 cleanup:
1956
1957     mpi_free( &Y );
1958
1959     return( ret );
1960 }
1961
1962 #endif
1963
1964 #if defined(POLARSSL_SELF_TEST)
1965
1966 #define GCD_PAIR_COUNT  3
1967
1968 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
1969 {
1970     { 693, 609, 21 },
1971     { 1764, 868, 28 },
1972     { 768454923, 542167814, 1 }
1973 };
1974
1975 /*
1976  * Checkup routine
1977  */
1978 int mpi_self_test( int verbose )
1979 {
1980     int ret, i;
1981     mpi A, E, N, X, Y, U, V;
1982
1983     mpi_init( &A ); mpi_init( &E ); mpi_init( &N ); mpi_init( &X );
1984     mpi_init( &Y ); mpi_init( &U ); mpi_init( &V );
1985
1986     MPI_CHK( mpi_read_string( &A, 16,
1987         "EFE021C2645FD1DC586E69184AF4A31E" \
1988         "D5F53E93B5F123FA41680867BA110131" \
1989         "944FE7952E2517337780CB0DB80E61AA" \
1990         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
1991
1992     MPI_CHK( mpi_read_string( &E, 16,
1993         "B2E7EFD37075B9F03FF989C7C5051C20" \
1994         "34D2A323810251127E7BF8625A4F49A5" \
1995         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
1996         "5B5C25763222FEFCCFC38B832366C29E" ) );
1997
1998     MPI_CHK( mpi_read_string( &N, 16,
1999         "0066A198186C18C10B2F5ED9B522752A" \
2000         "9830B69916E535C8F047518A889A43A5" \
2001         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2002
2003     MPI_CHK( mpi_mul_mpi( &X, &A, &N ) );
2004
2005     MPI_CHK( mpi_read_string( &U, 16,
2006         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2007         "9E857EA95A03512E2BAE7391688D264A" \
2008         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2009         "8001B72E848A38CAE1C65F78E56ABDEF" \
2010         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2011         "ECF677152EF804370C1A305CAF3B5BF1" \
2012         "30879B56C61DE584A0F53A2447A51E" ) );
2013
2014     if( verbose != 0 )
2015         printf( "  MPI test #1 (mul_mpi): " );
2016
2017     if( mpi_cmp_mpi( &X, &U ) != 0 )
2018     {
2019         if( verbose != 0 )
2020             printf( "failed\n" );
2021
2022         return( 1 );
2023     }
2024
2025     if( verbose != 0 )
2026         printf( "passed\n" );
2027
2028     MPI_CHK( mpi_div_mpi( &X, &Y, &A, &N ) );
2029
2030     MPI_CHK( mpi_read_string( &U, 16,
2031         "256567336059E52CAE22925474705F39A94" ) );
2032
2033     MPI_CHK( mpi_read_string( &V, 16,
2034         "6613F26162223DF488E9CD48CC132C7A" \
2035         "0AC93C701B001B092E4E5B9F73BCD27B" \
2036         "9EE50D0657C77F374E903CDFA4C642" ) );
2037
2038     if( verbose != 0 )
2039         printf( "  MPI test #2 (div_mpi): " );
2040
2041     if( mpi_cmp_mpi( &X, &U ) != 0 ||
2042         mpi_cmp_mpi( &Y, &V ) != 0 )
2043     {
2044         if( verbose != 0 )
2045             printf( "failed\n" );
2046
2047         return( 1 );
2048     }
2049
2050     if( verbose != 0 )
2051         printf( "passed\n" );
2052
2053     MPI_CHK( mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2054
2055     MPI_CHK( mpi_read_string( &U, 16,
2056         "36E139AEA55215609D2816998ED020BB" \
2057         "BD96C37890F65171D948E9BC7CBAA4D9" \
2058         "325D24D6A3C12710F10A09FA08AB87" ) );
2059
2060     if( verbose != 0 )
2061         printf( "  MPI test #3 (exp_mod): " );
2062
2063     if( mpi_cmp_mpi( &X, &U ) != 0 )
2064     {
2065         if( verbose != 0 )
2066             printf( "failed\n" );
2067
2068         return( 1 );
2069     }
2070
2071     if( verbose != 0 )
2072         printf( "passed\n" );
2073
2074 #if defined(POLARSSL_GENPRIME)
2075     MPI_CHK( mpi_inv_mod( &X, &A, &N ) );
2076
2077     MPI_CHK( mpi_read_string( &U, 16,
2078         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2079         "C3DBA76456363A10869622EAC2DD84EC" \
2080         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2081
2082     if( verbose != 0 )
2083         printf( "  MPI test #4 (inv_mod): " );
2084
2085     if( mpi_cmp_mpi( &X, &U ) != 0 )
2086     {
2087         if( verbose != 0 )
2088             printf( "failed\n" );
2089
2090         return( 1 );
2091     }
2092
2093     if( verbose != 0 )
2094         printf( "passed\n" );
2095 #endif
2096
2097     if( verbose != 0 )
2098         printf( "  MPI test #5 (simple gcd): " );
2099
2100     for ( i = 0; i < GCD_PAIR_COUNT; i++)
2101     {
2102         MPI_CHK( mpi_lset( &X, gcd_pairs[i][0] ) );
2103         MPI_CHK( mpi_lset( &Y, gcd_pairs[i][1] ) );
2104
2105             MPI_CHK( mpi_gcd( &A, &X, &Y ) );
2106
2107             if( mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2108             {
2109                     if( verbose != 0 )
2110                             printf( "failed at %d\n", i );
2111
2112                     return( 1 );
2113             }
2114     }
2115
2116     if( verbose != 0 )
2117         printf( "passed\n" );
2118
2119 cleanup:
2120
2121     if( ret != 0 && verbose != 0 )
2122         printf( "Unexpected error, return code = %08X\n", ret );
2123
2124     mpi_free( &A ); mpi_free( &E ); mpi_free( &N ); mpi_free( &X );
2125     mpi_free( &Y ); mpi_free( &U ); mpi_free( &V );
2126
2127     if( verbose != 0 )
2128         printf( "\n" );
2129
2130     return( ret );
2131 }
2132
2133 #endif
2134
2135 #endif