in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplMul.java [43:506]
public static long bid64_mul(final long x, final long y) {
long P_w0, P_w1, C128_w0, C128_w1, Q_high_w0, Q_high_w1, Q_low_w0, Q_low_w1;
long C64, remainder_h;
int extra_digits, bin_expon_cx, bin_expon_cy, bin_expon_product;
int digits_p, bp, amount, amount2, final_exponent, round_up;
//valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
final long sign_x;
int exponent_x;
long coefficient_x;
final long valid_x;
{
sign_x = x & 0x8000000000000000L;
if ((x & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
long tmp = x >>> EXPONENT_SHIFT_SMALL64;
exponent_x = (int) (tmp & EXPONENT_MASK64);
// coefficient
coefficient_x = (x & SMALL_COEFF_MASK64);
valid_x = coefficient_x;
} else {
// special encodings
if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
exponent_x = 0;
coefficient_x = x & 0xfe03ffffffffffffL;
if ((x & 0x0003ffffffffffffL) >= 1000000000000000L)
coefficient_x = x & 0xfe00000000000000L;
if ((x & NAN_MASK64) == INFINITY_MASK64)
coefficient_x = x & SINFINITY_MASK64;
valid_x = 0; // NaN or Infinity
} else {
// coefficient
long coeff = (x & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_x = coeff;
// get exponent
long tmp = x >>> EXPONENT_SHIFT_LARGE64;
exponent_x = (int) (tmp & EXPONENT_MASK64);
valid_x = coeff;
}
}
}
//valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
long sign_y;
int exponent_y;
long coefficient_y;
final long valid_y;
{
sign_y = y & 0x8000000000000000L;
if ((y & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
long tmp = y >>> EXPONENT_SHIFT_SMALL64;
exponent_y = (int) (tmp & EXPONENT_MASK64);
// coefficient
coefficient_y = (y & SMALL_COEFF_MASK64);
valid_y = coefficient_y;
} else {
// special encodings
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
exponent_y = 0;
coefficient_y = y & 0xfe03ffffffffffffL;
if ((y & 0x0003ffffffffffffL) >= 1000000000000000L)
coefficient_y = y & 0xfe00000000000000L;
if ((y & NAN_MASK64) == INFINITY_MASK64)
coefficient_y = y & SINFINITY_MASK64;
valid_y = 0; // NaN or Infinity
} else {
// coefficient
long coeff = (y & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_y = coeff;
// get exponent
long tmp = y >>> EXPONENT_SHIFT_LARGE64;
exponent_y = (int) (tmp & EXPONENT_MASK64);
valid_y = coeff;
}
}
}
// unpack arguments, check for NaN or Infinity
if (valid_x == 0) {
// x is Inf. or NaN
// test if x is NaN
if ((x & NAN_MASK64) == NAN_MASK64) {
return coefficient_x & QUIET_MASK64;
}
// x is Infinity?
if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
// check if y is 0
if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y == 0) {
// y==0 , return NaN
return NAN_MASK64;
}
// check if y is NaN
if ((y & NAN_MASK64) == NAN_MASK64)
// y==NaN , return NaN
return coefficient_y & QUIET_MASK64;
// otherwise return +/-Inf
return (((x ^ y) & 0x8000000000000000L) | INFINITY_MASK64);
}
// x is 0
if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
exponent_y = ((int) (y >>> 51)) & 0x3ff;
else
exponent_y = ((int) (y >>> 53)) & 0x3ff;
sign_y = y & 0x8000000000000000L;
exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
if (exponent_x > DECIMAL_MAX_EXPON_64)
exponent_x = DECIMAL_MAX_EXPON_64;
else if (exponent_x < 0)
exponent_x = 0;
return (sign_x ^ sign_y) | (((long) exponent_x) << 53);
}
}
if (valid_y == 0) {
// y is Inf. or NaN
// test if y is NaN
if ((y & NAN_MASK64) == NAN_MASK64) {
return coefficient_y & QUIET_MASK64;
}
// y is Infinity?
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
// check if x is 0
if (coefficient_x == 0) {
// x==0, return NaN
return NAN_MASK64;
}
// otherwise return +/-Inf
return (((x ^ y) & 0x8000000000000000L) | INFINITY_MASK64);
}
// y is 0
exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
if (exponent_x > DECIMAL_MAX_EXPON_64)
exponent_x = DECIMAL_MAX_EXPON_64;
else if (exponent_x < 0)
exponent_x = 0;
return ((sign_x ^ sign_y) | (((long) exponent_x) << 53));
}
//--- get number of bits in the coefficients of x and y ---
// version 2 (original)
long tempxi = Double.doubleToRawLongBits((double) coefficient_x);
bin_expon_cx = (int) ((tempxi & MASK_BINARY_EXPONENT) >>> 52);
long tempyi = Double.doubleToRawLongBits((double) coefficient_y);
bin_expon_cy = (int) ((tempyi & MASK_BINARY_EXPONENT) >>> 52);
// magnitude estimate for coefficient_x*coefficient_y is
// 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
bin_expon_product = bin_expon_cx + bin_expon_cy;
// check if coefficient_x*coefficient_y<2^(10*k+3)
// equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {
// easy multiply
C64 = coefficient_x * coefficient_y;
return get_BID64_small_mantissa(sign_x ^ sign_y,
exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64);
} else {
// get 128-bit product: coefficient_x*coefficient_y
//__mul_64x64_to_128(P, coefficient_x, coefficient_y);
P_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_x, coefficient_y);
P_w0 = coefficient_x * coefficient_y;
// tighten binary range of P: leading bit is 2^bp
// unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;
//__tight_bin_range_128(bp, P, bin_expon_product); // tighten exponent range
{
final long __P_w0 = P_w0;
final long __P_w1 = P_w1;
final int __bin_expon = bin_expon_product;
long M = 1;
bp = __bin_expon;
if (bp < 63) {
M <<= bp + 1;
if ((UnsignedLong.isGreaterOrEqual(__P_w0, M)))
bp++;
} else if (bp > 64) {
M <<= bp + 1 - 64;
if (((UnsignedLong.isGreater(__P_w1, M))) || (__P_w1 == M && __P_w0 != 0))
bp++;
} else if (__P_w1 != 0)
bp++;
}
// get number of decimal digits in the product
digits_p = bid_estimate_decimal_digits[bp];
if (!(__unsigned_compare_gt_128(
bid_power10_table_128_BID_UINT128[(digits_p << 1) + 0],
bid_power10_table_128_BID_UINT128[(digits_p << 1) + 1], P_w0, P_w1)))
digits_p++; // if bid_power10_table_128[digits_p] <= P
// determine number of decimal digits to be rounded out
extra_digits = digits_p - MAX_FORMAT_DIGITS;
final_exponent =
exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;
round_up = 0;
if ((/*UnsignedInteger.compare*/(final_exponent) + Integer.MIN_VALUE >= (3 * 256) + Integer.MIN_VALUE)) {
if (final_exponent < 0) {
// underflow
if (final_exponent + 16 < 0) {
return sign_x ^ sign_y;
}
extra_digits -= final_exponent;
final_exponent = 0;
if (extra_digits > 17) {
//__mul_128x128_full(Q_high, Q_low, P, bid_reciprocals10_128[16]);
{
final long _A_w0 = P_w0;
final long _A_w1 = P_w1;
final long _B_w0 = bid_reciprocals10_128_BID_UINT128[(16 << 1) /*+ 0*/];
final long _B_w1 = bid_reciprocals10_128_BID_UINT128[(16 << 1) + 1];
long _ALBL_w0, _ALBH_w0, _AHBL_w0, _AHBH_w0, _QM_w0, _QM2_w0;
long _ALBL_w1, _ALBH_w1, _AHBL_w1, _AHBH_w1, _QM_w1, _QM2_w1;
//__mul_64x64_to_128(ALBH, (A)_w0, (B)_w1);
_ALBH_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w0, _B_w1);
_ALBH_w0 = _A_w0 * _B_w1;
//__mul_64x64_to_128(AHBL, (B)_w0, (A)_w1);
_AHBL_w1 = Mul64Impl.unsignedMultiplyHigh(_B_w0, _A_w1);
_AHBL_w0 = _B_w0 * _A_w1;
//__mul_64x64_to_128(ALBL, (A)_w0, (B)_w0);
_ALBL_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w0, _B_w0);
_ALBL_w0 = _A_w0 * _B_w0;
//__mul_64x64_to_128(AHBH, (A)_w1,(B)_w1);
_AHBH_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w1, _B_w1);
_AHBH_w0 = _A_w1 * _B_w1;
//__add_128_128(QM, ALBH, AHBL); // add 128-bit value to 128-bit assume no carry-out
{
final long __A128_w0 = _ALBH_w0;
final long __A128_w1 = _ALBH_w1;
final long __B128_w0 = _AHBL_w0;
final long __B128_w1 = _AHBL_w1;
long Q128_w0, Q128_w1;
Q128_w1 = __A128_w1 + __B128_w1;
Q128_w0 = __B128_w0 + __A128_w0;
if ((UnsignedLong.isLess(Q128_w0, __B128_w0)))
Q128_w1++;
_QM_w1 = Q128_w1;
_QM_w0 = Q128_w0;
}
Q_low_w0 = _ALBL_w0;
//__add_128_64(_QM2, QM, ALBL_w1); // add 64-bit value to 128-bit
{
final long __A128_w0 = _QM_w0;
final long __A128_w1 = _QM_w1;
final long __B64 = _ALBL_w1;
long __R64H;
__R64H = __A128_w1;
_QM2_w0 = __B64 + __A128_w0;
if ((UnsignedLong.isLess(_QM2_w0, __B64)))
__R64H++;
_QM2_w1 = __R64H;
}
//__add_128_64((_Qh), AHBH, QM2_w1); // add 64-bit value to 128-bit
{
final long __A128_w0 = _AHBH_w0;
final long __A128_w1 = _AHBH_w1;
final long __B64 = _QM2_w1;
long __R64H;
__R64H = __A128_w1;
Q_high_w0 = __B64 + __A128_w0;
if ((UnsignedLong.isLess(Q_high_w0, __B64)))
__R64H++;
Q_high_w1 = __R64H;
}
Q_low_w1 = _QM2_w0;
}
amount = bid_recip_scale[16];
//__shr_128(P, Q_high, amount);
{
final long __A_w0 = Q_high_w0;
final long __A_w1 = Q_high_w1;
final long __k = amount;
P_w0 = __A_w0 >>> __k;
P_w0 |= __A_w1 << (64 - __k);
P_w1 = __A_w1 >>> __k;
}
// get sticky bits
amount2 = 64 - amount;
remainder_h = 0;
remainder_h--;
remainder_h >>>= amount2;
remainder_h = remainder_h & Q_high_w0;
extra_digits -= 16;
if (remainder_h != 0 || ((UnsignedLong.isGreater(Q_low_w1, bid_reciprocals10_128_BID_UINT128[(16 << 1) + 1]))
|| (Q_low_w1 == bid_reciprocals10_128_BID_UINT128[(16 << 1) + 1]
&& (UnsignedLong.isGreaterOrEqual(Q_low_w0, bid_reciprocals10_128_BID_UINT128[(16 << 1) + 0]))))) {
round_up = 1;
// __set_status_flags(pfpsf, BID_UNDERFLOW_EXCEPTION | BID_INEXACT_EXCEPTION);
P_w0 = (P_w0 << 3) + (P_w0 << 1);
P_w0 |= 1;
extra_digits++;
}
}
} else {
return fast_get_BID64_check_OF(sign_x ^ sign_y, final_exponent,
1000000000000000L/*, BID_ROUNDING_TO_NEAREST, pfpsf*/);
}
}
if (extra_digits > 0) {
// will divide by 10^(digits_p - 16)
// add a constant to P, depending on rounding mode
// 0.5*10^(digits_p - 16) for round-to-nearest
//__add_128_64(P, P, bid_round_const_table[rmode][extra_digits]); // add 64-bit value to 128-bit
{
final long __A128_w0 = P_w0;
final long __A128_w1 = P_w1;
final long __B64 = bid_round_const_table_nearest[extra_digits];
long __R64H;
__R64H = __A128_w1;
P_w0 = __B64 + __A128_w0;
if ((UnsignedLong.isLess(P_w0, __B64)))
__R64H++;
P_w1 = __R64H;
}
// get P*(2^M[extra_digits])/10^extra_digits
// __mul_128x128_full(Q_high, Q_low, P, bid_reciprocals10_128[extra_digits]);
{
final long _A_w0 = P_w0;
final long _A_w1 = P_w1;
final long _B_w0 = bid_reciprocals10_128_BID_UINT128[(extra_digits << 1) + 0];
final long _B_w1 = bid_reciprocals10_128_BID_UINT128[(extra_digits << 1) + 1];
long _ALBL_w0, _ALBH_w0, _AHBL_w0, _AHBH_w0, _QM_w0, _QM2_w0;
long _ALBL_w1, _ALBH_w1, _AHBL_w1, _AHBH_w1, _QM_w1, _QM2_w1;
//__mul_64x64_to_128(ALBH, (A)_w0, (B)_w1);
_ALBH_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w0, _B_w1);
_ALBH_w0 = _A_w0 * _B_w1;
//__mul_64x64_to_128(AHBL, (B)_w0, (A)_w1);
_AHBL_w1 = Mul64Impl.unsignedMultiplyHigh(_B_w0, _A_w1);
_AHBL_w0 = _B_w0 * _A_w1;
//__mul_64x64_to_128(ALBL, (A)_w0, (B)_w0);
_ALBL_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w0, _B_w0);
_ALBL_w0 = _A_w0 * _B_w0;
//__mul_64x64_to_128(AHBH, (A)_w1,(B)_w1);
_AHBH_w1 = Mul64Impl.unsignedMultiplyHigh(_A_w1, _B_w1);
_AHBH_w0 = _A_w1 * _B_w1;
//__add_128_128(QM, ALBH, AHBL); // add 128-bit value to 128-bit assume no carry-out
{
final long __A128_w0 = _ALBH_w0;
final long __A128_w1 = _ALBH_w1;
final long __B128_w0 = _AHBL_w0;
final long __B128_w1 = _AHBL_w1;
long Q128_w0, Q128_w1;
Q128_w1 = __A128_w1 + __B128_w1;
Q128_w0 = __B128_w0 + __A128_w0;
if ((UnsignedLong.isLess(Q128_w0, __B128_w0)))
Q128_w1++;
_QM_w1 = Q128_w1;
_QM_w0 = Q128_w0;
}
Q_low_w0 = _ALBL_w0;
//__add_128_64(_QM2, QM, ALBL_w1); // add 64-bit value to 128-bit
{
final long __A128_w0 = _QM_w0;
final long __A128_w1 = _QM_w1;
final long __B64 = _ALBL_w1;
long __R64H;
__R64H = __A128_w1;
_QM2_w0 = __B64 + __A128_w0;
if ((UnsignedLong.isLess(_QM2_w0, __B64)))
__R64H++;
_QM2_w1 = __R64H;
}
//__add_128_64((_Qh), AHBH, QM2_w1); // add 64-bit value to 128-bit
{
final long __A128_w0 = _AHBH_w0;
final long __A128_w1 = _AHBH_w1;
final long __B64 = _QM2_w1;
long __R64H;
__R64H = __A128_w1;
Q_high_w0 = __B64 + __A128_w0;
if ((UnsignedLong.isLess(Q_high_w0, __B64)))
__R64H++;
Q_high_w1 = __R64H;
}
Q_low_w1 = _QM2_w0;
}
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
amount = bid_recip_scale[extra_digits];
//__shr_128(C128, Q_high, amount);
{
final long __A_w0 = Q_high_w0;
final long __A_w1 = Q_high_w1;
final long __k = amount;
C128_w0 = __A_w0 >>> __k;
C128_w0 |= __A_w1 << (64 - __k);
C128_w1 = __A_w1 >>> __k;
}
C64 = C128_w0;
/*if (BID_ROUNDING_TO_NEAREST == 0)*/ //BID_ROUNDING_TO_NEAREST
if ((C64 & 1) != 0 && round_up == 0) {
// check whether fractional part of initial_P/10^extra_digits
// is exactly .5
// this is the same as fractional part of
// (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
// get remainder
remainder_h = Q_high_w0 << (64 - amount);
// test whether fractional part is 0
if (remainder_h == 0
&& ((UnsignedLong.isLess(Q_low_w1, bid_reciprocals10_128_BID_UINT128[(extra_digits << 1) + 1]))
|| (Q_low_w1 == bid_reciprocals10_128_BID_UINT128[(extra_digits << 1) + 1]
&& (UnsignedLong.isLess(Q_low_w0, bid_reciprocals10_128_BID_UINT128[(extra_digits << 1) + 0]))))) {
C64--;
}
}
// convert to BID and return
return fast_get_BID64_check_OF(sign_x ^ sign_y, final_exponent, C64/*, BID_ROUNDING_TO_NEAREST*/);
}
// go to convert_format and exit
C64 = P_w0;
return get_BID64(sign_x ^ sign_y,
exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64);
}
}