7a7abfb00d262d1cdbe14bc2ca0283222ebd93b2
[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 static int mpi_cmp_abs_limbs (size_t n, const t_uint *p0, const t_uint *p1)
650 {
651     size_t i, j;
652
653     p0 += n;
654     for( i = n; i > 0; i-- )
655         if( *--p0 != 0 )
656             break;
657
658     p1 += n;
659     for( j = n; j > 0; j-- )
660         if( *--p1 != 0 )
661             break;
662
663     if( i == 0 && j == 0 )
664         return( 0 );
665
666     if( i > j ) return(  1 );
667     if( j > i ) return( -1 );
668
669     for( ; i > 0; i-- )
670     {
671         if( *p0 > *p1 ) return(  1 );
672         if( *p0 < *p1 ) return( -1 );
673         p0--; p1--;
674     }
675
676     return( 0 );
677 }
678
679 /*
680  * Compare unsigned values
681  */
682 int mpi_cmp_abs( const mpi *X, const mpi *Y )
683 {
684     size_t i, j;
685
686     for( i = X->n; i > 0; i-- )
687         if( X->p[i - 1] != 0 )
688             break;
689
690     for( j = Y->n; j > 0; j-- )
691         if( Y->p[j - 1] != 0 )
692             break;
693
694     if( i == 0 && j == 0 )
695         return( 0 );
696
697     if( i > j ) return(  1 );
698     if( j > i ) return( -1 );
699
700     for( ; i > 0; i-- )
701     {
702         if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
703         if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
704     }
705
706     return( 0 );
707 }
708
709 /*
710  * Compare signed values
711  */
712 int mpi_cmp_mpi( const mpi *X, const mpi *Y )
713 {
714     size_t i, j;
715
716     for( i = X->n; i > 0; i-- )
717         if( X->p[i - 1] != 0 )
718             break;
719
720     for( j = Y->n; j > 0; j-- )
721         if( Y->p[j - 1] != 0 )
722             break;
723
724     if( i == 0 && j == 0 )
725         return( 0 );
726
727     if( i > j ) return(  X->s );
728     if( j > i ) return( -Y->s );
729
730     if( X->s > 0 && Y->s < 0 ) return(  1 );
731     if( Y->s > 0 && X->s < 0 ) return( -1 );
732
733     for( ; i > 0; i-- )
734     {
735         if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
736         if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
737     }
738
739     return( 0 );
740 }
741
742 /*
743  * Compare signed values
744  */
745 int mpi_cmp_int( const mpi *X, t_sint z )
746 {
747     mpi Y;
748     t_uint p[1];
749
750     *p  = ( z < 0 ) ? -z : z;
751     Y.s = ( z < 0 ) ? -1 : 1;
752     Y.n = 1;
753     Y.p = p;
754
755     return( mpi_cmp_mpi( X, &Y ) );
756 }
757
758 /*
759  * Unsigned addition: X = |A| + |B|  (HAC 14.7)
760  */
761 int mpi_add_abs( mpi *X, const mpi *A, const mpi *B )
762 {
763     int ret;
764     size_t i, j;
765     t_uint *o, *p, c;
766
767     if( X == B )
768     {
769         const mpi *T = A; A = X; B = T;
770     }
771
772     if( X != A )
773         MPI_CHK( mpi_copy( X, A ) );
774    
775     /*
776      * X should always be positive as a result of unsigned additions.
777      */
778     X->s = 1;
779
780     for( j = B->n; j > 0; j-- )
781         if( B->p[j - 1] != 0 )
782             break;
783
784     MPI_CHK( mpi_grow( X, j ) );
785
786     o = B->p; p = X->p; c = 0;
787
788     for( i = 0; i < j; i++, o++, p++ )
789     {
790         *p +=  c; c  = ( *p <  c );
791         *p += *o; c += ( *p < *o );
792     }
793
794     while( c != 0 )
795     {
796         if( i >= X->n )
797         {
798             MPI_CHK( mpi_grow( X, i + 1 ) );
799             p = X->p + i;
800         }
801
802         *p += c; c = ( *p < c ); i++; p++;
803     }
804
805 cleanup:
806
807     return( ret );
808 }
809
810 /*
811  * Helper for mpi substraction
812  */
813 static t_uint mpi_sub_hlp( size_t n, const t_uint *s, t_uint *d )
814 {
815     size_t i;
816     t_uint c, z;
817
818     for( i = c = 0; i < n; i++, s++, d++ )
819     {
820         z = ( *d <  c );     *d -=  c;
821         c = ( *d < *s ) + z; *d -= *s;
822     }
823
824     return c;
825 }
826
827 /*
828  * Unsigned substraction: X = |A| - |B|  (HAC 14.9)
829  */
830 int mpi_sub_abs( mpi *X, const mpi *A, const mpi *B )
831 {
832     mpi TB;
833     int ret;
834     size_t n;
835     t_uint *d;
836     t_uint c, z;
837
838     if( mpi_cmp_abs( A, B ) < 0 )
839         return( POLARSSL_ERR_MPI_NEGATIVE_VALUE );
840
841     mpi_init( &TB );
842
843     if( X == B )
844     {
845         MPI_CHK( mpi_copy( &TB, B ) );
846         B = &TB;
847     }
848
849     if( X != A )
850         MPI_CHK( mpi_copy( X, A ) );
851
852     /*
853      * X should always be positive as a result of unsigned substractions.
854      */
855     X->s = 1;
856
857     ret = 0;
858
859     for( n = B->n; n > 0; n-- )
860         if( B->p[n - 1] != 0 )
861             break;
862
863     c = mpi_sub_hlp( n, B->p, X->p );
864     d = X->p + n;
865
866     while( c != 0 )
867     {
868         z = ( *d < c ); *d -= c;
869         c = z; d++;
870     }
871
872 cleanup:
873
874     mpi_free( &TB );
875
876     return( ret );
877 }
878
879 /*
880  * Signed addition: X = A + B
881  */
882 int mpi_add_mpi( mpi *X, const mpi *A, const mpi *B )
883 {
884     int ret, s = A->s;
885
886     if( A->s * B->s < 0 )
887     {
888         if( mpi_cmp_abs( A, B ) >= 0 )
889         {
890             MPI_CHK( mpi_sub_abs( X, A, B ) );
891             X->s =  s;
892         }
893         else
894         {
895             MPI_CHK( mpi_sub_abs( X, B, A ) );
896             X->s = -s;
897         }
898     }
899     else
900     {
901         MPI_CHK( mpi_add_abs( X, A, B ) );
902         X->s = s;
903     }
904
905 cleanup:
906
907     return( ret );
908 }
909
910 /*
911  * Signed substraction: X = A - B
912  */
913 int mpi_sub_mpi( mpi *X, const mpi *A, const mpi *B )
914 {
915     int ret, s = A->s;
916
917     if( A->s * B->s > 0 )
918     {
919         if( mpi_cmp_abs( A, B ) >= 0 )
920         {
921             MPI_CHK( mpi_sub_abs( X, A, B ) );
922             X->s =  s;
923         }
924         else
925         {
926             MPI_CHK( mpi_sub_abs( X, B, A ) );
927             X->s = -s;
928         }
929     }
930     else
931     {
932         MPI_CHK( mpi_add_abs( X, A, B ) );
933         X->s = s;
934     }
935
936 cleanup:
937
938     return( ret );
939 }
940
941 /*
942  * Signed addition: X = A + b
943  */
944 int mpi_add_int( mpi *X, const mpi *A, t_sint b )
945 {
946     mpi _B;
947     t_uint p[1];
948
949     p[0] = ( b < 0 ) ? -b : b;
950     _B.s = ( b < 0 ) ? -1 : 1;
951     _B.n = 1;
952     _B.p = p;
953
954     return( mpi_add_mpi( X, A, &_B ) );
955 }
956
957 /*
958  * Signed substraction: X = A - b
959  */
960 int mpi_sub_int( mpi *X, const mpi *A, t_sint b )
961 {
962     mpi _B;
963     t_uint p[1];
964
965     p[0] = ( b < 0 ) ? -b : b;
966     _B.s = ( b < 0 ) ? -1 : 1;
967     _B.n = 1;
968     _B.p = p;
969
970     return( mpi_sub_mpi( X, A, &_B ) );
971 }
972
973 /*
974  * Helper for mpi multiplication
975  */
976 static
977 #if defined(__APPLE__) && defined(__arm__)
978 /*
979  * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
980  * appears to need this to prevent bad ARM code generation at -O3.
981  */
982 __attribute__ ((noinline))
983 #endif
984 t_uint mpi_mul_hlp( size_t i, const t_uint *s, t_uint *d, t_uint b )
985 {
986     t_uint c = 0, t = 0;
987
988 #if defined(MULADDC_1024_LOOP)
989     MULADDC_1024_LOOP
990
991     for( ; i > 0; i-- )
992     {
993         MULADDC_INIT
994         MULADDC_CORE
995         MULADDC_STOP
996     }
997 #elif defined(MULADDC_HUIT)
998     for( ; i >= 8; i -= 8 )
999     {
1000         MULADDC_INIT
1001         MULADDC_HUIT
1002         MULADDC_STOP
1003     }
1004
1005     for( ; i > 0; i-- )
1006     {
1007         MULADDC_INIT
1008         MULADDC_CORE
1009         MULADDC_STOP
1010     }
1011 #else
1012     for( ; i >= 16; i -= 16 )
1013     {
1014         MULADDC_INIT
1015         MULADDC_CORE   MULADDC_CORE
1016         MULADDC_CORE   MULADDC_CORE
1017         MULADDC_CORE   MULADDC_CORE
1018         MULADDC_CORE   MULADDC_CORE
1019
1020         MULADDC_CORE   MULADDC_CORE
1021         MULADDC_CORE   MULADDC_CORE
1022         MULADDC_CORE   MULADDC_CORE
1023         MULADDC_CORE   MULADDC_CORE
1024         MULADDC_STOP
1025     }
1026
1027     for( ; i >= 8; i -= 8 )
1028     {
1029         MULADDC_INIT
1030         MULADDC_CORE   MULADDC_CORE
1031         MULADDC_CORE   MULADDC_CORE
1032
1033         MULADDC_CORE   MULADDC_CORE
1034         MULADDC_CORE   MULADDC_CORE
1035         MULADDC_STOP
1036     }
1037
1038     for( ; i > 0; i-- )
1039     {
1040         MULADDC_INIT
1041         MULADDC_CORE
1042         MULADDC_STOP
1043     }
1044 #endif
1045
1046     t++;
1047
1048     *d += c; c = ( *d < c );
1049     return c;
1050 }
1051
1052 /*
1053  * Help function of square: X = A * A  (HAC 14.16)
1054  * A is stored at the upper half of X.
1055  * Note that there was an error in the 1st printing of HAC.
1056  * Still, 5th printing has an error at step 2.3.
1057  *   The step 2.3 should be w[i+t] += c and handle overflow to w[i+t+1].
1058  *
1059  * N is the size (the number of limbs) of A.
1060  * D is a pointer to the limb of X.
1061  */
1062 static void
1063 mpi_sqr_hlp (size_t n, t_uint *d)
1064 {
1065   size_t i;
1066   t_uint c = 0;
1067
1068   for (i = 0; i < n; i++)
1069     {
1070       t_uint *wij = &d[i*2];
1071       t_uint *xj = &d[i+n];
1072       t_uint u, x_i;
1073
1074       x_i = *xj;
1075       *xj++ = c;
1076       asm ("mov    r8, #0\n\t"
1077            /* R8 := 0, the constant ZERO from here.  */
1078            "ldr    r9, [%[wij]]\n\t"          /* R9 := w_i_i; */
1079            "umull  r11, r12, %[x_i], %[x_i]\n\t"
1080            "adds   r9, r9, r11\n\t"
1081            "adc    %[u], r8, r12\n\t"
1082            "str    r9, [%[wij]], #4\n\t"
1083            "mov    %[c], r8\n\t"
1084            /* (C,U,R9) := w_i_i + x_i*x_i; w_i_i := R9; */
1085            "cmp    %[xj], %[xj_max]\n\t"
1086            "bcs    1f\n"
1087    "0:\n\t"
1088            /* (C,U,R9) := (C,U) + w_i_j + 2*x_i*x_j; */
1089            "ldr    r9, [%[wij]]\n\t"
1090            "adds   r9, r9, %[u]\n\t"
1091            "adc    %[u], %[c], r8\n\t"
1092            "ldr    r10, [%[xj]], #4\n\t"
1093            "umull  r11, r12, %[x_i], r10\n\t"
1094            "adds   r9, r9, r11\n\t"
1095            "adcs   %[u], %[u], r12\n\t"
1096            "adc    %[c], r8, r8\n\t"
1097            "adds   r9, r9, r11\n\t"
1098            "adcs   %[u], %[u], r12\n\t"
1099            "adc    %[c], %[c], r8\n\t"
1100            "str    r9, [%[wij]], #4\n\t"
1101            /**/
1102            "cmp    %[xj], %[xj_max]\n\t"
1103            "bcc    0b\n"
1104    "1:\n\t"
1105            "ldr    r9, [%[wij]]\n\t"
1106            "adds   %[u], %[u], r9\n\t"
1107            "adc    %[c], %[c], r8\n\t"
1108            "str    %[u], [%[wij]]"
1109            : [c] "=&r" (c), [u] "=&r" (u), [wij] "=r" (wij), [xj] "=r" (xj)
1110            : [x_i] "r" (x_i), [xj_max] "r" (&d[n*2]),
1111              "[wij]" (wij), "[xj]" (xj)
1112            : "r8", "r9", "r10", "r11", "r12", "memory", "cc" );
1113     }
1114 }
1115
1116 /*
1117  * Baseline multiplication: X = A * B  (HAC 14.12)
1118  */
1119 int mpi_mul_mpi( mpi *X, const mpi *A, const mpi *B )
1120 {
1121     int ret;
1122     size_t i, j, k;
1123     mpi TA, TB;
1124
1125     mpi_init( &TA ); mpi_init( &TB );
1126
1127     if( X == A ) { MPI_CHK( mpi_copy( &TA, A ) ); A = &TA; }
1128     if( X == B ) { MPI_CHK( mpi_copy( &TB, B ) ); B = &TB; }
1129
1130     for( i = A->n; i > 0; i-- )
1131         if( A->p[i - 1] != 0 )
1132             break;
1133
1134     for( j = B->n; j > 0; j-- )
1135         if( B->p[j - 1] != 0 )
1136             break;
1137
1138     MPI_CHK( mpi_grow( X, i + j ) );
1139     MPI_CHK( mpi_lset( X, 0 ) );
1140
1141     for(k = 0; k < j; k++ )
1142         mpi_mul_hlp( i, A->p, X->p + k, B->p[k]);
1143
1144     X->s = A->s * B->s;
1145
1146 cleanup:
1147
1148     mpi_free( &TB ); mpi_free( &TA );
1149
1150     return( ret );
1151 }
1152
1153 /*
1154  * Baseline multiplication: X = A * b
1155  */
1156 int mpi_mul_int( mpi *X, const mpi *A, t_sint b )
1157 {
1158     mpi _B;
1159     t_uint p[1];
1160
1161     _B.s = 1;
1162     _B.n = 1;
1163     _B.p = p;
1164     p[0] = b;
1165
1166     return( mpi_mul_mpi( X, A, &_B ) );
1167 }
1168
1169 /*
1170  * Division by mpi: A = Q * B + R  (HAC 14.20)
1171  */
1172 int mpi_div_mpi( mpi *Q, mpi *R, const mpi *A, const mpi *B )
1173 {
1174     int ret;
1175     size_t i, n, t, k;
1176     mpi X, Y, Z, T1, T2;
1177
1178     if( mpi_cmp_int( B, 0 ) == 0 )
1179         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1180
1181     mpi_init( &X ); mpi_init( &Y ); mpi_init( &Z );
1182     mpi_init( &T1 ); mpi_init( &T2 );
1183
1184     if( mpi_cmp_abs( A, B ) < 0 )
1185     {
1186         if( Q != NULL ) MPI_CHK( mpi_lset( Q, 0 ) );
1187         if( R != NULL ) MPI_CHK( mpi_copy( R, A ) );
1188         return( 0 );
1189     }
1190
1191     MPI_CHK( mpi_copy( &X, A ) );
1192     MPI_CHK( mpi_copy( &Y, B ) );
1193     X.s = Y.s = 1;
1194
1195     MPI_CHK( mpi_grow( &Z, A->n + 2 ) );
1196     MPI_CHK( mpi_lset( &Z,  0 ) );
1197     MPI_CHK( mpi_grow( &T1, 2 ) );
1198     MPI_CHK( mpi_grow( &T2, 3 ) );
1199
1200     k = mpi_msb( &Y ) % biL;
1201     if( k < biL - 1 )
1202     {
1203         k = biL - 1 - k;
1204         MPI_CHK( mpi_shift_l( &X, k ) );
1205         MPI_CHK( mpi_shift_l( &Y, k ) );
1206     }
1207     else k = 0;
1208
1209     n = X.n - 1;
1210     t = Y.n - 1;
1211     MPI_CHK( mpi_shift_l( &Y, biL * (n - t) ) );
1212
1213     while( mpi_cmp_mpi( &X, &Y ) >= 0 )
1214     {
1215         Z.p[n - t]++;
1216         mpi_sub_mpi( &X, &X, &Y );
1217     }
1218     mpi_shift_r( &Y, biL * (n - t) );
1219
1220     for( i = n; i > t ; i-- )
1221     {
1222         if( X.p[i] >= Y.p[t] )
1223             Z.p[i - t - 1] = ~0;
1224         else
1225         {
1226 #if defined(POLARSSL_HAVE_UDBL)
1227             t_udbl r;
1228
1229             r  = (t_udbl) X.p[i] << biL;
1230             r |= (t_udbl) X.p[i - 1];
1231             r /= Y.p[t];
1232             if( r > ((t_udbl) 1 << biL) - 1)
1233                 r = ((t_udbl) 1 << biL) - 1;
1234
1235             Z.p[i - t - 1] = (t_uint) r;
1236 #else
1237             /*
1238              * __udiv_qrnnd_c, from gmp/longlong.h
1239              */
1240             t_uint q0, q1, r0, r1;
1241             t_uint d0, d1, d, m;
1242
1243             d  = Y.p[t];
1244             d0 = ( d << biH ) >> biH;
1245             d1 = ( d >> biH );
1246
1247             q1 = X.p[i] / d1;
1248             r1 = X.p[i] - d1 * q1;
1249             r1 <<= biH;
1250             r1 |= ( X.p[i - 1] >> biH );
1251
1252             m = q1 * d0;
1253             if( r1 < m )
1254             {
1255                 q1--, r1 += d;
1256                 while( r1 >= d && r1 < m )
1257                     q1--, r1 += d;
1258             }
1259             r1 -= m;
1260
1261             q0 = r1 / d1;
1262             r0 = r1 - d1 * q0;
1263             r0 <<= biH;
1264             r0 |= ( X.p[i - 1] << biH ) >> biH;
1265
1266             m = q0 * d0;
1267             if( r0 < m )
1268             {
1269                 q0--, r0 += d;
1270                 while( r0 >= d && r0 < m )
1271                     q0--, r0 += d;
1272             }
1273             r0 -= m;
1274
1275             Z.p[i - t - 1] = ( q1 << biH ) | q0;
1276 #endif
1277         }
1278
1279         Z.p[i - t - 1]++;
1280         do
1281         {
1282             Z.p[i - t - 1]--;
1283
1284             MPI_CHK( mpi_lset( &T1, 0 ) );
1285             T1.p[0] = (t < 1) ? 0 : Y.p[t - 1];
1286             T1.p[1] = Y.p[t];
1287             MPI_CHK( mpi_mul_int( &T1, &T1, Z.p[i - t - 1] ) );
1288
1289             MPI_CHK( mpi_lset( &T2, 0 ) );
1290             T2.p[0] = (i < 2) ? 0 : X.p[i - 2];
1291             T2.p[1] = (i < 1) ? 0 : X.p[i - 1];
1292             T2.p[2] = X.p[i];
1293         }
1294         while( mpi_cmp_mpi( &T1, &T2 ) > 0 );
1295
1296         MPI_CHK( mpi_mul_int( &T1, &Y, Z.p[i - t - 1] ) );
1297         MPI_CHK( mpi_shift_l( &T1,  biL * (i - t - 1) ) );
1298         MPI_CHK( mpi_sub_mpi( &X, &X, &T1 ) );
1299
1300         if( mpi_cmp_int( &X, 0 ) < 0 )
1301         {
1302             MPI_CHK( mpi_copy( &T1, &Y ) );
1303             MPI_CHK( mpi_shift_l( &T1, biL * (i - t - 1) ) );
1304             MPI_CHK( mpi_add_mpi( &X, &X, &T1 ) );
1305             Z.p[i - t - 1]--;
1306         }
1307     }
1308
1309     if( Q != NULL )
1310     {
1311         mpi_copy( Q, &Z );
1312         Q->s = A->s * B->s;
1313     }
1314
1315     if( R != NULL )
1316     {
1317         mpi_shift_r( &X, k );
1318         X.s = A->s;
1319         mpi_copy( R, &X );
1320
1321         if( mpi_cmp_int( R, 0 ) == 0 )
1322             R->s = 1;
1323     }
1324
1325 cleanup:
1326
1327     mpi_free( &X ); mpi_free( &Y ); mpi_free( &Z );
1328     mpi_free( &T1 ); mpi_free( &T2 );
1329
1330     return( ret );
1331 }
1332
1333 /*
1334  * Division by int: A = Q * b + R
1335  */
1336 int mpi_div_int( mpi *Q, mpi *R, const mpi *A, t_sint b )
1337 {
1338     mpi _B;
1339     t_uint p[1];
1340
1341     p[0] = ( b < 0 ) ? -b : b;
1342     _B.s = ( b < 0 ) ? -1 : 1;
1343     _B.n = 1;
1344     _B.p = p;
1345
1346     return( mpi_div_mpi( Q, R, A, &_B ) );
1347 }
1348
1349 /*
1350  * Modulo: R = A mod B
1351  */
1352 int mpi_mod_mpi( mpi *R, const mpi *A, const mpi *B )
1353 {
1354     int ret;
1355
1356     if( mpi_cmp_int( B, 0 ) < 0 )
1357         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1358
1359     MPI_CHK( mpi_div_mpi( NULL, R, A, B ) );
1360
1361     while( mpi_cmp_int( R, 0 ) < 0 )
1362       MPI_CHK( mpi_add_mpi( R, R, B ) );
1363
1364     while( mpi_cmp_mpi( R, B ) >= 0 )
1365       MPI_CHK( mpi_sub_mpi( R, R, B ) );
1366
1367 cleanup:
1368
1369     return( ret );
1370 }
1371
1372 /*
1373  * Modulo: r = A mod b
1374  */
1375 int mpi_mod_int( t_uint *r, const mpi *A, t_sint b )
1376 {
1377     size_t i;
1378     t_uint x, y, z;
1379
1380     if( b == 0 )
1381         return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1382
1383     if( b < 0 )
1384         return POLARSSL_ERR_MPI_NEGATIVE_VALUE;
1385
1386     /*
1387      * handle trivial cases
1388      */
1389     if( b == 1 )
1390     {
1391         *r = 0;
1392         return( 0 );
1393     }
1394
1395     if( b == 2 )
1396     {
1397         *r = A->p[0] & 1;
1398         return( 0 );
1399     }
1400
1401     /*
1402      * general case
1403      */
1404     for( i = A->n, y = 0; i > 0; i-- )
1405     {
1406         x  = A->p[i - 1];
1407         y  = ( y << biH ) | ( x >> biH );
1408         z  = y / b;
1409         y -= z * b;
1410
1411         x <<= biH;
1412         y  = ( y << biH ) | ( x >> biH );
1413         z  = y / b;
1414         y -= z * b;
1415     }
1416
1417     /*
1418      * If A is negative, then the current y represents a negative value.
1419      * Flipping it to the positive side.
1420      */
1421     if( A->s < 0 && y != 0 )
1422         y = b - y;
1423
1424     *r = y;
1425
1426     return( 0 );
1427 }
1428
1429 /*
1430  * Fast Montgomery initialization (thanks to Tom St Denis)
1431  */
1432 static void mpi_montg_init( t_uint *mm, const mpi *N )
1433 {
1434     t_uint x, m0 = N->p[0];
1435
1436     x  = m0;
1437     x += ( ( m0 + 2 ) & 4 ) << 1;
1438     x *= ( 2 - ( m0 * x ) );
1439
1440     if( biL >= 16 ) x *= ( 2 - ( m0 * x ) );
1441     if( biL >= 32 ) x *= ( 2 - ( m0 * x ) );
1442     if( biL >= 64 ) x *= ( 2 - ( m0 * x ) );
1443
1444     *mm = ~x + 1;
1445 }
1446
1447 /*
1448  * Montgomery multiplication: A = A * B * R^-1 mod N  (HAC 14.36)
1449  * A is placed at the upper half of T.
1450  */
1451 static void mpi_montmul( const mpi *B, const mpi *N, t_uint mm, mpi *T )
1452 {
1453     size_t i, n, m;
1454     t_uint u0, u1, *d, c = 0;
1455
1456     d = T->p;
1457     n = N->n;
1458     m = ( B->n < n ) ? B->n : n;
1459
1460     for( i = 0; i < n; i++ )
1461     {
1462         /*
1463          * T = (T + u0*B + u1*N) / 2^biL
1464          */
1465         u0 = d[n];
1466         d[n] = c;
1467         u1 = ( d[0] + u0 * B->p[0] ) * mm;
1468
1469         mpi_mul_hlp( m, B->p, d, u0 );
1470         c = mpi_mul_hlp( n, N->p, d, u1 );
1471         d++;
1472     }
1473
1474     /* prevent timing attacks */
1475     if( ((mpi_cmp_abs_limbs ( n, d, N->p ) >= 0) | c) )
1476         mpi_sub_hlp( n, N->p, d );
1477     else
1478         mpi_sub_hlp( n, T->p, T->p);
1479 }
1480
1481 /*
1482  * Montgomery reduction: A = A * R^-1 mod N
1483  * A is placed at the upper half of T.
1484  */
1485 static void mpi_montred( const mpi *N, t_uint mm, mpi *T )
1486 {
1487     size_t i, j, n;
1488     t_uint u0, u1, *d, c = 0;
1489
1490     d = T->p;
1491     n = N->n;
1492
1493     for( i = 0; i < n; i++ )
1494     {
1495         /*
1496          * T = (T + u0 + u1*N) / 2^biL
1497          */
1498         u0 = d[n];
1499         d[n] = c;
1500         u1 = (d[0] + u0) * mm;
1501
1502         d[0] += u0;
1503         c = (d[0] < u0);
1504         for (j = 1; j < n; j++)
1505           {
1506             d[j] += c; c = ( d[j] < c );
1507           }
1508
1509         c = mpi_mul_hlp( n, N->p, d, u1 );
1510         d++;
1511     }
1512
1513     /* prevent timing attacks */
1514     if( ((mpi_cmp_abs_limbs ( n, d, N->p ) >= 0) | c) )
1515         mpi_sub_hlp( n, N->p, d );
1516     else
1517         mpi_sub_hlp( n, T->p, T->p);
1518 }
1519
1520 static void mpi_montsqr( mpi *X, const mpi *N, t_uint mm, mpi *T )
1521 {
1522 #if 1
1523     size_t n;
1524     t_uint c = 0, *d;
1525
1526     d = T->p;
1527     n = N->n;
1528     mpi_sqr_hlp (n, d);
1529
1530     while (d < T->p + n)
1531     {
1532         t_uint u, *d1;
1533
1534         u = *d * mm;
1535         c = mpi_mul_hlp( n, N->p, d, u );
1536         d1 = d + n + 1;
1537         while (d1 < T->p + n*2)
1538           {
1539             *d1 += c; c = ( *d1 < c );
1540             d1++;
1541           }
1542
1543         d++;
1544     }
1545
1546     /* prevent timing attacks */
1547     if( ((mpi_cmp_abs_limbs ( n, d, N->p ) >= 0) | c) )
1548         mpi_sub_hlp( n, N->p, d );
1549     else
1550         mpi_sub_hlp( n, T->p, T->p);
1551 #else
1552     memcpy ( X->p, T->p + N->n, N->n * ciL);
1553     mpi_montmul( X, N, mm, T );
1554 #endif
1555 }
1556
1557 /*
1558  * Sliding-window exponentiation: X = A^E mod N  (HAC 14.85)
1559  */
1560 int mpi_exp_mod( mpi *X, const mpi *A, const mpi *E, const mpi *N, mpi *_RR )
1561 {
1562     int ret;
1563     size_t wbits, wsize, one = 1;
1564     size_t i, j, nblimbs;
1565     size_t bufsize, nbits;
1566     t_uint ei, mm, state;
1567     mpi RR, T, W[ 2 << POLARSSL_MPI_WINDOW_SIZE ], Apos;
1568     int neg;
1569
1570     if( mpi_cmp_int( N, 0 ) < 0 || ( N->p[0] & 1 ) == 0 )
1571         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1572
1573     if( mpi_cmp_int( E, 0 ) < 0 )
1574         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1575
1576     /*
1577      * Init temps and window size
1578      */
1579     mpi_montg_init( &mm, N );
1580     mpi_init( &RR ); mpi_init( &T );
1581     memset( W, 0, sizeof( W ) );
1582
1583     i = mpi_msb( E );
1584
1585     wsize = ( i > 671 ) ? 6 : ( i > 239 ) ? 5 :
1586             ( i >  79 ) ? 4 : ( i >  23 ) ? 3 : 1;
1587
1588     if( wsize > POLARSSL_MPI_WINDOW_SIZE )
1589         wsize = POLARSSL_MPI_WINDOW_SIZE;
1590
1591     j = N->n;
1592     MPI_CHK( mpi_grow( X, N->n ) );
1593     MPI_CHK( mpi_grow( &W[1],  N->n ) );
1594     MPI_CHK( mpi_grow( &T, N->n * 2 ) ); /* T = 0 here.  */
1595
1596     /*
1597      * Compensate for negative A (and correct at the end)
1598      */
1599     neg = ( A->s == -1 );
1600
1601     mpi_init( &Apos );
1602     if( neg )
1603     {
1604         MPI_CHK( mpi_copy( &Apos, A ) );
1605         Apos.s = 1;
1606         A = &Apos;
1607     }
1608
1609     /*
1610      * If 1st call, pre-compute R^2 mod N
1611      */
1612     if( _RR == NULL || _RR->p == NULL )
1613     {
1614         /* T->p is all zero here. */
1615         mpi_sub_hlp( N->n, N->p, T.p + N->n);
1616         MPI_CHK( mpi_mod_mpi( &RR, &T, N ) );
1617
1618         if( _RR != NULL )
1619             memcpy( _RR, &RR, sizeof( mpi ) );
1620
1621         /* The condition of "the lower half of T is all zero" is kept. */
1622     }
1623     else
1624         memcpy( &RR, _RR, sizeof( mpi ) );
1625
1626     /*
1627      * W[1] = A * R^2 * R^-1 mod N = A * R mod N
1628      */
1629     if( mpi_cmp_mpi( A, N ) >= 0 )
1630         mpi_mod_mpi( &W[1], A, N );
1631     else   mpi_copy( &W[1], A );
1632
1633     memcpy ( T.p + N->n, W[1].p, N->n * ciL);
1634     mpi_montmul( &RR, N, mm, &T );
1635     memcpy ( W[1].p, T.p + N->n, N->n * ciL);
1636
1637     if( wsize > 1 )
1638     {
1639         /*
1640          * W[1 << (wsize - 1)] = W[1] ^ ( 2 ^ (wsize - 1) )
1641          */
1642         j =  one << (wsize - 1);
1643
1644         MPI_CHK( mpi_grow( &W[j], N->n  ) );
1645
1646         for( i = 0; i < wsize - 1; i++ )
1647             mpi_montsqr( X, N, mm, &T );
1648         memcpy ( W[j].p, T.p + N->n, N->n * ciL);
1649
1650         /*
1651          * W[i] = W[i - 1] * W[1]
1652          */
1653         for( i = j + 1; i < (one << wsize); i++ )
1654         {
1655             MPI_CHK( mpi_grow( &W[i], N->n      ) );
1656
1657             mpi_montmul( &W[1], N, mm, &T );
1658             memcpy ( W[i].p, T.p + N->n, N->n * ciL);
1659         }
1660     }
1661
1662     /*
1663      * X = R^2 * R^-1 mod N = R mod N
1664      */
1665     memcpy ( T.p + N->n, RR.p, N->n * ciL);
1666     mpi_montred( N, mm, &T );
1667
1668     nblimbs = E->n;
1669     bufsize = 0;
1670     nbits   = 0;
1671     wbits   = 0;
1672     state   = 0;
1673
1674     while( 1 )
1675     {
1676         if( bufsize == 0 )
1677         {
1678             if( nblimbs-- == 0 )
1679                 break;
1680
1681             bufsize = sizeof( t_uint ) << 3;
1682         }
1683
1684         bufsize--;
1685
1686         ei = (E->p[nblimbs] >> bufsize) & 1;
1687
1688         /*
1689          * skip leading 0s
1690          */
1691         if( ei == 0 && state == 0 )
1692             continue;
1693
1694         if( ei == 0 && state == 1 )
1695         {
1696             /*
1697              * out of window, square X
1698              */
1699             mpi_montsqr( X, N, mm, &T );
1700             continue;
1701         }
1702
1703         /*
1704          * add ei to current window
1705          */
1706         state = 2;
1707
1708         nbits++;
1709         wbits |= (ei << (wsize - nbits));
1710
1711         if( nbits == wsize )
1712         {
1713             /*
1714              * X = X^wsize R^-1 mod N
1715              */
1716             for( i = 0; i < wsize; i++ )
1717                 mpi_montsqr( X, N, mm, &T );
1718
1719             /*
1720              * X = X * W[wbits] R^-1 mod N
1721              */
1722             mpi_montmul( &W[wbits], N, mm, &T );
1723
1724             state--;
1725             nbits = 0;
1726             wbits = 0;
1727         }
1728     }
1729
1730     /*
1731      * process the remaining bits
1732      */
1733     for( i = 0; i < nbits; i++ )
1734     {
1735         mpi_montsqr( X, N, mm, &T );
1736
1737         wbits <<= 1;
1738
1739         if( (wbits & (one << wsize)) != 0 )
1740             mpi_montmul( &W[1], N, mm, &T );
1741     }
1742
1743     /*
1744      * X = A^E * R * R^-1 mod N = A^E mod N
1745      */
1746     mpi_montred( N, mm, &T );
1747     memcpy ( X->p, T.p + N->n, N->n * ciL);
1748
1749     if( neg )
1750     {
1751         X->s = -1;
1752         mpi_add_mpi( X, N, X );
1753     }
1754
1755 cleanup:
1756
1757     for( i = (one << (wsize - 1)); i < (one << wsize); i++ )
1758         mpi_free( &W[i] );
1759
1760     mpi_free( &W[1] ); mpi_free( &T ); mpi_free( &Apos );
1761
1762     if( _RR == NULL )
1763         mpi_free( &RR );
1764
1765     return( ret );
1766 }
1767
1768 /*
1769  * Greatest common divisor: G = gcd(A, B)  (HAC 14.54)
1770  */
1771 int mpi_gcd( mpi *G, const mpi *A, const mpi *B )
1772 {
1773     int ret;
1774     size_t lz, lzt;
1775     mpi TG, TA, TB;
1776
1777     mpi_init( &TG ); mpi_init( &TA ); mpi_init( &TB );
1778
1779     MPI_CHK( mpi_copy( &TA, A ) );
1780     MPI_CHK( mpi_copy( &TB, B ) );
1781
1782     lz = mpi_lsb( &TA );
1783     lzt = mpi_lsb( &TB );
1784
1785     if ( lzt < lz )
1786         lz = lzt;
1787
1788     MPI_CHK( mpi_shift_r( &TA, lz ) );
1789     MPI_CHK( mpi_shift_r( &TB, lz ) );
1790
1791     TA.s = TB.s = 1;
1792
1793     while( mpi_cmp_int( &TA, 0 ) != 0 )
1794     {
1795         MPI_CHK( mpi_shift_r( &TA, mpi_lsb( &TA ) ) );
1796         MPI_CHK( mpi_shift_r( &TB, mpi_lsb( &TB ) ) );
1797
1798         if( mpi_cmp_mpi( &TA, &TB ) >= 0 )
1799         {
1800             MPI_CHK( mpi_sub_abs( &TA, &TA, &TB ) );
1801             MPI_CHK( mpi_shift_r( &TA, 1 ) );
1802         }
1803         else
1804         {
1805             MPI_CHK( mpi_sub_abs( &TB, &TB, &TA ) );
1806             MPI_CHK( mpi_shift_r( &TB, 1 ) );
1807         }
1808     }
1809
1810     MPI_CHK( mpi_shift_l( &TB, lz ) );
1811     MPI_CHK( mpi_copy( G, &TB ) );
1812
1813 cleanup:
1814
1815     mpi_free( &TG ); mpi_free( &TA ); mpi_free( &TB );
1816
1817     return( ret );
1818 }
1819
1820 int mpi_fill_random( mpi *X, size_t size,
1821                      int (*f_rng)(void *, unsigned char *, size_t),
1822                      void *p_rng )
1823 {
1824     int ret;
1825
1826     MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( size ) ) );
1827     MPI_CHK( mpi_lset( X, 0 ) );
1828
1829     MPI_CHK( f_rng( p_rng, (unsigned char *) X->p, size ) );
1830
1831 cleanup:
1832     return( ret );
1833 }
1834
1835 /*
1836  * Modular inverse: X = A^-1 mod N  (HAC 14.61 / 14.64)
1837  */
1838 int mpi_inv_mod( mpi *X, const mpi *A, const mpi *N )
1839 {
1840     int ret;
1841     mpi G, TA, TU, U1, U2, TB, TV, V1, V2;
1842
1843     if( mpi_cmp_int( N, 0 ) <= 0 )
1844         return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
1845
1846     mpi_init( &TA ); mpi_init( &TU ); mpi_init( &U1 ); mpi_init( &U2 );
1847     mpi_init( &G ); mpi_init( &TB ); mpi_init( &TV );
1848     mpi_init( &V1 ); mpi_init( &V2 );
1849
1850     MPI_CHK( mpi_gcd( &G, A, N ) );
1851
1852     if( mpi_cmp_int( &G, 1 ) != 0 )
1853     {
1854         ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
1855         goto cleanup;
1856     }
1857
1858     MPI_CHK( mpi_mod_mpi( &TA, A, N ) );
1859     MPI_CHK( mpi_copy( &TU, &TA ) );
1860     MPI_CHK( mpi_copy( &TB, N ) );
1861     MPI_CHK( mpi_copy( &TV, N ) );
1862
1863     MPI_CHK( mpi_lset( &U1, 1 ) );
1864     MPI_CHK( mpi_lset( &U2, 0 ) );
1865     MPI_CHK( mpi_lset( &V1, 0 ) );
1866     MPI_CHK( mpi_lset( &V2, 1 ) );
1867
1868     do
1869     {
1870         while( ( TU.p[0] & 1 ) == 0 )
1871         {
1872             MPI_CHK( mpi_shift_r( &TU, 1 ) );
1873
1874             if( ( U1.p[0] & 1 ) != 0 || ( U2.p[0] & 1 ) != 0 )
1875             {
1876                 MPI_CHK( mpi_add_mpi( &U1, &U1, &TB ) );
1877                 MPI_CHK( mpi_sub_mpi( &U2, &U2, &TA ) );
1878             }
1879
1880             MPI_CHK( mpi_shift_r( &U1, 1 ) );
1881             MPI_CHK( mpi_shift_r( &U2, 1 ) );
1882         }
1883
1884         while( ( TV.p[0] & 1 ) == 0 )
1885         {
1886             MPI_CHK( mpi_shift_r( &TV, 1 ) );
1887
1888             if( ( V1.p[0] & 1 ) != 0 || ( V2.p[0] & 1 ) != 0 )
1889             {
1890                 MPI_CHK( mpi_add_mpi( &V1, &V1, &TB ) );
1891                 MPI_CHK( mpi_sub_mpi( &V2, &V2, &TA ) );
1892             }
1893
1894             MPI_CHK( mpi_shift_r( &V1, 1 ) );
1895             MPI_CHK( mpi_shift_r( &V2, 1 ) );
1896         }
1897
1898         if( mpi_cmp_mpi( &TU, &TV ) >= 0 )
1899         {
1900             MPI_CHK( mpi_sub_mpi( &TU, &TU, &TV ) );
1901             MPI_CHK( mpi_sub_mpi( &U1, &U1, &V1 ) );
1902             MPI_CHK( mpi_sub_mpi( &U2, &U2, &V2 ) );
1903         }
1904         else
1905         {
1906             MPI_CHK( mpi_sub_mpi( &TV, &TV, &TU ) );
1907             MPI_CHK( mpi_sub_mpi( &V1, &V1, &U1 ) );
1908             MPI_CHK( mpi_sub_mpi( &V2, &V2, &U2 ) );
1909         }
1910     }
1911     while( mpi_cmp_int( &TU, 0 ) != 0 );
1912
1913     while( mpi_cmp_int( &V1, 0 ) < 0 )
1914         MPI_CHK( mpi_add_mpi( &V1, &V1, N ) );
1915
1916     while( mpi_cmp_mpi( &V1, N ) >= 0 )
1917         MPI_CHK( mpi_sub_mpi( &V1, &V1, N ) );
1918
1919     MPI_CHK( mpi_copy( X, &V1 ) );
1920
1921 cleanup:
1922
1923     mpi_free( &TA ); mpi_free( &TU ); mpi_free( &U1 ); mpi_free( &U2 );
1924     mpi_free( &G ); mpi_free( &TB ); mpi_free( &TV );
1925     mpi_free( &V1 ); mpi_free( &V2 );
1926
1927     return( ret );
1928 }
1929
1930 #if defined(POLARSSL_GENPRIME)
1931
1932 static const int small_prime[] =
1933 {
1934 #if 0
1935         3,    5,    7,   11,   13,   17,   19,   23,
1936        29,   31,   37,   41,   43,   47,   53,   59,
1937        61,   67,   71,   73,   79,   83,   89,   97,
1938       101,  103,  107,  109,  113,  127,  131,  137,
1939       139,  149,  151,  157,  163,  167,  173,  179,
1940       181,  191,  193,  197,  199,  211,  223,  227,
1941       229,  233,  239,  241,  251,  257,  263,  269,
1942       271,  277,  281,  283,  293,  307,  311,  313,
1943       317,  331,  337,  347,  349,  353,  359,  367,
1944       373,  379,  383,  389,  397,  401,  409,  419,
1945       421,  431,  433,  439,  443,  449,  457,  461,
1946       463,  467,  479,  487,  491,  499,  503,  509,
1947       521,  523,  541,  547,  557,  563,  569,  571,
1948       577,  587,  593,  599,  601,  607,  613,  617,
1949       619,  631,  641,  643,  647,  653,  659,  661,
1950       673,  677,  683,  691,  701,
1951 #else
1952        97,
1953 #endif
1954                                     709,  719,  727,
1955       733,  739,  743,  751,  757,  761,  769,  773,
1956       787,  
1957 #if 0
1958             797,
1959 #endif
1960                   809,  811,  821,  823,  827,  829,
1961       839,  853,  857,  859,  863,  877,  881,  883,
1962       887,  907,  911,  919,  929,  937,  941,  947,
1963       953,  967,  971,  977,  983,  991,  997, 
1964 #if 1
1965                                                1009, 
1966      1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051,
1967      1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103,
1968      1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171,
1969      1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229,
1970      1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289,
1971      1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327,
1972      1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427,
1973      1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471,
1974      1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523,
1975      1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
1976      1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621,
1977      1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697,
1978      1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753,
1979      1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823,
1980      1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879,
1981      1889,
1982 #endif
1983
1984      -103
1985 };
1986
1987 /*
1988  * From Public domain code of JKISS RNG.
1989  *
1990  * Reference: David Jones, UCL Bioinformatics Group
1991  * Good Practice in (Pseudo) Random Number Generation for
1992  * Bioinformatics Applications
1993  *
1994  */
1995 struct jkiss_state { uint32_t x, y, z, c; };
1996 static struct jkiss_state jkiss_state_v;
1997
1998 int prng_seed (int (*f_rng)(void *, unsigned char *, size_t),
1999                void *p_rng)
2000 {
2001   int ret;
2002
2003   struct jkiss_state *s = &jkiss_state_v;
2004
2005   MPI_CHK ( f_rng (p_rng, (unsigned char *)s, sizeof (struct jkiss_state)) );
2006   while (s->y == 0)
2007     MPI_CHK ( f_rng (p_rng, (unsigned char *)&s->y, sizeof (uint32_t)) );
2008   s->z |= 1;                    /* avoiding z=c=0 */
2009
2010 cleanup:
2011   return ret;
2012 }
2013
2014 static uint32_t
2015 jkiss (struct jkiss_state *s)
2016 {
2017   uint64_t t;
2018
2019   s->x = 314527869 * s->x + 1234567;
2020   s->y ^= s->y << 5;
2021   s->y ^= s->y >> 7;
2022   s->y ^= s->y << 22;
2023   t = 4294584393ULL * s->z + s->c;
2024   s->c = (uint32_t)(t >> 32);
2025   s->z = (uint32_t)t;
2026
2027   return s->x + s->y + s->z;
2028 }
2029
2030 static int mpi_fill_pseudo_random ( mpi *X, size_t size)
2031 {
2032   int ret;
2033   uint32_t *p;
2034
2035   MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( size ) ) );
2036   MPI_CHK( mpi_lset( X, 0 ) );
2037
2038   /* Assume little endian.  */
2039   p = X->p;
2040   while (p < X->p + (size/ciL))
2041     *p++ = jkiss (&jkiss_state_v);
2042   if ((size % ciL))
2043     *p = jkiss (&jkiss_state_v) & ((1 << (8*(size % ciL))) - 1);
2044
2045 cleanup:
2046   return ret;
2047 }
2048
2049 /*
2050  * Miller-Rabin primality test  (HAC 4.24)
2051  */
2052 static
2053 int mpi_is_prime( mpi *X)
2054 {
2055     int ret, xs;
2056     size_t i, j, n, s;
2057     mpi W, R, T, A, RR;
2058
2059     if( mpi_cmp_int( X, 0 ) == 0 ||
2060         mpi_cmp_int( X, 1 ) == 0 )
2061         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
2062
2063     if( mpi_cmp_int( X, 2 ) == 0 )
2064         return( 0 );
2065
2066     mpi_init( &W ); mpi_init( &R ); mpi_init( &T ); mpi_init( &A );
2067     mpi_init( &RR );
2068
2069     xs = X->s; X->s = 1;
2070     ret = 0;
2071
2072 #if 0
2073     /*
2074      * test trivial factors first
2075      */
2076     if( ( X->p[0] & 1 ) == 0 )
2077         return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
2078 #endif
2079
2080     for( i = 0; small_prime[i] > 0; i++ )
2081     {
2082         t_uint r;
2083
2084         if( mpi_cmp_int( X, small_prime[i] ) <= 0 )
2085             return( 0 );
2086
2087         MPI_CHK( mpi_mod_int( &r, X, small_prime[i] ) );
2088
2089         if( r == 0 )
2090             return( POLARSSL_ERR_MPI_NOT_ACCEPTABLE );
2091     }
2092
2093     /*
2094      * W = |X| - 1
2095      * R = W >> lsb( W )
2096      */
2097     MPI_CHK( mpi_sub_int( &W, X, 1 ) );
2098     s = mpi_lsb( &W );
2099     MPI_CHK( mpi_copy( &R, &W ) );
2100     MPI_CHK( mpi_shift_r( &R, s ) );
2101     i = mpi_msb( X );
2102
2103     /* Fermat primality test with 2.  */
2104     mpi_lset (&T, 2);
2105     MPI_CHK( mpi_exp_mod( &T, &T, &W, X, &RR ) );
2106     if ( mpi_cmp_int (&T, 1) != 0)
2107       {
2108         ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
2109         goto cleanup;
2110       }
2111
2112
2113     /*
2114      * HAC, table 4.4
2115      */
2116     n = ( ( i >= 1300 ) ?  2 : ( i >=  850 ) ?  3 :
2117           ( i >=  650 ) ?  4 : ( i >=  350 ) ?  8 :
2118           ( i >=  250 ) ? 12 : ( i >=  150 ) ? 18 : 27 );
2119
2120     for( i = 0; i < n; i++ )
2121     {
2122         /*
2123          * pick a random A, 1 < A < |X| - 1
2124          */
2125         MPI_CHK( mpi_fill_pseudo_random( &A, X->n * ciL ) );
2126
2127         if( mpi_cmp_mpi( &A, &W ) >= 0 )
2128         {
2129             j = mpi_msb( &A ) - mpi_msb( &W );
2130             MPI_CHK( mpi_shift_r( &A, j + 1 ) );
2131         }
2132         A.p[0] |= 3;
2133
2134         /*
2135          * A = A^R mod |X|
2136          */
2137         MPI_CHK( mpi_exp_mod( &A, &A, &R, X, &RR ) );
2138
2139         if( mpi_cmp_mpi( &A, &W ) == 0 ||
2140             mpi_cmp_int( &A,  1 ) == 0 )
2141             continue;
2142
2143         j = 1;
2144         while( j < s && mpi_cmp_mpi( &A, &W ) != 0 )
2145         {
2146             /*
2147              * A = A * A mod |X|
2148              */
2149             MPI_CHK( mpi_mul_mpi( &T, &A, &A ) );
2150             MPI_CHK( mpi_mod_mpi( &A, &T, X  ) );
2151
2152             if( mpi_cmp_int( &A, 1 ) == 0 )
2153                 break;
2154
2155             j++;
2156         }
2157
2158         /*
2159          * not prime if A != |X| - 1 or A == 1
2160          */
2161         if( mpi_cmp_mpi( &A, &W ) != 0 ||
2162             mpi_cmp_int( &A,  1 ) == 0 )
2163         {
2164             ret = POLARSSL_ERR_MPI_NOT_ACCEPTABLE;
2165             break;
2166         }
2167     }
2168
2169 cleanup:
2170
2171     X->s = xs;
2172
2173     mpi_free( &W ); mpi_free( &R ); mpi_free( &T ); mpi_free( &A );
2174     mpi_free( &RR );
2175
2176     return( ret );
2177 }
2178
2179
2180 /*
2181  * Value M: multiply all primes up to 701 (except 97) and 797
2182  * (so that MAX_A will be convenient value)
2183  */
2184 #define M_LIMBS 31
2185 #define M_SIZE 122
2186
2187 static const t_uint limbs_M[] = { /* Little endian */
2188   0x84EEB59E, 0x9344A6AB, 0xFF21529F, 0xEC855CDA,
2189   0x009BAB38, 0x477E991E, 0x9F5B86F3, 0x2EEA2357,
2190   0x41D6502F, 0xAC17D304, 0x0A468A6D, 0x38FF52B9,
2191   0xFD42E5EF, 0x63630419, 0x91DB2572, 0x48CE17D0,
2192   0xE3B57D0E, 0x708AB00A, 0xCD723598, 0xF8A9DE08,
2193   0x4432C93B, 0x73141137, 0x2779FAB3, 0x554DF261,
2194   0x953D2BA5, 0xDEEBDA58, 0x5F57D007, 0xD1D66F2F,
2195   0xE84E9F2B, 0xB85C9607, 0x0000401D
2196 };
2197
2198 static const mpi M[1] = {{ 1, M_LIMBS, (t_uint *)limbs_M }};
2199
2200 /*
2201  * MAX_A : 2^1024 / M - 1
2202  */
2203 #define MAX_A_LIMBS 2
2204 #define MAX_A_FILL_SIZE  6
2205 static const t_uint limbs_MAX_A[] = { /* Little endian */
2206   0x56A2B35F, 0x0003FE25
2207 };
2208
2209 static const mpi MAX_A[1] = {{ 1, MAX_A_LIMBS, (t_uint *)limbs_MAX_A }};
2210
2211 /*
2212  * Prime number generation
2213  *
2214  * Special version for 1024-bit only.  Ignores DH_FLAG.
2215  */
2216 int mpi_gen_prime( mpi *X, size_t nbits, int dh_flag,
2217                    int (*f_rng)(void *, unsigned char *, size_t),
2218                    void *p_rng )
2219 {
2220   int ret;
2221   mpi B[1], G[1];
2222
2223   (void)dh_flag;
2224   if (nbits != 1024)
2225     return POLARSSL_ERR_MPI_BAD_INPUT_DATA;
2226
2227   mpi_init ( B );  mpi_init ( G );
2228
2229   /*
2230    * Get random value 1 to M-1 avoiding bias, and proceed when it is
2231    * coprime to all small primes.
2232    */
2233   do
2234     {
2235       MPI_CHK ( mpi_fill_random ( B, M_SIZE, f_rng, p_rng ) );
2236       B->p[0] |= 0x1;
2237       B->p[M_LIMBS - 1] &= 0x00007FFF;
2238       if (mpi_cmp_abs (B, M) >= 0)
2239         continue;
2240
2241       MPI_CHK ( mpi_gcd ( G, B, M ) );
2242     }
2243   while (mpi_cmp_int ( G, 1 ) != 0);
2244
2245   /*
2246    * Get random value avoiding bias, comput P with the value,
2247    * check if it's big enough, lastly, check if it's prime.
2248    */
2249   while (1)
2250     {
2251       MPI_CHK( mpi_fill_random( X, MAX_A_FILL_SIZE, f_rng, p_rng ) );
2252       MPI_CHK ( mpi_sub_abs (X, MAX_A, X) );
2253
2254       MPI_CHK ( mpi_mul_mpi ( X, X, M ) );
2255       MPI_CHK ( mpi_add_abs ( X, X, B ) );
2256       if (X->n <= 31 || (X->p[31] & 0xc0000000) == 0)
2257         continue;
2258
2259       ret = mpi_is_prime ( X );
2260       if (ret == 0 || ret != POLARSSL_ERR_MPI_NOT_ACCEPTABLE)
2261         break;
2262     }
2263
2264 cleanup:
2265
2266   mpi_free ( B );  mpi_free ( G );
2267
2268   return ret;
2269 }
2270
2271 #endif
2272
2273 #if defined(POLARSSL_SELF_TEST)
2274
2275 #define GCD_PAIR_COUNT  3
2276
2277 static const int gcd_pairs[GCD_PAIR_COUNT][3] =
2278 {
2279     { 693, 609, 21 },
2280     { 1764, 868, 28 },
2281     { 768454923, 542167814, 1 }
2282 };
2283
2284 /*
2285  * Checkup routine
2286  */
2287 int mpi_self_test( int verbose )
2288 {
2289     int ret, i;
2290     mpi A, E, N, X, Y, U, V;
2291
2292     mpi_init( &A ); mpi_init( &E ); mpi_init( &N ); mpi_init( &X );
2293     mpi_init( &Y ); mpi_init( &U ); mpi_init( &V );
2294
2295     MPI_CHK( mpi_read_string( &A, 16,
2296         "EFE021C2645FD1DC586E69184AF4A31E" \
2297         "D5F53E93B5F123FA41680867BA110131" \
2298         "944FE7952E2517337780CB0DB80E61AA" \
2299         "E7C8DDC6C5C6AADEB34EB38A2F40D5E6" ) );
2300
2301     MPI_CHK( mpi_read_string( &E, 16,
2302         "B2E7EFD37075B9F03FF989C7C5051C20" \
2303         "34D2A323810251127E7BF8625A4F49A5" \
2304         "F3E27F4DA8BD59C47D6DAABA4C8127BD" \
2305         "5B5C25763222FEFCCFC38B832366C29E" ) );
2306
2307     MPI_CHK( mpi_read_string( &N, 16,
2308         "0066A198186C18C10B2F5ED9B522752A" \
2309         "9830B69916E535C8F047518A889A43A5" \
2310         "94B6BED27A168D31D4A52F88925AA8F5" ) );
2311
2312     MPI_CHK( mpi_mul_mpi( &X, &A, &N ) );
2313
2314     MPI_CHK( mpi_read_string( &U, 16,
2315         "602AB7ECA597A3D6B56FF9829A5E8B85" \
2316         "9E857EA95A03512E2BAE7391688D264A" \
2317         "A5663B0341DB9CCFD2C4C5F421FEC814" \
2318         "8001B72E848A38CAE1C65F78E56ABDEF" \
2319         "E12D3C039B8A02D6BE593F0BBBDA56F1" \
2320         "ECF677152EF804370C1A305CAF3B5BF1" \
2321         "30879B56C61DE584A0F53A2447A51E" ) );
2322
2323     if( verbose != 0 )
2324         printf( "  MPI test #1 (mul_mpi): " );
2325
2326     if( mpi_cmp_mpi( &X, &U ) != 0 )
2327     {
2328         if( verbose != 0 )
2329             printf( "failed\n" );
2330
2331         return( 1 );
2332     }
2333
2334     if( verbose != 0 )
2335         printf( "passed\n" );
2336
2337     MPI_CHK( mpi_div_mpi( &X, &Y, &A, &N ) );
2338
2339     MPI_CHK( mpi_read_string( &U, 16,
2340         "256567336059E52CAE22925474705F39A94" ) );
2341
2342     MPI_CHK( mpi_read_string( &V, 16,
2343         "6613F26162223DF488E9CD48CC132C7A" \
2344         "0AC93C701B001B092E4E5B9F73BCD27B" \
2345         "9EE50D0657C77F374E903CDFA4C642" ) );
2346
2347     if( verbose != 0 )
2348         printf( "  MPI test #2 (div_mpi): " );
2349
2350     if( mpi_cmp_mpi( &X, &U ) != 0 ||
2351         mpi_cmp_mpi( &Y, &V ) != 0 )
2352     {
2353         if( verbose != 0 )
2354             printf( "failed\n" );
2355
2356         return( 1 );
2357     }
2358
2359     if( verbose != 0 )
2360         printf( "passed\n" );
2361
2362     MPI_CHK( mpi_exp_mod( &X, &A, &E, &N, NULL ) );
2363
2364     MPI_CHK( mpi_read_string( &U, 16,
2365         "36E139AEA55215609D2816998ED020BB" \
2366         "BD96C37890F65171D948E9BC7CBAA4D9" \
2367         "325D24D6A3C12710F10A09FA08AB87" ) );
2368
2369     if( verbose != 0 )
2370         printf( "  MPI test #3 (exp_mod): " );
2371
2372     if( mpi_cmp_mpi( &X, &U ) != 0 )
2373     {
2374         if( verbose != 0 )
2375             printf( "failed\n" );
2376
2377         return( 1 );
2378     }
2379
2380     if( verbose != 0 )
2381         printf( "passed\n" );
2382
2383 #if defined(POLARSSL_GENPRIME)
2384     MPI_CHK( mpi_inv_mod( &X, &A, &N ) );
2385
2386     MPI_CHK( mpi_read_string( &U, 16,
2387         "003A0AAEDD7E784FC07D8F9EC6E3BFD5" \
2388         "C3DBA76456363A10869622EAC2DD84EC" \
2389         "C5B8A74DAC4D09E03B5E0BE779F2DF61" ) );
2390
2391     if( verbose != 0 )
2392         printf( "  MPI test #4 (inv_mod): " );
2393
2394     if( mpi_cmp_mpi( &X, &U ) != 0 )
2395     {
2396         if( verbose != 0 )
2397             printf( "failed\n" );
2398
2399         return( 1 );
2400     }
2401
2402     if( verbose != 0 )
2403         printf( "passed\n" );
2404 #endif
2405
2406     if( verbose != 0 )
2407         printf( "  MPI test #5 (simple gcd): " );
2408
2409     for ( i = 0; i < GCD_PAIR_COUNT; i++)
2410     {
2411         MPI_CHK( mpi_lset( &X, gcd_pairs[i][0] ) );
2412         MPI_CHK( mpi_lset( &Y, gcd_pairs[i][1] ) );
2413
2414             MPI_CHK( mpi_gcd( &A, &X, &Y ) );
2415
2416             if( mpi_cmp_int( &A, gcd_pairs[i][2] ) != 0 )
2417             {
2418                     if( verbose != 0 )
2419                             printf( "failed at %d\n", i );
2420
2421                     return( 1 );
2422             }
2423     }
2424
2425     if( verbose != 0 )
2426         printf( "passed\n" );
2427
2428 cleanup:
2429
2430     if( ret != 0 && verbose != 0 )
2431         printf( "Unexpected error, return code = %08X\n", ret );
2432
2433     mpi_free( &A ); mpi_free( &E ); mpi_free( &N ); mpi_free( &X );
2434     mpi_free( &Y ); mpi_free( &U ); mpi_free( &V );
2435
2436     if( verbose != 0 )
2437         printf( "\n" );
2438
2439     return( ret );
2440 }
2441
2442 #endif
2443
2444 #endif