src/src/pdkim/bignum.c

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