in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplAdd.java [143:532]
public static long bid64_add(final long x, final long y) {
long CA_w0, CA_w1, CT_w0, CT_w1, CT_new_w0, CT_new_w1;
long sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
long valid_x, valid_y;
long sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab, rem_a;
long saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
// long valid_x = unpack_BID64(&sign_x, &exponent_x, &coefficient_x, x);
{
sign_x = x & 0x8000000000000000L;
if ((x & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
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
// coefficient
long coeff = (x & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
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 {
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_x = coeff;
// get exponent
tmp = x >>> EXPONENT_SHIFT_LARGE64;
exponent_x = (int) (tmp & EXPONENT_MASK64);
valid_x = coeff;
}
}
}
// long valid_y = unpack_BID64(&sign_y, &exponent_y, &coefficient_y, y);
{
sign_y = y & 0x8000000000000000L;
if ((y & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
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
// coefficient
long coeff = (y & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
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 {
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_y = coeff;
// get exponent
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 Inf
if (((y & NAN_MASK64) == INFINITY_MASK64)) {
if (sign_x == (y & 0x8000000000000000L)) {
return coefficient_x;
}
// return NaN
{
return NAN_MASK64;
}
}
// check if y is NaN
if (((y & NAN_MASK64) == NAN_MASK64)) {
return coefficient_y & QUIET_MASK64;
}
// otherwise return +/-Inf
{
return coefficient_x;
}
}
// x is 0
{
if (((y & INFINITY_MASK64) != INFINITY_MASK64) && (coefficient_y != 0)) {
if (exponent_y <= exponent_x) {
return y;
}
}
}
}
if (valid_y == 0) {
// y is Inf. or NaN?
if (((y & INFINITY_MASK64) == INFINITY_MASK64)) {
return coefficient_y & QUIET_MASK64;
}
// y is 0
if (coefficient_x == 0) { // x==0
long res;
if (exponent_x <= exponent_y)
res = ((long) exponent_x) << 53;
else
res = ((long) exponent_y) << 53;
if (sign_x == sign_y)
res |= sign_x;
return res;
} else if (exponent_y >= exponent_x) {
return x;
}
}
// 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 ---
// version 2 (original)
long tempxi = Double.doubleToRawLongBits((double) coefficient_a);
bin_expon_ca = (int) ((tempxi & MASK_BINARY_EXPONENT) >>> 52) - 0x3ff;
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_w0[scale_ca])))
scale_ca++;
scale_k = 16 - scale_ca;
coefficient_a *= bid_power10_table_128_w0[scale_k];
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 ---
tempxi = Double.doubleToRawLongBits((double) coefficient_a);
bin_expon_ca = (int) ((tempxi & MASK_BINARY_EXPONENT) >>> 52) - 0x3ff;
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
// 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 fast_get_BID64_check_OF(sign_a, exponent_a, coefficient_a);
}
}
// 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_w0[diff_dec_expon];
// sign mask
sign_b = sign_b >> 63; // @AD: signed value shift
// apply sign to coeff. of b
coefficient_b = (coefficient_b + sign_b) ^ sign_b;
// apply sign to coefficient a
sign_a = sign_a >> 63; // @AD: signed value shift
coefficient_a = (coefficient_a + sign_a) ^ sign_a;
coefficient_a += coefficient_b;
// get sign
sign_s = coefficient_a >> 63; // @AD: signed value shift
coefficient_a = (coefficient_a + sign_s) ^ sign_s;
sign_s &= 0x8000000000000000L;
// coefficient_a < 10^16 ?
if ((UnsignedLong.isLess(coefficient_a, bid_power10_table_128_w0[MAX_FORMAT_DIGITS]))) {
return very_fast_get_BID64(sign_s, exponent_b, coefficient_a);
}
// otherwise rounding is necessary
// already know coefficient_a<10^19
// coefficient_a < 10^17 ?
if ((UnsignedLong.isLess(coefficient_a, bid_power10_table_128_w0[17])))
extra_digits = 1;
else if ((UnsignedLong.isLess(coefficient_a, bid_power10_table_128_w0[18])))
extra_digits = 2;
else
extra_digits = 3;
coefficient_a += bid_round_const_table_nearest[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;
// check whether we can take faster path
scale_ca = bid_estimate_decimal_digits[bin_expon_ca];
sign_ab = sign_a ^ sign_b;
sign_ab = sign_ab >> 63; // @AD: signed value shift
// T1 = 10^(16-diff_dec_expon)
T1 = bid_power10_table_128_w0[16 - diff_dec_expon];
// get number of digits in coefficient_a
if ((UnsignedLong.isGreaterOrEqual(coefficient_a, bid_power10_table_128_w0[scale_ca]))) {
scale_ca++;
}
scale_k = 16 - scale_ca;
// addition
saved_ca = coefficient_a - T1;
coefficient_a = saved_ca * bid_power10_table_128_w0[scale_k];
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_nearest[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
// this test 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;
//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 += rem_a * bid_power10_table_128_w0[diff_dec_expon];
} else
coefficient_a = (saved_ca - T1 - (T1 << 3)) * bid_power10_table_128_w0[scale_k - 1];
extra_digits++;
coefficient_b = saved_cb + 100000000000000000L + bid_round_const_table_nearest[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 = saved_ca * bid_power10_table_128_w0[scale_k + 1];
//extra_digits --;
exponent_b--;
coefficient_b = (saved_cb << 3) + (saved_cb << 1) + 100000000000000000L
+ bid_round_const_table_nearest[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 (BID_ROUNDING_TO_NEAREST == 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--;
}
}
return
fast_get_BID64_check_OF(sign_s, exponent_b + extra_digits, C64);
}