src/src/pdkim/bignum.c

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