in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplFma.java [566:930]
static long /*BID_UINT64*/ bid_get_add64(final long /*BID_UINT64*/ sign_x, final int exponent_x, final long /*BID_UINT64*/ coefficient_x,
final long /*BID_UINT64*/ sign_y, final int exponent_y, final long /*BID_UINT64*/ coefficient_y
/*, final int rounding_mode, final FloatingPointStatusFlag fpsc*/) {
long /*BID_UINT128*/ CA_w0, CA_w1, CT_w0, CT_w1, CT_new_w0, CT_new_w1;
long /*BID_UINT64*/ sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab, rem_a;
long /*BID_UINT64*/ saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp, C64_new;
long /*int_double*/ tempx_i;
int exponent_a, exponent_b, diff_dec_expon;
int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
/*unsigned*/
int rmode/*, status*/;
// sort arguments by exponent
if (exponent_x <= exponent_y) {
sign_a = sign_y;
exponent_a = exponent_y;
coefficient_a = coefficient_y;
sign_b = sign_x;
exponent_b = exponent_x;
coefficient_b = coefficient_x;
} else {
sign_a = sign_x;
exponent_a = exponent_x;
coefficient_a = coefficient_x;
sign_b = sign_y;
exponent_b = exponent_y;
coefficient_b = coefficient_y;
}
// exponent difference
diff_dec_expon = exponent_a - exponent_b;
/* get binary coefficients of x and y */
//--- get number of bits in the coefficients of x and y ---
tempx_i = UnsignedLong.longToDoubleRawBits(coefficient_a);
bin_expon_ca = (int) (((tempx_i & MASK_BINARY_EXPONENT) >>> 52) - 0x3ff);
if (coefficient_a == 0) {
return get_BID64(sign_b, exponent_b, coefficient_b/*, rounding_mode, fpsc*/);
}
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
// normalize a to a 16-digit coefficient
scale_ca = bid_estimate_decimal_digits[bin_expon_ca];
if (UnsignedLong.isGreaterOrEqual(coefficient_a, bid_power10_table_128_BID_UINT128[(scale_ca << 1) /*+ 0*/]))
scale_ca++;
scale_k = 16 - scale_ca;
coefficient_a *= bid_power10_table_128_BID_UINT128[(scale_k << 1) /*+ 0*/];
diff_dec_expon -= scale_k;
exponent_a -= scale_k;
/* get binary coefficients of x and y */
//--- get number of bits in the coefficients of x and y ---
tempx_i = UnsignedLong.longToDoubleRawBits(coefficient_a);
bin_expon_ca = (int) (((tempx_i & MASK_BINARY_EXPONENT) >>> 52) - 0x3ff);
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
// if (coefficient_b != 0) {
// __set_status_flags(fpsc, BID_INEXACT_EXCEPTION);
// }
// if (((rounding_mode) & 3) != 0 && coefficient_b != 0) // not BID_ROUNDING_TO_NEAREST
// {
// switch (rounding_mode) {
// case BID_ROUNDING_DOWN:
// if (sign_b != 0) {
// coefficient_a -= ((((BID_SINT64) sign_a) >> 63) | 1);
// if (UnsignedLong.isLess(coefficient_a, 1000000000000000L)) {
// exponent_a--;
// coefficient_a = 9999999999999999L;
// } else if (UnsignedLong.isGreaterOrEqual(coefficient_a, 10000000000000000L)) {
// exponent_a++;
// coefficient_a = 1000000000000000L;
// }
// }
// break;
// case BID_ROUNDING_UP:
// if (sign_b == 0) {
// coefficient_a += ((((BID_SINT64) sign_a) >> 63) | 1);
// if (UnsignedLong.isLess(coefficient_a, 1000000000000000L) {
// exponent_a--;
// coefficient_a = 9999999999999999L;
// } else if (UnsignedLong.isGreaterOrEqual(coefficient_a, 10000000000000000L) {
// exponent_a++;
// coefficient_a = 1000000000000000L;
// }
// }
// break;
// default: // RZ
// if (sign_a != sign_b) {
// coefficient_a--;
// if (UnsignedLong.isLess(coefficient_a, 1000000000000000L) {
// exponent_a--;
// coefficient_a = 9999999999999999L;
// }
// }
// break;
// }
// } else
// check special case here
if ((coefficient_a == 1000000000000000L)
&& (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
&& (sign_a ^ sign_b) != 0
&& (UnsignedLong.isGreater(coefficient_b, 5000000000000000L))) {
coefficient_a = 9999999999999999L;
exponent_a--;
}
return get_BID64(sign_a, exponent_a, coefficient_a/*, rounding_mode, fpsc*/);
}
}
// test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
if (bin_expon_ca + bid_estimate_bin_expon[diff_dec_expon] < 60) {
// coefficient_a*10^(exponent_a-exponent_b)<2^63
// multiply by 10^(exponent_a-exponent_b)
coefficient_a *= bid_power10_table_128_BID_UINT128[(diff_dec_expon << 1) /*+ 0*/];
// sign mask
sign_b = (/*(BID_SINT64)*/sign_b) >> 63;
// apply sign to coeff. of b
coefficient_b = (coefficient_b + sign_b) ^ sign_b;
// apply sign to coefficient a
sign_a = (/*(BID_SINT64)*/sign_a) >> 63;
coefficient_a = (coefficient_a + sign_a) ^ sign_a;
coefficient_a += coefficient_b;
// get sign
sign_s = (/*(BID_SINT64)*/coefficient_a) >> 63;
coefficient_a = (coefficient_a + sign_s) ^ sign_s;
sign_s &= 0x8000000000000000L;
// coefficient_a < 10^16 ?
if (UnsignedLong.isLess(coefficient_a, bid_power10_table_128_BID_UINT128[(MAX_FORMAT_DIGITS << 1) /*+ 0*/])) {
// if (rounding_mode == BID_ROUNDING_DOWN && (coefficient_a == 0) && sign_a != sign_b)
// sign_s = 0x8000000000000000L;
return get_BID64(sign_s, exponent_b, coefficient_a/*, rounding_mode, fpsc*/);
}
// otherwise rounding is necessary
// already know coefficient_a<10^19
// coefficient_a < 10^17 ?
if (UnsignedLong.isLess(coefficient_a, bid_power10_table_128_BID_UINT128[(17 << 1) /*+ 0*/]))
extra_digits = 1;
else if (UnsignedLong.isLess(coefficient_a, bid_power10_table_128_BID_UINT128[(18 << 1) /*+ 0*/]))
extra_digits = 2;
else
extra_digits = 3;
rmode = /*rounding_mode*/BID_ROUNDING_TO_NEAREST;
// if (sign_s != 0 && (rmode == BID_ROUNDING_DOWN || rmode == BID_ROUNDING_UP) /*(uint)(rmode - 1) < 2*/)
// rmode = 3 - rmode;
coefficient_a += bid_round_const_table[rmode][extra_digits];
// get P*(2^M[extra_digits])/10^extra_digits
//__mul_64x64_to_128(CT, coefficient_a, bid_reciprocals10_64[extra_digits]);
{
final long __CY = bid_reciprocals10_64[extra_digits];
CT_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_a, __CY);
CT_w0 = coefficient_a * __CY;
}
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[extra_digits];
C64 = CT_w1 >>> amount;
} else {
// coefficient_a*10^(exponent_a-exponent_b) is large
sign_s = sign_a;
rmode = /*rounding_mode*/BID_ROUNDING_TO_NEAREST;
// if (sign_s != 0 && (rmode == BID_ROUNDING_DOWN || rmode == BID_ROUNDING_UP) /*(uint)(rmode - 1) < 2*/)
// rmode = 3 - rmode;
// check whether we can take faster path
scale_ca = bid_estimate_decimal_digits[bin_expon_ca];
sign_ab = sign_a ^ sign_b;
sign_ab = (/*(BID_SINT64)*/ sign_ab) >> 63;
// T1 = 10^(16-diff_dec_expon)
T1 = bid_power10_table_128_BID_UINT128[((16 - diff_dec_expon) << 1) /*+ 0*/];
// get number of digits in coefficient_a
//P_ca = bid_power10_table_128[scale_ca].w[0];
//P_ca_m1 = bid_power10_table_128[scale_ca-1].w[0];
if (UnsignedLong.isGreaterOrEqual(coefficient_a, bid_power10_table_128_BID_UINT128[(scale_ca << 1) /*+ 0*/])) {
scale_ca++;
//P_ca_m1 = P_ca;
//P_ca = bid_power10_table_128[scale_ca].w[0];
}
scale_k = 16 - scale_ca;
// apply sign
//Ts = (T1 + sign_ab) ^ sign_ab;
// test range of ca
//X = coefficient_a + Ts - P_ca_m1;
// addition
saved_ca = coefficient_a - T1;
coefficient_a = (long /*BID_SINT64*/) saved_ca *
(long /*BID_SINT64*/) bid_power10_table_128_BID_UINT128[(scale_k << 1) /*+ 0*/];
extra_digits = diff_dec_expon - scale_k;
// apply sign
saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
// add 10^16 and rounding constant
coefficient_b = saved_cb + 10000000000000000L + bid_round_const_table[rmode][extra_digits];
// get P*(2^M[extra_digits])/10^extra_digits
//__mul_64x64_to_128(CT, coefficient_b, bid_reciprocals10_64[extra_digits]);
{
final long __CY = bid_reciprocals10_64[extra_digits];
CT_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_b, __CY);
CT_w0 = coefficient_b * __CY;
}
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[extra_digits];
C0_64 = CT_w1 >>> amount;
// result coefficient
C64 = C0_64 + coefficient_a;
// filter out difficult (corner) cases
// the following test is equivalent to
// ( (initial_coefficient_a + Ts) < P_ca &&
// (initial_coefficient_a + Ts) > P_ca_m1 ),
// which ensures the number of digits in coefficient_a does not change
// after adding (the appropriately scaled and rounded) coefficient_b
if (UnsignedLong.isGreater(C64 - 1000000000000000L - 1, 9000000000000000L - 2)) {
if (UnsignedLong.isGreaterOrEqual(C64, 10000000000000000L)) {
// result has more than 16 digits
if (scale_k == 0) {
// must divide coeff_a by 10
saved_ca = saved_ca + T1;
//__mul_64x64_to_128(CA, saved_ca, 0x3333333333333334L);
CA_w1 = Mul64Impl.unsignedMultiplyHigh(saved_ca, 0x3333333333333334L);
// CA_w0 = saved_ca * 0x3333333333333334L; // @optimization
//reciprocals10_64[1]);
coefficient_a = CA_w1 >>> 1;
rem_a =
saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
coefficient_a = coefficient_a - T1;
saved_cb += /*90000000000000000 */ +rem_a *
bid_power10_table_128_BID_UINT128[(diff_dec_expon << 1) /*+ 0*/];
} else
coefficient_a = (long /*BID_SINT64*/) (saved_ca - T1 - (T1 << 3)) *
(long /*BID_SINT64*/) bid_power10_table_128_BID_UINT128[((scale_k - 1) << 1) /*+ 0*/];
extra_digits++;
coefficient_b =
saved_cb + 100000000000000000L +
bid_round_const_table[rmode][extra_digits];
// get P*(2^M[extra_digits])/10^extra_digits
//__mul_64x64_to_128(CT, coefficient_b, bid_reciprocals10_64[extra_digits]);
{
final long __CY = bid_reciprocals10_64[extra_digits];
CT_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_b, __CY);
CT_w0 = coefficient_b * __CY;
}
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[extra_digits];
C0_64 = CT_w1 >>> amount;
// result coefficient
C64 = C0_64 + coefficient_a;
} else if (UnsignedLong.isLessOrEqual(C64, 1000000000000000L)) {
// less than 16 digits in result
coefficient_a = (long/*BID_SINT64*/) saved_ca *
(long/*BID_SINT64*/) bid_power10_table_128_BID_UINT128[((scale_k + 1) << 1) /*+ 0*/];
//extra_digits --;
exponent_b--;
coefficient_b =
(saved_cb << 3) + (saved_cb << 1) + 100000000000000000L +
bid_round_const_table[rmode][extra_digits];
// get P*(2^M[extra_digits])/10^extra_digits
//__mul_64x64_to_128(CT_new, coefficient_b, bid_reciprocals10_64[extra_digits]);
{
final long __CY = bid_reciprocals10_64[extra_digits];
CT_new_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_b, __CY);
CT_new_w0 = coefficient_b * __CY;
}
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[extra_digits];
C0_64 = CT_new_w1 >>> amount;
// result coefficient
C64_new = C0_64 + coefficient_a;
if (UnsignedLong.isLess(C64_new, 10000000000000000L)) {
C64 = C64_new;
CT_w0 = CT_new_w0;
CT_w1 = CT_new_w1;
} else
exponent_b++;
}
}
}
// if (rmode == 0) //BID_ROUNDING_TO_NEAREST
if ((C64 & 1) != 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 = CT_w1 << (64 - amount);
// test whether fractional part is 0
if (remainder_h == 0 && UnsignedLong.isLess(CT_w0, bid_reciprocals10_64[extra_digits])) {
C64--;
}
}
// status = BID_INEXACT_EXCEPTION;
// get remainder
remainder_h = CT_w1 << (64 - amount);
// switch (rmode) {
// case BID_ROUNDING_TO_NEAREST:
// case BID_ROUNDING_TIES_AWAY:
// // test whether fractional part is 0
// if ((remainder_h == 0x8000000000000000L) && UnsignedLong.isLess(CT_w0, bid_reciprocals10_64[extra_digits]))
// status = BID_EXACT_STATUS;
// break;
// case BID_ROUNDING_DOWN:
// case BID_ROUNDING_TO_ZERO:
// if (remainder_h == 0 && UnsignedLong.isLess(CT_w0, bid_reciprocals10_64[extra_digits]))
// status = BID_EXACT_STATUS;
// break;
// default:
// // round up
// //__add_carry_out(tmp, carry, CT_w0, bid_reciprocals10_64[extra_digits]);
// {
// tmp = CT_w0 + bid_reciprocals10_64[extra_digits];
// carry = UnsignedLong.isLess(tmp, CT_w0) ? 1 : 0;
// }
//
// if (UnsignedLong.isGreaterOrEqual((remainder_h >>> (64 - amount)) + carry, (((long) 1) << amount)))
// status = BID_EXACT_STATUS;
// break;
// }
// __set_status_flags(fpsc, status);
return get_BID64(sign_s, exponent_b + extra_digits, C64/*, rounding_mode, fpsc*/);
}