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