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