src/src/pdkim/bignum.c

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