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