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