in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplFma.java [43:537]
public static long /*BID_UINT64*/ bid64_fma(long /*BID_UINT64*/ x, long /*BID_UINT64*/ y, long /*BID_UINT64*/ z
/*, final int rnd_mode, final JavaImplParse.FloatingPointStatusFlag pfpsf*/) {
long /*BID_UINT128*/ P_w0, P_w1, CT_w0, CT_w1, CZ_w0, CZ_w1;
long /*BID_UINT64*/ sign_x, sign_y, coefficient_x, coefficient_y, sign_z, coefficient_z;
long /*BID_UINT64*/ C64, remainder_y, res;
long /*BID_UINT64*/ CYh, CY0L, T, valid_x, valid_y, valid_z;
//int_double tempx, tempy;
int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy, bin_expon_product;
int digits_p, bp, final_exponent, exponent_z, digits_z, ez, ey, scale_z/*, uf_status*/;
//valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
//valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
// 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
final 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
// 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
final long 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
final 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
// 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
final long tmp = y >>> EXPONENT_SHIFT_LARGE64;
exponent_y = (int) (tmp & EXPONENT_MASK64);
valid_y = coeff;
}
}
}
//valid_z = unpack_BID64 (&sign_z, &exponent_z, &coefficient_z, z);
{
sign_z = z & 0x8000000000000000L;
if ((z & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
final long tmp = z >>> EXPONENT_SHIFT_SMALL64;
exponent_z = (int) (tmp & EXPONENT_MASK64);
// coefficient
coefficient_z = (z & SMALL_COEFF_MASK64);
valid_z = coefficient_z;
} else {
// special encodings
// coefficient
long coeff = (z & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
if ((z & INFINITY_MASK64) == INFINITY_MASK64) {
exponent_z = 0;
coefficient_z = z & 0xfe03ffffffffffffL;
if ((z & 0x0003ffffffffffffL) >= 1000000000000000L)
coefficient_z = z & 0xfe00000000000000L;
if ((z & NAN_MASK64) == INFINITY_MASK64)
coefficient_z = z & SINFINITY_MASK64;
valid_z = 0; // NaN or Infinity
} else {
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_z = coeff;
// get exponent
final long tmp = z >>> EXPONENT_SHIFT_LARGE64;
exponent_z = (int) (tmp & EXPONENT_MASK64);
valid_z = coeff;
}
}
}
// unpack arguments, check for NaN, Infinity, or 0
if (valid_x == 0 || valid_y == 0 || valid_z == 0) {
if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
// if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
// check first for non-canonical NaN payload
y = y & 0xfe03ffffffffffffL; // clear G6-G12
if ((y & 0x0003ffffffffffffL) >= 1000000000000000L) {
y = y & 0xfe00000000000000L; // clear G6-G12 and the payload bits
}
if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
// set invalid flag
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// return quiet (y)
res = y & 0xfdffffffffffffffL;
} else { // y is QNaN
// return y
res = y;
// if z = SNaN or x = SNaN signal invalid exception
// if ((z & MASK_SNAN) == MASK_SNAN || (x & MASK_SNAN) == MASK_SNAN) {
// // set invalid flag
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// }
}
return res;
} else if ((z & MASK_NAN) == MASK_NAN) { // z is NAN
// if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
// check first for non-canonical NaN payload
z = z & 0xfe03ffffffffffffL; // clear G6-G12
if ((z & 0x0003ffffffffffffL) >= 1000000000000000L) {
z = z & 0xfe00000000000000L; // clear G6-G12 and the payload bits
}
if ((z & MASK_SNAN) == MASK_SNAN) { // z is SNAN
// set invalid flag
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// return quiet (z)
res = z & 0xfdffffffffffffffL;
} else { // z is QNaN
// return z
res = z;
// if x = SNaN signal invalid exception
// if ((x & MASK_SNAN) == MASK_SNAN) {
// // set invalid flag
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// }
}
return res;
} else if ((x & MASK_NAN) == MASK_NAN) { // x is NAN
// if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
// check first for non-canonical NaN payload
x = x & 0xfe03ffffffffffffL; // clear G6-G12
if ((x & 0x0003ffffffffffffL) >= 1000000000000000L) {
x = x & 0xfe00000000000000L; // clear G6-G12 and the payload bits
}
if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN
// set invalid flag
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// return quiet (x)
res = x & 0xfdffffffffffffffL;
} else { // x is QNaN
// return x
res = x; // clear out G[6]-G[16]
}
return res;
}
if (valid_x == 0) {
// x is Inf. or 0
// x is Infinity?
if ((x & 0x7800000000000000L) == 0x7800000000000000L) {
// check if y is 0
if (coefficient_y == 0) {
// y==0, return NaN
// if ((z & 0x7e00000000000000L) != 0x7c00000000000000L)
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
return 0x7c00000000000000L;
}
// test if z is Inf of oposite sign
if (((z & 0x7c00000000000000L) == 0x7800000000000000L)
&& (((x ^ y) ^ z) & 0x8000000000000000L) != 0) {
// return NaN
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
return 0x7c00000000000000L;
}
// otherwise return +/-Inf
return ((x ^ y) & 0x8000000000000000L) | 0x7800000000000000L;
}
// x is 0
if (((y & 0x7800000000000000L) != 0x7800000000000000L)
&& ((z & 0x7800000000000000L) != 0x7800000000000000L)) {
if (coefficient_z != 0) {
exponent_y = exponent_x - DECIMAL_EXPONENT_BIAS + exponent_y;
sign_z = z & 0x8000000000000000L;
if (exponent_y >= exponent_z)
return z;
return add_zero64(exponent_y, sign_z, exponent_z, coefficient_z/*, rnd_mode, pfpsf*/);
}
}
}
if (valid_y == 0) {
// y is Inf. or 0
// y is Infinity?
if ((y & 0x7800000000000000L) == 0x7800000000000000L) {
// check if x is 0
if (coefficient_x == 0) {
// y==0, return NaN
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
return 0x7c00000000000000L;
}
// test if z is Inf of oposite sign
if (((z & 0x7c00000000000000L) == 0x7800000000000000L)
&& (((x ^ y) ^ z) & 0x8000000000000000L) != 0) {
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
// return NaN
return 0x7c00000000000000L;
}
// otherwise return +/-Inf
return ((x ^ y) & 0x8000000000000000L) | 0x7800000000000000L;
}
// y is 0
if (((z & 0x7800000000000000L) != 0x7800000000000000L)) {
if (coefficient_z != 0) {
exponent_y += exponent_x - DECIMAL_EXPONENT_BIAS;
sign_z = z & 0x8000000000000000L;
if (exponent_y >= exponent_z)
return z;
return
add_zero64(exponent_y, sign_z, exponent_z, coefficient_z/*, rnd_mode, pfpsf*/);
}
}
}
if (valid_z == 0) {
// y is Inf. or 0
// test if y is NaN/Inf
if ((z & 0x7800000000000000L) == 0x7800000000000000L) {
return (coefficient_z & QUIET_MASK64);
}
// z is 0, return x*y
if ((coefficient_x == 0) || (coefficient_y == 0)) {
//0+/-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;
if (exponent_x <= exponent_z)
res = ((long) exponent_x) << 53;
else
res = ((long) exponent_z) << 53;
if ((sign_x ^ sign_y) == sign_z)
res |= sign_z;
// else if (rnd_mode == BID_ROUNDING_DOWN)
// res |= 0x8000000000000000L;
return res;
}
}
}
/* get binary coefficients of x and y */
//--- get number of bits in the coefficients of x and y ---
// version 2 (original)
long tempx_i = Double.doubleToRawLongBits((double) coefficient_x);
bin_expon_cx = (int) ((tempx_i & MASK_BINARY_EXPONENT) >>> 52);
final long tempy_i = Double.doubleToRawLongBits((double) coefficient_y);
bin_expon_cy = (int) ((tempy_i & 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;
final_exponent = exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS;
if ((final_exponent > 0) || (coefficient_z == 0)) {
return bid_get_add64(sign_x ^ sign_y,
final_exponent, C64, sign_z, exponent_z, coefficient_z/*, rnd_mode, pfpsf*/);
} else {
P_w0 = C64;
P_w1 = 0;
extra_digits = 0;
}
} else {
if (coefficient_z == 0) {
return bid64_mul(x, y/*, rnd_mode*/);
}
// 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];
final int bid_power10_table_128_index = digits_p << 1;
if (!(__unsigned_compare_gt_128(bid_power10_table_128_BID_UINT128[bid_power10_table_128_index],
bid_power10_table_128_BID_UINT128[bid_power10_table_128_index + 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;
}
if (UnsignedInteger.isGreaterOrEqual(final_exponent, 3 * 256)) {
if (final_exponent < 0) {
//--- get number of bits in the coefficients of z ---
tempx_i = Double.doubleToRawLongBits((double) coefficient_z);
bin_expon_cx = (int) ((tempx_i & MASK_BINARY_EXPONENT) >>> 52) - 0x3ff;
// get number of decimal digits in the coeff_x
digits_z = bid_estimate_decimal_digits[bin_expon_cx];
if (UnsignedLong.isGreaterOrEqual(coefficient_z, bid_power10_table_128_BID_UINT128[(digits_z << 1) /*+ 0*/]))
digits_z++;
// underflow
if ((final_exponent + 16 < 0)
|| (exponent_z + digits_z > 33 + final_exponent)) {
return BID_normalize(sign_z, exponent_z, coefficient_z, sign_x ^ sign_y, 1/*, rnd_mode, pfpsf*/);
}
ez = exponent_z + digits_z - 16;
if (ez < 0)
ez = 0;
scale_z = exponent_z - ez;
coefficient_z *= bid_power10_table_128_BID_UINT128[(scale_z << 1) /*+ 0*/];
ey = final_exponent - extra_digits;
extra_digits = ez - ey;
if (extra_digits > 17) {
CYh = __truncate(P_w0, P_w1, 16);
// get remainder
T = bid_power10_table_128_BID_UINT128[(16 << 1) /*+ 0*/];
//__mul_64x64_to_64(CY0L, CYh, T);
CY0L = CYh * T;
remainder_y = P_w0 - CY0L;
extra_digits -= 16;
P_w0 = CYh;
P_w1 = 0;
} else
remainder_y = 0;
// align coeff_x, CYh
//__mul_64x64_to_128(CZ, coefficient_z, bid_power10_table_128_flat[(extra_digits << 1) /*+ 0*/]);
{
final long __CY = bid_power10_table_128_BID_UINT128[(extra_digits << 1) /*+ 0*/];
CZ_w1 = Mul64Impl.unsignedMultiplyHigh(coefficient_z, __CY);
CZ_w0 = coefficient_z * __CY;
}
if (sign_z == (sign_y ^ sign_x)) {
//__add_128_128(CT, CZ, P);
{
long Q128_w0, Q128_w1;
Q128_w1 = CZ_w1 + P_w1;
Q128_w0 = P_w0 + CZ_w0;
if ((UnsignedLong.isLess(Q128_w0, P_w0)))
Q128_w1++;
CT_w1 = Q128_w1;
CT_w0 = Q128_w0;
}
final int bid_power10_table_128_index = (16 + extra_digits) << 1;
if (__unsigned_compare_ge_128(CT_w0, CT_w1,
bid_power10_table_128_BID_UINT128[bid_power10_table_128_index],
bid_power10_table_128_BID_UINT128[bid_power10_table_128_index + 1])) {
extra_digits++;
ez++;
}
} else {
if (remainder_y != 0 && (__unsigned_compare_ge_128(CZ_w0, CZ_w1, P_w0, P_w1))) {
P_w0++;
if (P_w0 == 0)
P_w1++;
}
//__sub_128_128(CT, CZ, P);
{
CT_w1 = CZ_w1 - P_w1;
CT_w0 = CZ_w0 - P_w0;
if (UnsignedLong.isLess(CZ_w0, P_w0))
CT_w1--;
}
if ((/*(BID_SINT64)*/ CT_w1) < 0) {
sign_z = sign_y ^ sign_x;
CT_w0 = -CT_w0;
CT_w1 = -CT_w1;
if (CT_w0 != 0)
CT_w1--;
} else if ((CT_w1 | CT_w0) == 0)
sign_z = /*(rnd_mode != BID_ROUNDING_DOWN) ?*/ 0 /*: 0x8000000000000000L*/;
final int bid_power10_table_128_index = (15 + extra_digits) << 1;
if (ez != 0 && __unsigned_compare_gt_128(bid_power10_table_128_BID_UINT128[bid_power10_table_128_index],
bid_power10_table_128_BID_UINT128[bid_power10_table_128_index + 1], CT_w0, CT_w1)) {
extra_digits--;
ez--;
}
}
// uf_status = 0;
// final int bid_power10_table_128_index = (extra_digits + 15) << 1;
// if (ez == 0 && __unsigned_compare_gt_128(bid_power10_table_128_BID_UINT128[bid_power10_table_128_index],
// bid_power10_table_128_BID_UINT128[bid_power10_table_128_index + 1], CT_w0, CT_w1)) {
// uf_status = BID_UNDERFLOW_EXCEPTION;
// }
return __bid_full_round64_remainder(sign_z, ez - extra_digits, CT_w0, CT_w1,
extra_digits, remainder_y/*, rnd_mode, pfpsf, uf_status*/);
} else {
if ((sign_z == (sign_x ^ sign_y)) || (final_exponent > 3 * 256 + 15)) {
return fast_get_BID64_check_OF(sign_x ^ sign_y, final_exponent,
1000000000000000L/*, rnd_mode, pfpsf*/);
}
}
}
if (extra_digits > 0) {
return bid_get_add128(sign_z, exponent_z, coefficient_z, sign_x ^ sign_y,
final_exponent, P_w0, P_w1, extra_digits/*, rnd_mode, pfpsf*/);
}
// go to convert_format and exit
else {
C64 = P_w0;
return bid_get_add64(sign_x ^ sign_y, exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
sign_z, exponent_z, coefficient_z/*, rnd_mode, pfpsf*/);
}
}