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