in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplCastBinary64.java [311:508]
public static double bid64_to_binary64(final long /*BID_UINT64*/ x, final int rnd_mode/*, final JavaImplParse.FloatingPointStatusFlag pfpsf*/) {
long /*BID_UINT64*/ c_prov;
long /*BID_UINT128*/ c_w0, c_w1;
long /*BID_UINT128*/ m_min_w0, m_min_w1;
int s, e, k, e_out;
long /*BID_UINT256*/ r_w0, r_w1, r_w2, r_w3;
long /*BID_UINT384*/ z_w0, z_w1, z_w2, z_w3, z_w4, z_w5;
//unpack_bid64 (x, s, e, k, (c_w1), return_binary64_zero (s), return_binary64_inf (s), return_binary64_nan);
//#define unpack_bid64(x,s,e,k,c,zero,inf,nan)
{
s = (int) (x >>> 63);
if ((x & (3L << 61)) == (3L << 61)) {
if ((x & (0xFL << 59)) == (0xFL << 59)) {
if ((x & (0x1FL << 58)) != (0x1FL << 58))
return return_binary64_inf(s);
// if ((x & (1L << 57)) != 0)
// __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
return return_binary64_nan(s, (x & 0x0003ffffffffffffL) >= 1000000000000000L ? 0 : (x << 14), 0L);
}
e = (int) (((x >>> 51) & ((1L << 10) - 1)) - 398);
c_w1 = (1L << 53) + (x & ((1L << 51) - 1));
if (UnsignedLong.isGreater(c_w1, 9999999999999999L))
return return_binary64_zero(s);
k = 0;
} else {
e = (int) (((x >>> 53) & ((1L << 10) - 1)) - 398);
c_w1 = x & ((1L << 53) - 1);
if (c_w1 == 0)
return return_binary64_zero(s);
k = clz64_nz(c_w1) - 10;
c_w1 = c_w1 << k;
}
}
// Correct to 2^112 <= c < 2^113 with corresponding exponent adding 113-54=59
// In fact shift a further 6 places ready for reciprocal multiplication
// Thus (113-54)+6=65, a shift of 1 given that we've already upacked in c.w[1]
c_w1 = c_w1 << 1;
c_w0 = 0;
k = k + 59;
// Check for "trivial" overflow, when 10^e * 1 > 2^{sci_emax+1}, just to
// keep tables smaller (it would be intercepted later otherwise).
//
// (Note that we may have normalized the coefficient, but we have a
// corresponding exponent postcorrection to account for; this can
// afford to be conservative anyway.)
//
// We actually check if e >= ceil((sci_emax + 1) * log_10(2))
// which in this case is 2 >= ceil(1024 * log_10(2)) = ceil(308.25) = 309
if (e >= 309) {
// __set_status_flags(pfpsf, BID_OVERFLOW_INEXACT_EXCEPTION);
return return_binary64_ovf(s, rnd_mode);
}
// Also check for "trivial" underflow, when 10^e * 2^113 <= 2^emin * 1/4,
// so test e <= floor((emin - 115) * log_10(2))
// In this case just fix ourselves at that value for uniformity.
//
// This is important not only to keep the tables small but to maintain the
// testing of the round/sticky words as a correct rounding method
if (e <= -358)
e = -358;
// Look up the breakpoint and approximate exponent
final int bid_breakpoints_binary64_index = (358 + e) << 1;
m_min_w0 = bid_breakpoints_binary64_BID_UINT128[bid_breakpoints_binary64_index];
m_min_w1 = bid_breakpoints_binary64_BID_UINT128[bid_breakpoints_binary64_index + 1];
e_out = bid_exponents_binary64[358 + e] - k;
// Choose provisional exponent and reciprocal multiplier based on breakpoint
final long[] bid_multipliers_binary64_BID_UINT256;
if (le128(c_w1, c_w0, m_min_w1, m_min_w0)) {
bid_multipliers_binary64_BID_UINT256 = bid_multipliers1_binary64_BID_UINT256;
} else {
bid_multipliers_binary64_BID_UINT256 = bid_multipliers2_binary64_BID_UINT256;
e_out = e_out + 1;
}
final int bid_multipliers_binary64_index = (358 + e) << 2;
r_w0 = bid_multipliers_binary64_BID_UINT256[bid_multipliers_binary64_index];
r_w1 = bid_multipliers_binary64_BID_UINT256[bid_multipliers_binary64_index + 1];
r_w2 = bid_multipliers_binary64_BID_UINT256[bid_multipliers_binary64_index + 2];
r_w3 = bid_multipliers_binary64_BID_UINT256[bid_multipliers_binary64_index + 3];
// Do the reciprocal multiplication
//__mul_64x256_to_320(z, c_w1, r);
{
long /*BID_UINT128*/ lP0_w0, lP0_w1, lP1_w0, lP1_w1, lP2_w0, lP2_w1, lP3_w0, lP3_w1;
long /*BID_UINT64*/ lC;
//__mul_64x64_to_128(lP0, c_w1, r_w0);
lP0_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w0);
lP0_w0 = c_w1 * r_w0;
//__mul_64x64_to_128(lP1, c_w1, r_w1);
lP1_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w1);
lP1_w0 = c_w1 * r_w1;
//__mul_64x64_to_128(lP2, c_w1, r_w2);
lP2_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w2);
lP2_w0 = c_w1 * r_w2;
//__mul_64x64_to_128(lP3, c_w1, r_w3);
lP3_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w3);
lP3_w0 = c_w1 * r_w3;
z_w0 = lP0_w0;
//__add_carry_out(P_w1,lC,lP1_w0,lP0_w1);
{
long /*BID_UINT64*/ X1 = lP1_w0;
z_w1 = lP1_w0 + lP0_w1;
lC = UnsignedLong.isLess(z_w1, X1) ? 1 : 0;
}
//__add_carry_in_out(P_w2,lC,lP2_w0,lP1_w1,lC);
{
long /*BID_UINT64*/ X1 = lP2_w0 + lC;
z_w2 = X1 + lP1_w1;
lC = (UnsignedLong.isLess(z_w2, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
//__add_carry_in_out(P_w3,lC,lP3_w0,lP2_w1,lC);
{
long /*BID_UINT64*/ X1 = lP3_w0 + lC;
z_w3 = X1 + lP2_w1;
lC = (UnsignedLong.isLess(z_w3, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
z_w4 = lP3_w1 + lC;
}
z_w5 = z_w4;
z_w4 = z_w3;
z_w3 = z_w2;
z_w2 = z_w1;
z_w1 = z_w0;
z_w0 = 0;
// Check for exponent underflow and compensate by shifting the product
// Cut off the process at precision+2, since we can't really shift further
if (e_out < 1) {
int d;
d = 1 - e_out;
if (d > 55)
d = 55;
e_out = 1;
//srl256_short(z_w5, z_w4, z_w3, z_w2, d);
{
z_w2 = (z_w3 << (64 - d)) + (z_w2 >>> d);
z_w3 = (z_w4 << (64 - d)) + (z_w3 >>> d);
z_w4 = (z_w5 << (64 - d)) + (z_w4 >>> d);
z_w5 = z_w5 >>> d;
}
}
c_prov = z_w5;
// Round using round-sticky words
// If we spill into the next binade, correct
// Flag underflow where it may be needed even for |result| = SNN
final int bid_roundbound_128_index = ((rnd_mode << 2) + ((s & 1) << 1) + (int) (c_prov & 1)) << 1;
if (lt128
(bid_roundbound_128_BID_UINT128[bid_roundbound_128_index + 1],
bid_roundbound_128_BID_UINT128[bid_roundbound_128_index], z_w4, z_w3)) {
c_prov = c_prov + 1;
if (c_prov == (1L << 53)) {
c_prov = 1L << 52;
e_out = e_out + 1;
}
// else if ((c_prov == (1L << 52)) && (e_out == 1)) {
// if (rnd_mode + (s & 1) == 2)
// __set_status_flags(pfpsf, BID_UNDERFLOW_EXCEPTION);
// }
}
// Check for overflow
if (e_out >= 2047) {
// __set_status_flags(pfpsf, BID_OVERFLOW_INEXACT_EXCEPTION);
return return_binary64_ovf(s, rnd_mode);
}
// Modify exponent for a tiny result, otherwise lop the implicit bit
if (UnsignedLong.isLess(c_prov, (1L << 52)))
e_out = 0;
else
c_prov = c_prov & ((1L << 52) - 1);
// Set the inexact and underflow flag as appropriate
// if (z_w4 != 0 || z_w3 != 0) {
// __set_status_flags(pfpsf, BID_INEXACT_EXCEPTION);
// if (e_out == 0)
// __set_status_flags(pfpsf, BID_UNDERFLOW_EXCEPTION);
// }
// Package up the result as a binary floating-point number
return return_binary64(s, e_out, c_prov);
}