mpi.cpp

RSA加密解密的算法例程

 * Right-shift: X >>= count
 *
 * Always returns 0.
 */
int mpi_shift_r( mpi *X, uint count )
{
    int i, v0, v1;
    t_int r0 = 0, r1;

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

    /*
     * 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 )
    {
        for( i = X->n - 1; i >= 0; i-- )
        {
            r1 = X->p[i] << (biL - v1);
            X->p[i] >>= v1;
            X->p[i] |= r0;
            r0 = r1;
        }
    }

    return( 0 );
}

/*
 * Compare absolute values
 *
 * Returns 1 if |X| is greater than |Y|
 *        -1 if |X| is lesser  than |Y|
 *         0 if |X| is equal to |Y|
 */
int mpi_cmp_abs( mpi *X, mpi *Y )
{
    int i, j;

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

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

    if( i < 0 && j < 0 )
        return( 0 );

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

    for( ; i >= 0; i-- )
    {
        if( X->p[i] > Y->p[i] ) return(  1 );
        if( X->p[i] < Y->p[i] ) return( -1 );
    }

    return( 0 );
}

/*
 * Compare signed values
 *
 * Returns 1 if X is greater than Y
 *        -1 if X is lesser  than Y
 *         0 if X is equal to Y
 */
int mpi_cmp_mpi( mpi *X, mpi *Y )
{
    int i, j;

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

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

    if( i < 0 && j < 0 )
        return( 0 );

    if( i > j ) return(  X->s );
    if( j > i ) return( -X->s );

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

    for( ; i >= 0; i-- )
    {
        if( X->p[i] > Y->p[i] ) return(  X->s );
        if( X->p[i] < Y->p[i] ) return( -X->s );
    }

    return( 0 );
}

/*
 * Compare signed values
 *
 * Returns 1 if X is greater than z
 *        -1 if X is lesser  than z
 *         0 if X is equal to z
 */
int mpi_cmp_int( mpi *X, int z )
{
    mpi Y;
    t_int p[1];

    *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)
 *
 * Returns 0 if successful
 *         1 if memory allocation failed
 */
int mpi_add_abs( mpi *X, mpi *A, mpi *B )
{
    int ret, i, j;
    t_int *o, *p, c;

    if( X == B )
    {
        mpi *T = A; A = X; B = T;
    }

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

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

    CHK( mpi_grow( X, j + 1 ) );

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

    for( i = 0; i <= j; i++, o++, p++ )
    {
        *p +=  c; c  = ( *p <  c );
        *p += *o; c += ( *p < *o );
    }

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

        *p += c; c = ( *p < c ); i++;
    }

cleanup:

    return( ret );
}

/*
 * Unsigned substraction: X = |A| - |B|  (HAC 14.9)
 *
 * Returns 0 if successful
 *         ERR_MPI_NEGATIVE_VALUE if B is greater than A
 */
int mpi_sub_abs( mpi *X, mpi *A, mpi *B )
{
    mpi TB;
    int ret, i, j;
    t_int *o, *p, c, z;

    if( mpi_cmp_abs( A, B ) < 0 )
        return( ERR_MPI_NEGATIVE_VALUE );

    mpi_init( &TB, NULL );

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

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

    ret = 0;

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

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

    for( i = 0; i <= j; i++, o++, p++ )
    {
        z = ( *p <  c );     *p -=  c;
        c = ( *p < *o ) + z; *p -= *o;
    }

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

cleanup:

    mpi_free( &TB, NULL );
    return( ret );
}

/*
 * Signed addition: X = A + B
 *
 * Returns 0 if successful
 *         1 if memory allocation failed
 */
int mpi_add_mpi( mpi *X, mpi *A, mpi *B )
{
    int ret, s = A->s;

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

cleanup:

    return( ret );
}

/*
 * Signed substraction: X = A - B
 *
 * Returns 0 if successful
 *         1 if memory allocation failed
 */
int mpi_sub_mpi( mpi *X, mpi *A, mpi *B )
{
    int ret, s = A->s;

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

cleanup:

    return( ret );
}

/*
 * Signed addition: X = A + b
 *
 * Returns 0 if successful
 *         1 if memory allocation failed
 */
int mpi_add_int( mpi *X, mpi *A, int b )
{
    mpi _B;
    t_int p[1];

    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
 *
 * Returns 0 if successful,
 *         1 if memory allocation failed
 */
int mpi_sub_int( mpi *X, mpi *A, int b )
{
    mpi _B;
    t_int p[1];

    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 ) );
}

/* 
 * Multiply source s with b, add result to
 * destination d and set carry c.
 */
#if defined _MSC_VER

#define MULADDC_INIT                    \
    __asm   mov     esi, s              \
    __asm   mov     edi, d              \
    __asm   mov     ecx, c              \
    __asm   mov     ebx, b

#define MULADDC_CORE                    \
    __asm   lodsd                       \
    __asm   mul     ebx                 \
    __asm   add     eax, ecx            \
    __asm   adc     edx, 0              \
    __asm   add     eax, [edi]          \
    __asm   adc     edx, 0              \
    __asm   mov     ecx, edx            \
    __asm   stosd

#define MULADDC_STOP                    \
    __asm   mov     c, ecx              \
    __asm   mov     d, edi              \
    __asm   mov     s, esi              \

#else
#if defined __i386__

#define MULADDC_INIT                    \
    asm( "movl  %0, %%esi" :: "m" (s)); \
    asm( "movl  %0, %%edi" :: "m" (d)); \
    asm( "movl  %0, %%ecx" :: "m" (c)); \
    asm( "movl  %0, %%ebx" :: "m" (b));

#define MULADDC_CORE                    \
    asm( "lodsl                 " );    \
    asm( "mull  %ebx            " );    \
    asm( "addl  %ecx,   %eax    " );    \
    asm( "adcl  $0,     %edx    " );    \
    asm( "addl  (%edi), %eax    " );    \
    asm( "adcl  $0,     %edx    " );    \
    asm( "movl  %edx,   %ecx    " );    \
    asm( "stosl                 " );

#define MULADDC_STOP                    \
    asm( "movl  %%ecx, %0" :: "m" (c)); \
    asm( "movl  %%edi, %0" :: "m" (d)); \
    asm( "movl  %%esi, %0" :: "m" (s) : \
         "eax", "ecx", "edx",           \
         "ebx", "esi", "edi" );

#else
#if defined __amd64__

#define MULADDC_INIT                    \
    asm( "movq  %0, %%rsi" :: "m" (s)); \
    asm( "movq  %0, %%rdi" :: "m" (d)); \
    asm( "movq  %0, %%rcx" :: "m" (c)); \
    asm( "movq  %0, %%rbx" :: "m" (b));

#define MULADDC_CORE                    \
    asm( "lodsq                 " );    \
    asm( "mulq  %rbx            " );    \
    asm( "addq  %rcx,   %rax    " );    \
    asm( "adcq  $0,     %rdx    " );    \
    asm( "addq  (%rdi), %rax    " );    \
    asm( "adcq  $0,     %rdx    " );    \
    asm( "movq  %rdx,   %rcx    " );    \
    asm( "stosq                 " );

#define MULADDC_STOP                    \
    asm( "movq  %%rcx, %0" :: "m" (c)); \
    asm( "movq  %%rdi, %0" :: "m" (d)); \
    asm( "movq  %%rsi, %0" :: "m" (s) : \
         "rax", "rcx", "rdx",           \
         "rbx", "rsi", "rdi" );

#else
#warning no muladdc assembly code for this cpu

#define MULADDC_INIT                    \
{                                       \
    t_int s0, s1, b0, b1;               \
    t_int r0, r1, rx, ry;               \
    b0 = ( b << biH ) >> biH;           \
    b1 = ( b >> biH );

#define MULADDC_CORE                    \
    s0 = ( *s << biH ) >> biH;          \
    s1 = ( *s >> biH ); s++;            \
    rx = s0 * b1; r0 = s0 * b0;         \
    ry = s1 * b0; r1 = s1 * b1;         \
    r1 += ( rx >> biH );                \
    r1 += ( ry >> biH );                \
    rx <<= biH; ry <<= biH;             \
    r0 += rx; r1 += (r0 < rx);          \
    r0 += ry; r1 += (r0 < ry);          \
    r0 +=  c; r1 += (r0 <  c);          \
    r0 += *d; r1 += (r0 < *d);          \
    c = r1; *(d++) = r0;

#define MULADDC_STOP                    \
}

#endif
#endif
#endif

void MULADDC( int i, t_int *s, t_int *d, t_int b )
{
    t_int c = 0;

    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 >= 4; i -= 4 )
    {
        MULADDC_INIT
        MULADDC_CORE  MULADDC_CORE
        MULADDC_CORE  MULADDC_CORE
        MULADDC_STOP
    }

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

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

/*
 * Baseline multiplication: X = A * B  (HAC 14.12)
 *
 * Returns 0 if successful
 *         1 if memory allocation failed
 */
int mpi_mul_mpi( mpi *X, mpi *A, mpi *B )
{
    int ret, i, j;
    mpi TA, TB;