bignum.c 44.1 KB
Newer Older
1 2 3
/*
 *  Multi-precision integer library
 *
4
 *  Copyright (C) 2006-2010, Brainspark B.V.
Paul Bakker's avatar
Paul Bakker committed
5 6
 *
 *  This file is part of PolarSSL (http://www.polarssl.org)
7
 *  Lead Maintainer: Paul Bakker <polarssl_maintainer at polarssl.org>
Paul Bakker's avatar
Paul Bakker committed
8
 *
9
 *  All rights reserved.
Paul Bakker's avatar
Paul Bakker committed
10
 *
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along
 *  with this program; if not, write to the Free Software Foundation, Inc.,
 *  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */
/*
 *  This MPI implementation is based on:
 *
 *  http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
 *  http://www.stillhq.com/extracted/gnupg-api/mpi/
 *  http://math.libtomcrypt.com/files/tommath.pdf
 */

33
#include "polarssl/config.h"
34

35
#if defined(POLARSSL_BIGNUM_C)
36

37 38
#include "polarssl/bignum.h"
#include "polarssl/bn_mul.h"
39

40 41 42 43 44 45 46
#if defined(POLARSSL_MEMORY_C)
#include "polarssl/memory.h"
#else
#define polarssl_malloc     malloc
#define polarssl_free       free
#endif

47 48
#include <stdlib.h>

49
#define ciL    (sizeof(t_uint))         /* chars in limb  */
50 51 52 53 54 55 56 57 58 59
#define biL    (ciL << 3)               /* bits  in limb  */
#define biH    (ciL << 2)               /* half limb size */

/*
 * Convert between bits/chars and number of limbs
 */
#define BITS_TO_LIMBS(i)  (((i) + biL - 1) / biL)
#define CHARS_TO_LIMBS(i) (((i) + ciL - 1) / ciL)

/*
60
 * Initialize one MPI
61
 */
62
void mpi_init( mpi *X )
63
{
64 65
    if( X == NULL )
        return;
66

67 68 69
    X->s = 1;
    X->n = 0;
    X->p = NULL;
70 71 72
}

/*
73
 * Unallocate one MPI
74
 */
75
void mpi_free( mpi *X )
76
{
77 78
    if( X == NULL )
        return;
79

80
    if( X->p != NULL )
81
    {
82
        memset( X->p, 0, X->n * ciL );
83
        polarssl_free( X->p );
84 85
    }

86 87 88
    X->s = 1;
    X->n = 0;
    X->p = NULL;
89 90 91 92 93
}

/*
 * Enlarge to the specified number of limbs
 */
94
int mpi_grow( mpi *X, size_t nblimbs )
95
{
96
    t_uint *p;
97

98
    if( nblimbs > POLARSSL_MPI_MAX_LIMBS )
99
        return( POLARSSL_ERR_MPI_MALLOC_FAILED );
100

101 102
    if( X->n < nblimbs )
    {
103
        if( ( p = (t_uint *) polarssl_malloc( nblimbs * ciL ) ) == NULL )
104
            return( POLARSSL_ERR_MPI_MALLOC_FAILED );
105 106 107 108 109 110 111

        memset( p, 0, nblimbs * ciL );

        if( X->p != NULL )
        {
            memcpy( p, X->p, X->n * ciL );
            memset( X->p, 0, X->n * ciL );
112
            polarssl_free( X->p );
113 114 115 116 117 118 119 120 121 122 123 124
        }

        X->n = nblimbs;
        X->p = p;
    }

    return( 0 );
}

/*
 * Copy the contents of Y into X
 */
125
int mpi_copy( mpi *X, const mpi *Y )
126
{
127 128
    int ret;
    size_t i;
129 130 131 132

    if( X == Y )
        return( 0 );

133 134 135 136 137 138
    if( Y->p == NULL )
    {
        mpi_free( X );
        return( 0 );
    }

139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
    for( i = Y->n - 1; i > 0; i-- )
        if( Y->p[i] != 0 )
            break;
    i++;

    X->s = Y->s;

    MPI_CHK( mpi_grow( X, i ) );

    memset( X->p, 0, X->n * ciL );
    memcpy( X->p, Y->p, i * ciL );

cleanup:

    return( ret );
}

/*
 * Swap the contents of X and Y
 */
void mpi_swap( mpi *X, mpi *Y )
{
    mpi T;

    memcpy( &T,  X, sizeof( mpi ) );
    memcpy(  X,  Y, sizeof( mpi ) );
    memcpy(  Y, &T, sizeof( mpi ) );
}

/*
 * Set value from integer
 */
171
int mpi_lset( mpi *X, t_sint z )
172 173 174 175 176 177 178 179 180 181 182 183 184 185
{
    int ret;

    MPI_CHK( mpi_grow( X, 1 ) );
    memset( X->p, 0, X->n * ciL );

    X->p[0] = ( z < 0 ) ? -z : z;
    X->s    = ( z < 0 ) ? -1 : 1;

cleanup:

    return( ret );
}

186 187 188
/*
 * Get a specific bit
 */
189
int mpi_get_bit( const mpi *X, size_t pos )
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
{
    if( X->n * biL <= pos )
        return( 0 );

    return ( X->p[pos / biL] >> ( pos % biL ) ) & 0x01;
}

/*
 * Set a bit to a specific value of 0 or 1
 */
int mpi_set_bit( mpi *X, size_t pos, unsigned char val )
{
    int ret = 0;
    size_t off = pos / biL;
    size_t idx = pos % biL;

    if( val != 0 && val != 1 )
        return POLARSSL_ERR_MPI_BAD_INPUT_DATA;
        
    if( X->n * biL <= pos )
    {
        if( val == 0 )
            return ( 0 );

        MPI_CHK( mpi_grow( X, off + 1 ) );
    }

    X->p[off] = ( X->p[off] & ~( 0x01 << idx ) ) | ( val << idx );

cleanup:
    
    return( ret );
}

224 225 226
/*
 * Return the number of least significant bits
 */
227
size_t mpi_lsb( const mpi *X )
228
{
229
    size_t i, j, count = 0;
230 231

    for( i = 0; i < X->n; i++ )
232
        for( j = 0; j < biL; j++, count++ )
233 234 235 236 237 238 239 240 241
            if( ( ( X->p[i] >> j ) & 1 ) != 0 )
                return( count );

    return( 0 );
}

/*
 * Return the number of most significant bits
 */
242
size_t mpi_msb( const mpi *X )
243
{
244
    size_t i, j;
245 246 247 248 249

    for( i = X->n - 1; i > 0; i-- )
        if( X->p[i] != 0 )
            break;

250 251
    for( j = biL; j > 0; j-- )
        if( ( ( X->p[i] >> ( j - 1 ) ) & 1 ) != 0 )
252 253
            break;

254
    return( ( i * biL ) + j );
255 256 257 258 259
}

/*
 * Return the total size in bytes
 */
260
size_t mpi_size( const mpi *X )
261 262 263 264 265 266 267
{
    return( ( mpi_msb( X ) + 7 ) >> 3 );
}

/*
 * Convert an ASCII character to digit value
 */
268
static int mpi_get_digit( t_uint *d, int radix, char c )
269 270 271 272 273 274 275
{
    *d = 255;

    if( c >= 0x30 && c <= 0x39 ) *d = c - 0x30;
    if( c >= 0x41 && c <= 0x46 ) *d = c - 0x37;
    if( c >= 0x61 && c <= 0x66 ) *d = c - 0x57;

276
    if( *d >= (t_uint) radix )
277
        return( POLARSSL_ERR_MPI_INVALID_CHARACTER );
278 279 280 281 282 283 284

    return( 0 );
}

/*
 * Import from an ASCII string
 */
285
int mpi_read_string( mpi *X, int radix, const char *s )
286
{
287 288
    int ret;
    size_t i, j, slen, n;
289
    t_uint d;
290 291 292
    mpi T;

    if( radix < 2 || radix > 16 )
293
        return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
294

295
    mpi_init( &T );
296

297 298
    slen = strlen( s );

299 300
    if( radix == 16 )
    {
301
        n = BITS_TO_LIMBS( slen << 2 );
302 303 304 305

        MPI_CHK( mpi_grow( X, n ) );
        MPI_CHK( mpi_lset( X, 0 ) );

306
        for( i = slen, j = 0; i > 0; i--, j++ )
307
        {
308
            if( i == 1 && s[i - 1] == '-' )
309 310 311 312 313
            {
                X->s = -1;
                break;
            }

314
            MPI_CHK( mpi_get_digit( &d, radix, s[i - 1] ) );
315 316 317 318 319 320 321
            X->p[j / (2 * ciL)] |= d << ( (j % (2 * ciL)) << 2 );
        }
    }
    else
    {
        MPI_CHK( mpi_lset( X, 0 ) );

322
        for( i = 0; i < slen; i++ )
323 324 325 326 327 328 329 330 331
        {
            if( i == 0 && s[i] == '-' )
            {
                X->s = -1;
                continue;
            }

            MPI_CHK( mpi_get_digit( &d, radix, s[i] ) );
            MPI_CHK( mpi_mul_int( &T, X, radix ) );
332 333 334 335 336 337 338 339 340

            if( X->s == 1 )
            {
                MPI_CHK( mpi_add_int( X, &T, d ) );
            }
            else
            {
                MPI_CHK( mpi_sub_int( X, &T, d ) );
            }
341 342 343 344 345
        }
    }

cleanup:

346
    mpi_free( &T );
347 348 349 350 351 352 353 354 355 356

    return( ret );
}

/*
 * Helper to write the digits high-order first
 */
static int mpi_write_hlp( mpi *X, int radix, char **p )
{
    int ret;
357
    t_uint r;
358 359

    if( radix < 2 || radix > 16 )
360
        return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380

    MPI_CHK( mpi_mod_int( &r, X, radix ) );
    MPI_CHK( mpi_div_int( X, NULL, X, radix ) );

    if( mpi_cmp_int( X, 0 ) != 0 )
        MPI_CHK( mpi_write_hlp( X, radix, p ) );

    if( r < 10 )
        *(*p)++ = (char)( r + 0x30 );
    else
        *(*p)++ = (char)( r + 0x37 );

cleanup:

    return( ret );
}

/*
 * Export into an ASCII string
 */
381
int mpi_write_string( const mpi *X, int radix, char *s, size_t *slen )
382
{
383 384
    int ret = 0;
    size_t n;
385 386 387 388
    char *p;
    mpi T;

    if( radix < 2 || radix > 16 )
389
        return( POLARSSL_ERR_MPI_BAD_INPUT_DATA );
390 391 392 393 394 395 396 397 398

    n = mpi_msb( X );
    if( radix >=  4 ) n >>= 1;
    if( radix >= 16 ) n >>= 1;
    n += 3;

    if( *slen < n )
    {
        *slen = n;
399
        return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );
400 401 402
    }

    p = s;
403
    mpi_init( &T );
404 405 406 407 408 409

    if( X->s == -1 )
        *p++ = '-';

    if( radix == 16 )
    {
410 411
        int c;
        size_t i, j, k;
412

413
        for( i = X->n, k = 0; i > 0; i-- )
414
        {
415
            for( j = ciL; j > 0; j-- )
416
            {
417
                c = ( X->p[i - 1] >> ( ( j - 1 ) << 3) ) & 0xFF;
418

419
                if( c == 0 && k == 0 && ( i + j + 3 ) != 0 )
420 421
                    continue;

422
                *(p++) = "0123456789ABCDEF" [c / 16];
Paul Bakker's avatar
Paul Bakker committed
423
                *(p++) = "0123456789ABCDEF" [c % 16];
424 425 426 427 428 429 430
                k = 1;
            }
        }
    }
    else
    {
        MPI_CHK( mpi_copy( &T, X ) );
431 432 433 434

        if( T.s == -1 )
            T.s = 1;

435 436 437 438 439 440 441 442
        MPI_CHK( mpi_write_hlp( &T, radix, &p ) );
    }

    *p++ = '\0';
    *slen = p - s;

cleanup:

443
    mpi_free( &T );
444 445 446 447

    return( ret );
}

448
#if defined(POLARSSL_FS_IO)
449 450 451 452 453
/*
 * Read X from an opened file
 */
int mpi_read_file( mpi *X, int radix, FILE *fin )
{
454
    t_uint d;
455
    size_t slen;
456
    char *p;
457
    /*
458 459
     * Buffer should have space for (short) label and decimal formatted MPI,
     * newline characters and '\0'
460
     */
461
    char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
462 463 464

    memset( s, 0, sizeof( s ) );
    if( fgets( s, sizeof( s ) - 1, fin ) == NULL )
465
        return( POLARSSL_ERR_MPI_FILE_IO_ERROR );
466 467

    slen = strlen( s );
468 469 470
    if( slen == sizeof( s ) - 2 )
        return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );

471 472 473 474 475 476 477 478 479 480 481 482 483 484
    if( s[slen - 1] == '\n' ) { slen--; s[slen] = '\0'; }
    if( s[slen - 1] == '\r' ) { slen--; s[slen] = '\0'; }

    p = s + slen;
    while( --p >= s )
        if( mpi_get_digit( &d, radix, *p ) != 0 )
            break;

    return( mpi_read_string( X, radix, p + 1 ) );
}

/*
 * Write X into an opened file (or stdout if fout == NULL)
 */
485
int mpi_write_file( const char *p, const mpi *X, int radix, FILE *fout )
486
{
487 488
    int ret;
    size_t n, slen, plen;
489
    /*
490 491
     * Buffer should have space for (short) label and decimal formatted MPI,
     * newline characters and '\0'
492
     */
493
    char s[ POLARSSL_MPI_RW_BUFFER_SIZE ];
494 495 496 497 498

    n = sizeof( s );
    memset( s, 0, n );
    n -= 2;

499
    MPI_CHK( mpi_write_string( X, radix, s, (size_t *) &n ) );
500 501 502 503 504 505 506 507 508 509 510 511

    if( p == NULL ) p = "";

    plen = strlen( p );
    slen = strlen( s );
    s[slen++] = '\r';
    s[slen++] = '\n';

    if( fout != NULL )
    {
        if( fwrite( p, 1, plen, fout ) != plen ||
            fwrite( s, 1, slen, fout ) != slen )
512
            return( POLARSSL_ERR_MPI_FILE_IO_ERROR );
513 514 515 516 517 518 519 520
    }
    else
        printf( "%s%s", p, s );

cleanup:

    return( ret );
}
521
#endif /* POLARSSL_FS_IO */
522 523 524 525

/*
 * Import X from unsigned binary data, big endian
 */
526
int mpi_read_binary( mpi *X, const unsigned char *buf, size_t buflen )
527
{
528 529
    int ret;
    size_t i, j, n;
530 531 532 533 534 535 536 537

    for( n = 0; n < buflen; n++ )
        if( buf[n] != 0 )
            break;

    MPI_CHK( mpi_grow( X, CHARS_TO_LIMBS( buflen - n ) ) );
    MPI_CHK( mpi_lset( X, 0 ) );

538
    for( i = buflen, j = 0; i > n; i--, j++ )
539
        X->p[j / ciL] |= ((t_uint) buf[i - 1]) << ((j % ciL) << 3);
540 541 542 543 544 545 546 547 548

cleanup:

    return( ret );
}

/*
 * Export X into unsigned binary data, big endian
 */
549
int mpi_write_binary( const mpi *X, unsigned char *buf, size_t buflen )
550
{
551
    size_t i, j, n;
552 553 554 555

    n = mpi_size( X );

    if( buflen < n )
556
        return( POLARSSL_ERR_MPI_BUFFER_TOO_SMALL );
557 558 559 560 561 562 563 564 565 566 567 568

    memset( buf, 0, buflen );

    for( i = buflen - 1, j = 0; n > 0; i--, j++, n-- )
        buf[i] = (unsigned char)( X->p[j / ciL] >> ((j % ciL) << 3) );

    return( 0 );
}

/*
 * Left-shift: X <<= count
 */
569
int mpi_shift_l( mpi *X, size_t count )
570
{
571 572
    int ret;
    size_t i, v0, t1;
573
    t_uint r0 = 0, r1;
574 575 576 577 578 579

    v0 = count / (biL    );
    t1 = count & (biL - 1);

    i = mpi_msb( X ) + count;

580
    if( X->n * biL < i )
581 582 583 584 585 586 587 588 589
        MPI_CHK( mpi_grow( X, BITS_TO_LIMBS( i ) ) );

    ret = 0;

    /*
     * shift by count / limb_size
     */
    if( v0 > 0 )
    {
590 591
        for( i = X->n; i > v0; i-- )
            X->p[i - 1] = X->p[i - v0 - 1];
592

593 594
        for( ; i > 0; i-- )
            X->p[i - 1] = 0;
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618
    }

    /*
     * shift by count % limb_size
     */
    if( t1 > 0 )
    {
        for( i = v0; i < X->n; i++ )
        {
            r1 = X->p[i] >> (biL - t1);
            X->p[i] <<= t1;
            X->p[i] |= r0;
            r0 = r1;
        }
    }

cleanup:

    return( ret );
}

/*
 * Right-shift: X >>= count
 */
619
int mpi_shift_r( mpi *X, size_t count )
620
{
621
    size_t i, v0, v1;
622
    t_uint r0 = 0, r1;
623 624 625 626

    v0 = count /  biL;
    v1 = count & (biL - 1);

627 628 629
    if( v0 > X->n || ( v0 == X->n && v1 > 0 ) )
        return mpi_lset( X, 0 );

630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646
    /*
     * shift by count / limb_size
     */
    if( v0 > 0 )
    {
        for( i = 0; i < X->n - v0; i++ )
            X->p[i] = X->p[i + v0];

        for( ; i < X->n; i++ )
            X->p[i] = 0;
    }

    /*
     * shift by count % limb_size
     */
    if( v1 > 0 )
    {
647
        for( i = X->n; i > 0; i-- )
648
        {
649 650 651
            r1 = X->p[i - 1] << (biL - v1);
            X->p[i - 1] >>= v1;
            X->p[i - 1] |= r0;
652 653 654 655 656 657 658 659 660 661
            r0 = r1;
        }
    }

    return( 0 );
}

/*
 * Compare unsigned values
 */
662
int mpi_cmp_abs( const mpi *X, const mpi *Y )
663
{
664
    size_t i, j;
665

666 667
    for( i = X->n; i > 0; i-- )
        if( X->p[i - 1] != 0 )
668 669
            break;

670 671
    for( j = Y->n; j > 0; j-- )
        if( Y->p[j - 1] != 0 )
672 673
            break;

674
    if( i == 0 && j == 0 )
675 676 677 678 679
        return( 0 );

    if( i > j ) return(  1 );
    if( j > i ) return( -1 );

680
    for( ; i > 0; i-- )
681
    {
682 683
        if( X->p[i - 1] > Y->p[i - 1] ) return(  1 );
        if( X->p[i - 1] < Y->p[i - 1] ) return( -1 );
684 685 686 687 688 689 690 691
    }

    return( 0 );
}

/*
 * Compare signed values
 */
692
int mpi_cmp_mpi( const mpi *X, const mpi *Y )
693
{
694
    size_t i, j;
695

696 697
    for( i = X->n; i > 0; i-- )
        if( X->p[i - 1] != 0 )
698 699
            break;

700 701
    for( j = Y->n; j > 0; j-- )
        if( Y->p[j - 1] != 0 )
702 703
            break;

704
    if( i == 0 && j == 0 )
705 706 707
        return( 0 );

    if( i > j ) return(  X->s );
708
    if( j > i ) return( -Y->s );
709 710 711 712

    if( X->s > 0 && Y->s < 0 ) return(  1 );
    if( Y->s > 0 && X->s < 0 ) return( -1 );

713
    for( ; i > 0; i-- )
714
    {
715 716
        if( X->p[i - 1] > Y->p[i - 1] ) return(  X->s );
        if( X->p[i - 1] < Y->p[i - 1] ) return( -X->s );
717 718 719 720 721 722 723 724
    }

    return( 0 );
}

/*
 * Compare signed values
 */
725
int mpi_cmp_int( const mpi *X, t_sint z )
726 727
{
    mpi Y;
728
    t_uint p[1];
729 730 731 732 733 734 735 736 737 738 739 740

    *p  = ( z < 0 ) ? -z : z;
    Y.s = ( z < 0 ) ? -1 : 1;
    Y.n = 1;
    Y.p = p;

    return( mpi_cmp_mpi( X, &Y ) );
}

/*
 * Unsigned addition: X = |A| + |B|  (HAC 14.7)
 */
741
int mpi_add_abs( mpi *X, const mpi *A, const mpi *B )
742
{
743 744
    int ret;
    size_t i, j;
745
    t_uint *o, *p, c;
746 747 748

    if( X == B )
    {
749
        const mpi *T = A; A = X; B = T;
750 751 752 753
    }

    if( X != A )
        MPI_CHK( mpi_copy( X, A ) );
754 755 756 757 758
   
    /*
     * X should always be positive as a result of unsigned additions.
     */
    X->s = 1;
759

760 761
    for( j = B->n; j > 0; j-- )
        if( B->p[j - 1] != 0 )
762 763
            break;

764
    MPI_CHK( mpi_grow( X, j ) );
765 766 767

    o = B->p; p = X->p; c = 0;

768
    for( i = 0; i < j; i++, o++, p++ )
769 770 771 772 773 774 775 776 777 778 779 780 781
    {
        *p +=  c; c  = ( *p <  c );
        *p += *o; c += ( *p < *o );
    }

    while( c != 0 )
    {
        if( i >= X->n )
        {
            MPI_CHK( mpi_grow( X, i + 1 ) );
            p = X->p + i;
        }

782
        *p += c; c = ( *p < c ); i++; p++;
783 784 785 786 787 788 789 790 791 792
    }

cleanup:

    return( ret );
}

/*
 * Helper for mpi substraction
 */
793
static void mpi_sub_hlp( size_t n, t_uint *s, t_uint *d )
794
{
795
    size_t i;
796
    t_uint c, z;
797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813

    for( i = c = 0; i < n; i++, s++, d++ )
    {
        z = ( *d <  c );     *d -=  c;
        c = ( *d < *s ) + z; *d -= *s;
    }

    while( c != 0 )
    {
        z = ( *d < c ); *d -= c;
        c = z; i++; d++;
    }
}

/*
 * Unsigned substraction: X = |A| - |B|  (HAC 14.9)
 */
814
int mpi_sub_abs( mpi *X, const mpi *A, const mpi *B )
815 816
{
    mpi TB;
817 818
    int ret;
    size_t n;
819 820

    if( mpi_cmp_abs( A, B ) < 0 )
821
        return( POLARSSL_ERR_MPI_NEGATIVE_VALUE );
822

823
    mpi_init( &TB );
824 825 826 827 828 829 830 831 832 833

    if( X == B )
    {
        MPI_CHK( mpi_copy( &TB, B ) );
        B = &TB;
    }

    if( X != A )
        MPI_CHK( mpi_copy( X, A ) );

834 835 836 837 838
    /*
     * X should always be positive as a result of unsigned substractions.
     */
    X->s = 1;

839 840
    ret = 0;

841 842
    for( n = B->n; n > 0; n-- )
        if( B->p[n - 1] != 0 )
843 844
            break;

845
    mpi_sub_hlp( n, B->p, X->p );
846 847 848

cleanup:

849
    mpi_free( &TB );
850 851 852 853 854 855 856

    return( ret );
}

/*
 * Signed addition: X = A + B
 */
857
int mpi_add_mpi( mpi *X, const mpi *A, const mpi *B )
858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887
{
    int ret, s = A->s;

    if( A->s * B->s < 0 )
    {
        if( mpi_cmp_abs( A, B ) >= 0 )
        {
            MPI_CHK( mpi_sub_abs( X, A, B ) );
            X->s =  s;
        }
        else
        {
            MPI_CHK( mpi_sub_abs( X, B, A ) );
            X->s = -s;
        }
    }
    else
    {
        MPI_CHK( mpi_add_abs( X, A, B ) );
        X->s = s;
    }

cleanup:

    return( ret );
}

/*
 * Signed substraction: X = A - B
 */
888
int mpi_sub_mpi( mpi *X, const mpi *A, const mpi *B )
889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918
{
    int ret, s = A->s;

    if( A->s * B->s > 0 )
    {
        if( mpi_cmp_abs( A, B ) >= 0 )
        {
            MPI_CHK( mpi_sub_abs( X, A, B ) );
            X->s =  s;
        }
        else
        {
            MPI_CHK( mpi_sub_abs( X, B, A ) );
            X->s = -s;
        }
    }
    else
    {
        MPI_CHK( mpi_add_abs( X, A, B ) );
        X->s = s;
    }

cleanup:

    return( ret );
}

/*
 * Signed addition: X = A + b
 */
919
int mpi_add_int( mpi *X, const mpi *A, t_sint b )
920 921
{
    mpi _B;
922
    t_uint p[1];
923 924 925 926 927 928 929 930 931 932 933 934

    p[0] = ( b < 0 ) ? -b : b;
    _B.s = ( b < 0 ) ? -1 : 1;
    _B.n = 1;
    _B.p = p;

    return( mpi_add_mpi( X, A, &_B ) );
}

/*
 * Signed substraction: X = A - b
 */
935
int mpi_sub_int( mpi *X, const mpi *A, t_sint b )
936 937
{
    mpi _B;
938
    t_uint p[1];
939 940 941 942 943 944 945 946 947 948 949

    p[0] = ( b < 0 ) ? -b : b;
    _B.s = ( b < 0 ) ? -1 : 1;
    _B.n = 1;
    _B.p = p;

    return( mpi_sub_mpi( X, A, &_B ) );
}

/*
 * Helper for mpi multiplication
950 951 952 953 954 955 956 957 958 959
 */
static
#if defined(__APPLE__) && defined(__arm__)
/*
 * Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
 * appears to need this to prevent bad ARM code generation at -O3.
 */
__attribute__ ((noinline))
#endif
void mpi_mul_hlp( size_t i, t_uint *s, t_uint *d, t_uint b )
960
{
961
    t_uint c = 0, t = 0;
962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022

#if defined(MULADDC_HUIT)
    for( ; i >= 8; i -= 8 )
    {
        MULADDC_INIT
        MULADDC_HUIT
        MULADDC_STOP
    }

    for( ; i > 0; i-- )
    {
        MULADDC_INIT
        MULADDC_CORE
        MULADDC_STOP
    }
#else
    for( ; i >= 16; i -= 16 )
    {
        MULADDC_INIT
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE

        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_STOP
    }

    for( ; i >= 8; i -= 8 )
    {
        MULADDC_INIT
        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE

        MULADDC_CORE   MULADDC_CORE
        MULADDC_CORE   MULADDC_CORE
        MULADDC_STOP
    }

    for( ; i > 0; i-- )
    {
        MULADDC_INIT
        MULADDC_CORE
        MULADDC_STOP
    }
#endif

    t++;

    do {
        *d += c; c = ( *d < c ); d++;
    }
    while( c != 0 );
}

/*
 * Baseline multiplication: X = A * B  (HAC 14.12)
 */
1023
int mpi_mul_mpi( mpi *X, const mpi *A, const mpi *B )
1024
{
1025 1026
    int ret;
    size_t i, j;
1027 1028
    mpi TA, TB;

1029
    mpi_init( &TA ); mpi_init( &TB );
1030 1031 1032 1033

    if( X == A ) { MPI_CHK( mpi_copy( &TA, A ) ); A = &TA; }
    if( X == B ) { MPI_CHK( mpi_copy( &TB, B ) ); B = &TB; }

1034 1035
    for( i = A->n; i > 0; i-- )
        if( A->p[i - 1] != 0 )
1036 1037
            break;

1038 1039
    for( j = B->n; j > 0; j-- )
        if( B->p[j - 1] != 0 )
1040 1041
            break;

1042
    MPI_CHK( mpi_grow( X, i + j ) );
1043 1044
    MPI_CHK( mpi_lset( X, 0 ) );

1045 1046
    for( i++; j > 0; j-- )
        mpi_mul_hlp( i - 1, A->p, X->p + j - 1, B->p[j - 1] );
1047 1048 1049 1050 1051

    X->s = A->s * B->s;

cleanup:

1052
    mpi_free( &TB ); mpi_free( &TA );
1053 1054 1055 1056 1057 1058 1059

    return( ret );
}

/*
 * Baseline multiplication: X = A * b
 */
1060
int mpi_mul_int( mpi *X, const mpi *A, t_sint b )
1061 1062
{
    mpi _B;
1063
    t_uint p[1];
1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075

    _B.s = 1;
    _B.n = 1;
    _B.p = p;
    p[0] = b;

    return( mpi_mul_mpi( X, A, &_B ) );
}

/*
 * Division by mpi: A = Q * B + R  (HAC 14.20)
 */
1076
int mpi_div_mpi( mpi *Q, mpi *R, const mpi *A, const mpi *B )
1077
{
1078 1079
    int ret;
    size_t i, n, t, k;
1080 1081 1082
    mpi X, Y, Z, T1, T2;

    if( mpi_cmp_int( B, 0 ) == 0 )
1083
        return( POLARSSL_ERR_MPI_DIVISION_BY_ZERO );
1084

1085 1086
    mpi_init( &X ); mpi_init( &Y ); mpi_init( &Z );
    mpi_init( &T1 ); mpi_init( &T2 );
1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104

    if( mpi_cmp_abs( A, B ) < 0 )
    {
        if( Q != NULL ) MPI_CHK( mpi_lset( Q, 0 ) );
        if( R != NULL ) MPI_CHK( mpi_copy( R, A ) );
        return( 0 );
    }

    MPI_CHK( mpi_copy( &X, A ) );
    MPI_CHK( mpi_copy( &Y, B ) );
    X.s = Y.s = 1;

    MPI_CHK( mpi_grow( &Z, A->n + 2 ) );
    MPI_CHK( mpi_lset( &Z,  0 ) );
    MPI_CHK( mpi_grow( &T1, 2 ) );
    MPI_CHK( mpi_grow( &T2, 3 ) );

    k = mpi_msb( &Y ) % biL;
1105
    if( k < biL - 1 )
1106 1107 1108 1109 1110 1111 1112 1113 1114
    {
        k = biL - 1 - k;
        MPI_CHK( mpi_shift_l( &X, k ) );
        MPI_CHK( mpi_shift_l( &Y, k ) );
    }
    else k = 0;

    n = X.n - 1;
    t = Y.n - 1;
1115
    MPI_CHK( mpi_shift_l( &Y, biL * (n - t) ) );
1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129

    while( mpi_cmp_mpi( &X, &Y ) >= 0 )
    {
        Z.p[n - t]++;
        mpi_sub_mpi( &X, &X, &Y );
    }
    mpi_shift_r( &Y, biL * (n - t) );

    for( i = n; i > t ; i-- )
    {
        if( X.p[i] >= Y.p[t] )
            Z.p[i - t - 1] = ~0;
        else
        {
1130
#if defined(POLARSSL_HAVE_UDBL)
1131
            t_udbl r;
1132

1133 1134
            r  = (t_udbl) X.p[i] << biL;
            r |= (t_udbl) X.p[i - 1];
1135
            r /= Y.p[t];
1136 1137
            if( r > ((t_udbl) 1 << biL) - 1)
                r = ((t_udbl) 1 << biL) - 1;
1138

1139
            Z.p[i - t - 1] = (t_uint) r;
1140 1141 1142 1143
#else
            /*
             * __udiv_qrnnd_c, from gmp/longlong.h
             */
1144 1145
            t_uint q0, q1, r0, r1;
            t_uint d0, d1, d, m