mp.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 1,082 行 · 第 1/2 页

    }

    /* check if dividend is zero or x is less than y */
    if( mp_eq( &q, x ) || mp_lt( x, y ) ) {
        rc = MP_NO_ERR;
        goto done;
    }

    /* make sure n and t are true lengths */
    while( x->num[n-1] == 0 ) {
        n--;
        if( n == 0 ) goto done;
    }
    while( y->num[t-1] == 0 ) {
        t--;
        if( t == 0 ) goto done;
    }

    /* make sure first digit of divisor is >= base / 2 to ensure we get a
       good guess */
    shift = 31 - ((mp_bitsize( y ) - 1) % 32);
    if( ((mp_bitsize( &r ) + 1) % 32 - 1) + shift > 32 ) n++;
    mp_shiftleft( y, y, shift );
    mp_shiftleft( &r, &r, shift );

    /* step 2 */
    if( temp.allocated < n - t +1 ) MP_GROW( &temp, n - t - temp.allocated + 1 );
    temp.num[n-t] = 1;
    temp.len = n - t + 1;
    rc = mp_mul( &temp, &temp, y );
    if( rc != MP_NO_ERR ) goto done;
    for(;;) {
        if( mp_lt( &r, &temp ) ) break;
        q.num[n-t]++;
        mp_sub( &r, &r, &temp );
    }

    /* step 3 */
    for( i = n; i >= t + 1; i-- ) {
        /* 3.1 - make guess using two digits */
        if( r.num[i-1] == y->num[t-1] ) {
            q.num[i-t-1] = 0xFFFFFFFF;
        } else {
            tmp = ((uint64)(r.num[i-1])<<32  | (uint64)r.num[i-2]) / y->num[t-1];
            q.num[i-t-1] = (uint32) tmp;
        }
        /* 3.2 - improve guess using three digits */
        if( t > 1 ) {
            tmp = ((uint64)(r.num[i-1])<<32) | (uint64)r.num[i-2];
            MP_INIT( &mptmp1, tmp );
            MP_INIT( &mptmp3, 0 );
            mp_shiftleft( &mptmp1, &mptmp1, 32 );
            mp_addsc( &mptmp1, &mptmp1, r.num[i-3] );
            tmp = ((uint64)(y->num[t-1])<<32) | (uint64)y->num[t-2];
            MP_INIT( &mptmp2, tmp );
            mp_mulsc( &mptmp3, &mptmp2, q.num[i-t-1] );
            for(;;) {
                if( mp_lte( &mptmp3, &mptmp1 ) ) break;
                q.num[i-t-1]--;
                mp_sub( &mptmp3, &mptmp3, &mptmp2 );
            }
            mp_free( &mptmp1 );
            mp_free( &mptmp2 );
            mp_free( &mptmp3 );
        }
        /* 3.3 + 3.4 */
        mp_zero( &temp );
        mp_zero( &temp2 );
        temp.len = i - t;
        temp.num[i-t-1] = 1;
        rc = mp_mul( &temp, &temp, y );
        if( rc != MP_NO_ERR ) goto done;
        rc = mp_mulsc( &temp2, &temp, q.num[i-t-1] );
        if( rc != MP_NO_ERR ) goto done;
        if( mp_lt( &r, &temp2 ) ) {
            /* the guess was off by one */
            rc = mp_sub( &temp2, &temp2, &temp );
            if( rc != MP_NO_ERR ) goto done;
            q.num[i-t-1]--;
        }
        rc = mp_sub( &r, &r, &temp2 );
        if( rc != MP_NO_ERR ) goto done;
    }

    rc = MP_NO_ERR;
done:
    if( rc== MP_NO_ERR ) {
        mp_shiftright( y, y, shift );
        MP_COPY( qdst, &q );
        MP_COPY( rdst, &r );
        mp_shiftright( rdst, rdst, shift );
    }
    mp_free( &temp );
    mp_free( &temp2 );
    mp_free( &q );
    mp_free( &r );
    return rc;
}

/* sqaure src */
int mp_sqr( mpnum *dst, mpnum *src )
{
    return mp_mul( dst, src, src );
}

/* calculate src to the exp power using the binary method */
int mp_pow( mpnum *dst, mpnum *src, uint32 exp )
{
    mpnum temp;
    uint32 bit = 1;
    int rc;

    MP_INIT( &temp, 0 );
    MP_COPY( &temp, src );
    mp_free( dst );
    MP_INIT( dst, 1 );

    for( ;; ) {
        if( exp & bit ) {
            rc = mp_mul( dst, dst, &temp );
            if( rc != MP_NO_ERR ) goto done;
        }
        bit <<= 1;
        if( bit > exp ) break;
        rc = mp_sqr( &temp, &temp );
        if( rc != MP_NO_ERR ) goto done;
    }

    rc = MP_NO_ERR;

done:
    mp_free( &temp );

    return rc;
}

/* Do a binary rounding at the 'bit'th bit;
 * Round to nearest rounding mode is used and if there are two value equally
 * close to the number, it rounds to the even one. */
int mp_binround( mpnum *dst, mpnum *src, uint64 bit )
{
    uint32 bitval = two_to_the_31;
    uint32 unit = bit / 32;
    int isFirstBit = TRUE;
    int i;
    int round = 0;  /* -1 => down, 0 => even, 1 => up */
    uint32 valoffirstbit;
    uint32 firstunit = unit;
    bit = bit % 32;

    if( dst->allocated < src->len + 1 ) {
        MP_GROW( dst, src->len + 1 - dst->allocated );
    }

    MP_COPY( dst, src );

    /* handle special case where bit is 0 */
    if( unit == 0 && bit == 0 ) {
        return MP_NO_ERR;
    }

    for( i = 31; i > bit; i-- ) {
        bitval >>= 1;
    }
    valoffirstbit = bitval;

    /* determine rounding */
    for(;;) {
        /* decrement bit */
        if( bitval == 1 ) {
            if( unit == 0 ) break;
            unit--;
            bitval = two_to_the_31;
        } else {
            bitval >>= 1;
        }

        if( isFirstBit ) {
            if( !(src->num[unit] & bitval) ) {
                round = -1;  /* round down */
                break;
            }
            isFirstBit = FALSE;
        } else {
            if( src->num[unit] & bitval ) {
                round = 1; /* round up */
                break;
            }
        }
    }

    /* add 1 if necessary */
    if( (round == 1) || (round == 0 && (src->num[firstunit] & valoffirstbit) ) ) {
        mpnum temp;
        MP_INIT( &temp, 0 );
        if( firstunit >= temp.allocated ) {
            MP_GROW( &temp, firstunit - temp.allocated + 1);
        }
        mp_zero( &temp );
        temp.num[firstunit] = valoffirstbit;
        temp.len = firstunit + 1;
        mp_add( dst, dst, &temp );
        mp_free( &temp );
    }

    /* truncate */
    dst->num[firstunit] = dst->num[firstunit] & (ULONG_MAX - (valoffirstbit - 1));
    for( unit = firstunit;; unit--) {
        if( unit == 0 ) break;
        dst->num[unit-1] = 0;
    }

    return MP_NO_ERR;
}

/* Convert src to a real 4 byte number using mp_binround to round to the
 * nearest value. */
int mp_tofloat( uint8 *dst, mpnum *src )
{
    uint32 firstunit = src->len - 1;
    uint8  firstbit = 31;
    uint32 bit = two_to_the_31;
    uint64 exp;
    uint8 *temp;
    uint64 temp2;
    uint64 sig = 0;
    uint32 left = 0;
    uint8 i;
    mpnum rounded_num;

    memset( dst, 0, 4 );

    /* find first non-zero bit */
    if( src->len == 0 ) {
        return MP_NO_ERR;
    }
    for( ;; firstunit-- ) {
        if( src->num[firstunit] > 0 ) {
            break;
        }
        if( firstunit == 0 ) return MP_NO_ERR; /* src = 0 */
    }
    for( ;; firstbit-- ) {
        if( src->num[firstunit] & bit ) {
            break;
        }
        bit >>= 1;
    }

    /* round this number so that we have 24 significant bits */
    MP_INIT( &rounded_num, 0 );
    MP_COPY( &rounded_num, src );
    if( firstunit * 32 + firstbit > 23 ) {
        mp_binround( &rounded_num, &rounded_num, firstunit * 32 + firstbit - 23 );
        /* if it rounded up, the first bit may be up one */
        if( firstbit == 31 ) {
            if( rounded_num.num[firstunit+1] > 0 ) {
                firstunit++;
                firstbit = 0;
            }
        } else {
            if( rounded_num.num[firstunit] & (bit * 2) ) {
                firstbit++;
            }
        }
    }

    /* generate exponent */
    exp = (uint64)firstunit * 32 + firstbit + float_bias;
    if( exp > 255 ) exp = 255;
    temp = (uint8*)&exp;
    dst[3] = temp[0]>>1;
    dst[2] = temp[0]<<7;
    if( exp == 255 ) {
        /* too large => use +infinity */
        mp_free( &rounded_num );
        return MP_LOSS_OF_PRECISION;
    }

    /* generate significand */
    if( firstbit > 22 ) {
        left = 0;
        sig = (uint64)rounded_num.num[firstunit] >> (firstbit - 23);
    } else {
        left = 23 - firstbit;
        sig = (uint64)rounded_num.num[firstunit] << left;
    }
    if( firstunit > 0 && left > 0 ) {
        temp2 = (uint64)rounded_num.num[firstunit-1] >> (32-left);
        sig = sig | temp2;
        left = 0;
    }
    temp = (uint8*)&sig;
    dst[2] = dst[2] | (temp[2] & 127);
    for( i = 0; i < 2; i++ ) {
        dst[i] = temp[i];
    }

    mp_free( &rounded_num );
    return MP_NO_ERR;
}

/* Convert src to a real 8 byte number using mp_binround to round to the
 * nearest value. */
int mp_todouble( uint8 *dst, mpnum *src )
{
    uint32 firstunit = src->len - 1;
    uint8  firstbit = 31;
    uint32 bit = two_to_the_31;
    uint64 exp;
    uint8 *temp;
    uint64 temp2;
    uint64 sig = 0;
    uint32 left = 0;
    uint8 i;
    mpnum rounded_num;

    memset( dst, 0, 8 );
    if( src->len == 0 ) {
        return MP_NO_ERR;
    }
    for( ;; firstunit-- ) {
        if( src->num[firstunit] > 0 ) {
            break;
        }
        if( firstunit == 0 ) return MP_NO_ERR; /* src = 0 */
    }
    for( ;; firstbit-- ) {
        if( src->num[firstunit] & bit ) {
            break;
        }
        bit >>= 1;
    }

    /* round this number so that we have 53 significant bits */
    MP_INIT( &rounded_num, 0 );
    MP_COPY( &rounded_num, src );
    if( firstunit * 32 + firstbit > 52 ) {
        mp_binround( &rounded_num, &rounded_num, firstunit * 32 + firstbit - 52 );
        /* if it rounded up, the first bit may be up one */
        if( firstbit == 31 ) {
            if( rounded_num.num[firstunit+1] > 0 ) {
                firstunit++;
                firstbit = 0;
            }
        } else {
            if( rounded_num.num[firstunit] & (bit * 2) ) {
                firstbit++;
            }
        }
    }

    /* generate exponent */
    exp = (uint64)firstunit * 32 + firstbit + double_bias;
    if( exp > 2047 ) exp = 2047;
    temp = (uint8*)&exp;
    dst[7] = temp[1] << 4;
    dst[7] = dst[7] | ((temp[0] & 240) >> 4);
    dst[6] = (temp[0] & 15) << 4;
    if( exp == 2047 ) {
        /* too large => use +infinity */
        mp_free( &rounded_num );
        return MP_LOSS_OF_PRECISION;
    }

    /* generate significand */
    left = 52 - firstbit;
    sig = (uint64)rounded_num.num[firstunit] << left;
    if( firstunit > 0 ) {
        if( left > 32 ) {
            left -= 32;
            temp2 = (uint64)rounded_num.num[firstunit-1] << left;
            sig = sig | temp2;
        } else if( left < 32 ) {
            temp2 = (uint64)rounded_num.num[firstunit-1] >> (32-left);
            sig = sig | temp2;
            left = 0;
        } else {
            temp2 = (uint64)rounded_num.num[firstunit-1];
            sig = sig | temp2;
            left = 0;
        }
    }
    if( left > 0 && firstunit > 1 ) {
        temp2 = (uint64)rounded_num.num[firstunit-2] >> (32-left);
        sig = sig | temp2;
    }
    temp = (uint8*)&sig;
    for( i = 0; i < 6; i++ ) {
        dst[i] = temp[i];
    }
    dst[6] = dst[6] | (temp[6] & 15);

    mp_free( &rounded_num );
    return MP_NO_ERR;
}

/* Convert src to a real 10 byte number using mp_binround to round to the
 * nearest value. */
int mp_toextended( uint8 *dst, mpnum *src )
{
    uint32 firstunit = src->len - 1;
    uint8  firstbit = 31;
    uint32 bit = two_to_the_31;
    uint64 exp;
    uint8 *temp;
    uint64 temp2;
    uint64 sig = 0;
    uint32 left = 0;
    uint8 i;
    mpnum rounded_num;

    memset( dst, 0, 10 );

    /* find first non-zero bit */
    if( src->len == 0 ) {
        return MP_NO_ERR;
    }
    for( ;; firstunit-- ) {
        if( src->num[firstunit] > 0 ) {
            break;
        }
        if( firstunit == 0 ) return MP_NO_ERR; /* src = 0 */
    }
    for( ;; firstbit-- ) {
        if( src->num[firstunit] & bit ) {
            break;
        }
        bit >>= 1;
    }

    /* round this number so that we have 64 significant bits */
    MP_INIT( &rounded_num, 0 );
    MP_COPY( &rounded_num, src );
    if( firstunit * 32 + firstbit > 63 ) {
        mp_binround( &rounded_num, &rounded_num, firstunit * 32 + firstbit - 63 );
        /* if it rounded up, the first bit may be up one */
        if( firstbit == 31 ) {
            if( rounded_num.num[firstunit+1] > 0 ) {
                firstunit++;
                firstbit = 0;
            }
        } else {
            if( rounded_num.num[firstunit] & (bit * 2) ) {
                firstbit++;
            }
        }
    }

    /* generate exponent */
    exp = (uint64)firstunit * 32 + firstbit + extended_bias;
    if( exp > 32767 ) exp = 32767;
    temp = (uint8*)&exp;
    dst[9] = temp[1];
    dst[8] = temp[0];
    if( exp == 32767 ) {
        /* too large => use +infinity */
        mp_free( &rounded_num );
        return MP_LOSS_OF_PRECISION;
    }

    /* generate significand */
    left = 63 - firstbit;
    sig = (uint64)rounded_num.num[firstunit] << left;
    if( firstunit > 0 ) {
        if( left > 32 ) {
            left -= 32;
            temp2 = (uint64)rounded_num.num[firstunit-1] << left;
            sig = sig | temp2;
        } else if( left < 32 ) {
            temp2 = (uint64)rounded_num.num[firstunit-1] >> (32-left);
            sig = sig | temp2;
            left = 0;
        } else {
            temp2 = (uint64)rounded_num.num[firstunit-1];
            sig = sig | temp2;
            left = 0;
        }
    }
    if( left > 0 && firstunit > 1 ) {
        temp2 = (uint64)rounded_num.num[firstunit-2] >> (32-left);
        sig = sig | temp2;
    }
    temp = (uint8*)&sig;
    for( i = 0; i < 8; i++ ) {
        dst[i] = temp[i];
    }

    mp_free( &rounded_num );
    return MP_NO_ERR;
}

int mp_touint64( uint64 *dst, mpnum *src )
{
    mpnum limit;
    uint32 *temp = (uint32*)dst;

    MP_INIT( &limit, (uint64)0xffffffffffffffff );
    if( mp_gt( src, &limit ) ) {
        temp[0] = limit.num[0];
        temp[1] = limit.num[1];
        mp_free( &limit );
        return MP_TOO_LARGE;
    }
    temp[0] = src->num[0];
    if( src->len > 1 ) {
        temp[1] = src->num[1];
    } else {
        temp[1] = 0;
    }
    mp_free( &limit );
    return MP_NO_ERR;
}

/* returns size of num in bits */
int mp_bitsize( mpnum *num )
{
    uint32 i = num->len - 1;
    uint32 j = 32;
    uint32 bit = two_to_the_31;

    /* find first non-zero bit */
    if( num->len == 0 ) {
        return 0;
    }
    for( ;; i-- ) {
        if( num->num[i] > 0 ) {
            break;
        }
        if( i == 0 ) return 0; /* src = 0 */
    }
    for( ;; j-- ) {
        if( num->num[i] & bit ) {
            break;
        }
        bit >>= 1;
    }
    return i * 32 + j;
}