in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplCastBinary64.java [15:308]
public static long binary64_to_bid64(double x, final int rnd_mode/*, final JavaImplParse.FloatingPointStatusFlag pfpsf*/) {
long /*BID_UINT128*/ c_w0, c_w1;
long /*BID_UINT64*/ c_prov;
long /*BID_UINT128*/ m_min_w0, m_min_w1;
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;
int e, s, t, e_out;
// Unpack the input
//unpack_binary64 (x, s, e, c.w[1], t, return_bid64_zero (s), return_bid64_inf (s), return_bid64_nan);
{
//union { BID_UINT64 i; double f; } x_in;
//x_in.f = x;
//c = x_in.i;
c_w1 = Double.doubleToRawLongBits(x);
e = (int) ((c_w1 >>> 52) & ((1L << 11) - 1));
s = (int) (c_w1 >>> 63);
c_w1 = c_w1 & ((1L << 52) - 1);
if (e == 0) {
int l;
if (c_w1 == 0) return return_bid64_zero(s);
l = clz64(c_w1) - (64 - 53);
c_w1 = c_w1 << l;
e = -(l + 1074);
t = 0;
//__set_status_flags(pfpsf, BID_DENORMAL_EXCEPTION);
} else if (e == ((1L << 11) - 1)) {
if (c_w1 == 0) return return_bid64_inf(s);
//if ((c_w1 & (1L << 51)) == 0) __set_status_flags(pfpsf, BID_INVALID_EXCEPTION);
return return_bid64_nan(s, c_w1 << 13, 0L);
} else {
c_w1 += 1L << 52;
t = ctz64(c_w1);
e -= 1075;
}
}
// Now -1126<=e<=971 (971 for max normal, -1074 for min normal, -1126 for min denormal)
// Treat like a quad input for uniformity, so (2^{113-53} * c * r) >> 312
// (312 is the shift value for these tables) which can be written as
// (2^68 c * r) >> 320, lopping off exactly 320 bits = 5 words. Thus we put
// input coefficient as the high part of c (<<64) shifted by 4 bits (<<68)
//
// Remember to compensate for the fact that exponents are integer for quad
c_w1 = c_w1 << 4;
c_w0 = 0;
t += (113 - 53);
e -= (113 - 53); // Now e belongs [-1186;911].
// Check for "trivial" overflow, when 2^e * 2^112 > 10^emax * 10^d.
// We actually check if e >= ceil((emax + d) * log_2(10) - 112)
// This could be intercepted later, but it's convenient to keep tables smaller
// Now filter out all the exact cases where we need to specially force
// the exponent to 0. We can let through inexact cases and those where the
// main path will do the right thing anyway, e.g. integers outside coeff range.
//
// First check that e <= 0, because if e > 0, the input must be >= 2^113,
// which is too large for the coefficient of any target decimal format.
// We write a = -(e + t)
//
// (1) If e + t >= 0 <=> a <= 0 the input is an integer; treat it specially
// iff it fits in the coefficient range. Shift c' = c >> -e, and
// compare with the coefficient range; if it's in range then c' is
// our coefficient, exponent is 0. Otherwise we pass through.
//
// (2) If a > 0 then we have a non-integer input. The special case would
// arise as c' / 2^a where c' = c >> t, i.e. 10^-a * (5^a c'). Now
// if a > 48 we can immediately forget this, since 5^49 > 10^34.
// Otherwise we determine whether we're in range by a table based on
// a, and if so get the multiplier also from a table based on a.
//
// Note that when we shift, we need to take into account the fact that
// c is already 8 places to the left in preparation for the reciprocal
// multiplication; thus we add 8 to all the shift counts
if (e <= 0) {
long /*BID_UINT128*/ cint_w1 = c_w1, cint_w0 = c_w0;
int a = -(e + t);
if (a <= 0) {
//srl128 (cint_w1, cint_w0, 8 - e);
{
final int __c = 8 - e;
// if (__c != 0) // No need for condition check: it is always true
{
if (__c >= 64) {
cint_w0 = cint_w1 >>> (__c - 64);
cint_w1 = 0;
} else {
//srl128_short(cint_w1, cint_w0, __c);
cint_w0 = (cint_w1 << (64 - __c)) + (cint_w0 >>> __c);
cint_w1 = cint_w1 >>> __c;
}
}
}
if ((cint_w1 == 0) && UnsignedLong.isLess(cint_w0, 10000000000000000L))
return return_bid64(s, 398, cint_w0);
} else if (a <= 48) {
final int bid_coefflimits_bid64_index = a << 1;
long /*BID_UINT128*/ pow5_w0 = bid_coefflimits_bid64_BID_UINT128[bid_coefflimits_bid64_index],
pow5_w1 = bid_coefflimits_bid64_BID_UINT128[bid_coefflimits_bid64_index + 1];
//srl128 (cint_w1, cint_w0, 8 + t);
{
final int __c = 8 + t;
if (__c != 0) {
if (__c >= 64) {
cint_w0 = cint_w1 >>> (__c - 64);
cint_w1 = 0;
} else {
//srl128_short(cint_w1, cint_w0, __c);
cint_w0 = (cint_w1 << (64 - __c)) + (cint_w0 >>> __c);
cint_w1 = cint_w1 >>> __c;
}
}
}
if (le128(cint_w1, cint_w0, pow5_w1, pow5_w0)) {
long /*BID_UINT128*/ cc_w1 = cint_w1, cc_w0 = cint_w0;
final int bid_power_five_index = a << 1;
pow5_w0 = bid_power_five_BID_UINT128[bid_power_five_index];
pow5_w1 = bid_power_five_BID_UINT128[bid_power_five_index + 1];
//__mul_128x128_low (cc, cc, pow5);
{
long /*BID_UINT128*/ ALBL_w0, ALBL_w1;
long /*BID_UINT64*/ QM64;
//__mul_64x64_to_128(ALBL, cc_w0, pow5_w0);
ALBL_w1 = Mul64Impl.unsignedMultiplyHigh(cc_w0, pow5_w0);
ALBL_w0 = cc_w0 * pow5_w0;
QM64 = pow5_w0 * cc_w1 + cc_w0 * pow5_w1;
cc_w0 = ALBL_w0;
cc_w1 = QM64 + ALBL_w1;
}
return return_bid64(s, 398 - a, cc_w0);
}
}
}
// Check for "trivial" underflow, when 2^e * 2^113 <= 10^emin * 1/4,
// so test e <= floor(emin * log_2(10) - 115)
// 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
// Now look up our exponent e, and the breakpoint between e and e+1
final int bid_breakpoints_bid64_index = (1437 + e) << 1;
m_min_w0 = bid_breakpoints_bid64_BID_UINT128[bid_breakpoints_bid64_index];
m_min_w1 = bid_breakpoints_bid64_BID_UINT128[bid_breakpoints_bid64_index + 1];
e_out = bid_exponents_bid64[1437 + e];
// Choose exponent and reciprocal multiplier based on breakpoint
final long[] bid_multipliers_bid64_BID_UINT256;
if (le128(c_w1, c_w0, m_min_w1, m_min_w0)) {
bid_multipliers_bid64_BID_UINT256 = bid_multipliers1_bid64_BID_UINT256;
} else {
bid_multipliers_bid64_BID_UINT256 = bid_multipliers2_bid64_BID_UINT256;
e_out = e_out + 1;
}
final int bid_multipliers_bid64_index = (1437 + e) << 2;
r_w0 = bid_multipliers_bid64_BID_UINT256[bid_multipliers_bid64_index];
r_w1 = bid_multipliers_bid64_BID_UINT256[bid_multipliers_bid64_index + 1];
r_w2 = bid_multipliers_bid64_BID_UINT256[bid_multipliers_bid64_index + 2];
r_w3 = bid_multipliers_bid64_BID_UINT256[bid_multipliers_bid64_index + 3];
// Do the reciprocal multiplication
//__mul_128x256_to_384 (z, c, r)
{
long /*BID_UINT512*/ P0_w0, P0_w1, P0_w2, P0_w3, P0_w4, P0_w5, P0_w6, P0_w7,
P1_w0, P1_w1, P1_w2, P1_w3, P1_w4, P1_w5, P1_w6, P1_w7;
long /*BID_UINT64*/ CY;
{
long /*BID_UINT128*/ lP0_w0, lP0_w1, lP1_w0, lP1_w1, lP2_w0, lP2_w1, lP3_w0, lP3_w1;
long /*BID_UINT64*/ lC;
lP0_w1 = Mul64Impl.unsignedMultiplyHigh(c_w0, r_w0);
lP0_w0 = c_w0 * r_w0;
lP1_w1 = Mul64Impl.unsignedMultiplyHigh(c_w0, r_w1);
lP1_w0 = c_w0 * r_w1;
lP2_w1 = Mul64Impl.unsignedMultiplyHigh(c_w0, r_w2);
lP2_w0 = c_w0 * r_w2;
lP3_w1 = Mul64Impl.unsignedMultiplyHigh(c_w0, r_w3);
lP3_w0 = c_w0 * r_w3;
P0_w0 = lP0_w0;
{
long /*BID_UINT64*/ X1 = lP1_w0;
P0_w1 = lP1_w0 + lP0_w1;
lC = UnsignedLong.isLess(P0_w1, X1) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = lP2_w0 + lC;
P0_w2 = X1 + lP1_w1;
lC = (UnsignedLong.isLess(P0_w2, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = lP3_w0 + lC;
P0_w3 = X1 + lP2_w1;
lC = (UnsignedLong.isLess(P0_w3, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
P0_w4 = lP3_w1 + lC;
}
{
long /*BID_UINT128*/ lP0_w0, lP0_w1, lP1_w0, lP1_w1, lP2_w0, lP2_w1, lP3_w0, lP3_w1;
long /*BID_UINT64*/ lC;
lP0_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w0);
lP0_w0 = c_w1 * r_w0;
lP1_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w1);
lP1_w0 = c_w1 * r_w1;
lP2_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w2);
lP2_w0 = c_w1 * r_w2;
lP3_w1 = Mul64Impl.unsignedMultiplyHigh(c_w1, r_w3);
lP3_w0 = c_w1 * r_w3;
P1_w0 = lP0_w0;
{
long /*BID_UINT64*/ X1 = lP1_w0;
P1_w1 = lP1_w0 + lP0_w1;
lC = UnsignedLong.isLess(P1_w1, X1) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = lP2_w0 + lC;
P1_w2 = X1 + lP1_w1;
lC = (UnsignedLong.isLess(P1_w2, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = lP3_w0 + lC;
P1_w3 = X1 + lP2_w1;
lC = (UnsignedLong.isLess(P1_w3, X1) || UnsignedLong.isLess(X1, lC)) ? 1 : 0;
}
P1_w4 = lP3_w1 + lC;
}
z_w0 = P0_w0;
{
long /*BID_UINT64*/ X1 = P1_w0;
z_w1 = P1_w0 + P0_w1;
CY = UnsignedLong.isLess(z_w1, X1) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = P1_w1 + CY;
z_w2 = X1 + P0_w2;
CY = (UnsignedLong.isLess(z_w2, X1) || UnsignedLong.isLess(X1, CY)) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = P1_w2 + CY;
z_w3 = X1 + P0_w3;
CY = (UnsignedLong.isLess(z_w3, X1) || UnsignedLong.isLess(X1, CY)) ? 1 : 0;
}
{
long /*BID_UINT64*/ X1 = P1_w3 + CY;
z_w4 = X1 + P0_w4;
CY = (UnsignedLong.isLess(z_w4, X1) || UnsignedLong.isLess(X1, CY)) ? 1 : 0;
}
z_w5 = P1_w4 + CY;
}
c_prov = z_w5;
// Round using round-sticky words
// If we spill over into the next decade, correct
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 == 10000000000000000L) {
c_prov = 1000000000000000L;
e_out = e_out + 1;
}
}
// Check for overflow
// Set the inexact flag as appropriate and check underflow
// It's no doubt superfluous to check inexactness, but anyway...
// if ((z_w4 != 0) || (z_w3 != 0)) {
// __set_status_flags(pfpsf, BID_INEXACT_EXCEPTION);
// }
// Package up the result
return return_bid64(s, e_out, c_prov);
}