in java/dfp/src/main/java/com/epam/deltix/dfp/JavaImplDiv.java [50:451]
public static long bid64_div(final long x, final long y) {
/*BID_UINT128*/
long CA_w0, CA_w1, CT_w0, CT_w1;
/*BID_UINT64*/
long A, B, QX, PD;
/*BID_UINT64*/
long A2, Q, Q2, B2, B4, B5, R, T, DU, res;
/*BID_SINT64*/
long D;
/*int_double*/
long t_scale_i, tempq_i, temp_b_i;
/*int_float*/
int tempx_i, tempy_i;
double da, db, dq, da_h, da_l;
int bin_expon_cx;
int diff_expon, ed1, ed2, bin_index;
int amount;
int nzeros, i, j, k, d5;
/*BID_UINT32*/
int QX32, digit, digit_h, digit_low;
int /*tdigit[3]*/ tdigit0, tdigit1, tdigit2;
//valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
final long sign_x;
int exponent_x;
long coefficient_x;
final long valid_x;
{
sign_x = x & 0x8000000000000000L;
if ((x & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
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
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 {
// coefficient
long coeff = (x & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_x = coeff;
// get exponent
long tmp = x >>> EXPONENT_SHIFT_LARGE64;
exponent_x = (int) (tmp & EXPONENT_MASK64);
valid_x = coeff;
}
}
}
//valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
long sign_y;
int exponent_y;
long coefficient_y;
final long valid_y;
{
sign_y = y & 0x8000000000000000L;
if ((y & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) {
// exponent
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
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 {
// coefficient
long coeff = (y & LARGE_COEFF_MASK64) | LARGE_COEFF_HIGH_BIT64;
// check for non-canonical values
if (coeff >= 10000000000000000L)
coeff = 0;
coefficient_y = coeff;
// get exponent
long 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 or NaN
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
// y==Inf, return NaN
if ((y & NAN_MASK64) == INFINITY_MASK64) { // Inf/Inf
return NAN_MASK64;
}
} else {
// otherwise return +/-Inf
return ((x ^ y) & 0x8000000000000000L) | INFINITY_MASK64;
}
}
// x==0
if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y == 0) {
// y==0 , return NaN
return NAN_MASK64;
}
if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
exponent_y = ((int) (y >>> 51)) & 0x3ff;
else
exponent_y = ((int) (y >>> 53)) & 0x3ff;
sign_y = y & 0x8000000000000000L;
exponent_x = 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;
return (sign_x ^ sign_y) | (((long) exponent_x) << 53);
}
}
if (valid_y == 0) {
// y is Inf. or NaN
// test if y is NaN
if ((y & NAN_MASK64) == NAN_MASK64) {
return (coefficient_y & QUIET_MASK64);
}
// y is Infinity?
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
// return +/-0
return ((x ^ y) & 0x8000000000000000L);
}
// y is 0
return (sign_x ^ sign_y) | INFINITY_MASK64;
}
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
if (UnsignedLong.isLess(coefficient_x, coefficient_y)) {
// get number of decimal digits for c_x, c_y
//--- get number of bits in the coefficients of x and y ---
tempx_i = Float.floatToRawIntBits(/*The UnsignedLong.floatValue is not required, because coefficient_x is always positive*/coefficient_x);
tempy_i = Float.floatToRawIntBits(/*The UnsignedLong.floatValue is not required, because coefficient_y is always positive*/coefficient_y);
bin_index = (tempy_i - tempx_i) >>> 23;
A = coefficient_x * bid_power10_index_binexp[bin_index];
B = coefficient_y;
temp_b_i = Double.doubleToRawLongBits(/*The UnsignedLong.doubleValue is not required, because coefficient_y is always positive*/B);
// compare A, B
DU = (A - B) >>> 63;
ed1 = 15 + (int) DU;
ed2 = bid_estimate_decimal_digits[bin_index] + ed1;
T = bid_power10_table_128_w0[ed1];
//__mul_64x64_to_128 (CA, A, T);
CA_w1 = Mul64Impl.unsignedMultiplyHigh(A, T);
CA_w0 = A * T;
Q = 0;
diff_expon = diff_expon - ed2;
// adjust double precision db, to ensure that later A/B - (int)(da/db) > -1
if (UnsignedLong.isLess(coefficient_y, 0x0020000000000000L)) {
temp_b_i += 1;
db = Double.longBitsToDouble(temp_b_i);
} else
db = UnsignedLong.doubleValue(B + 2 + (B & 1));
} else {
// get c_x/c_y
// set last bit before conversion to DP
A2 = coefficient_x | 1;
da = UnsignedLong.doubleValue(A2);
db = UnsignedLong.doubleValue(coefficient_y);
final double tempq_d = da / db;
tempq_i = Double.doubleToRawLongBits(tempq_d);
Q = UnsignedLong.fromDouble(tempq_d);
R = coefficient_x - coefficient_y * Q;
// will use to get number of dec. digits of Q
bin_expon_cx = (int) (tempq_i >>> 52) - 0x3ff;
// R<0 ?
D = R >> 63; // @AD: signed operation
Q += D;
R += (coefficient_y & D);
// exact result ?
if (R <= 0) { // @AD: signed operation
// can have R==-1 for coeff_y==1
return get_BID64(sign_x ^ sign_y, diff_expon, (Q + R));
}
// get decimal digits of Q
DU = bid_power10_index_binexp[bin_expon_cx] - Q - 1;
DU >>>= 63;
ed2 = 16 - bid_estimate_decimal_digits[bin_expon_cx] - (int) DU;
T = bid_power10_table_128_w0[ed2];
//__mul_64x64_to_128 (CA, R, T);
CA_w1 = Mul64Impl.unsignedMultiplyHigh(R, T);
CA_w0 = R * T;
B = coefficient_y;
Q *= bid_power10_table_128_w0[ed2];
diff_expon -= ed2;
}
if (CA_w1 == 0) {
Q2 = UnsignedLong.divide(CA_w0, B);
B2 = B + B;
B4 = B2 + B2;
R = CA_w0 - Q2 * B;
Q += Q2;
} else {
// 2^64
t_scale_i = 0x43f0000000000000L;
// convert CA to DP
da_h = UnsignedLong.doubleValue(CA_w1);
da_l = UnsignedLong.doubleValue(CA_w0);
da = da_h * Double.longBitsToDouble(t_scale_i) + da_l;
// quotient
dq = da / db;
Q2 = UnsignedLong.fromDouble(dq);
// get w[0] remainder
R = CA_w0 - Q2 * B;
// R<0 ?
D = R >> 63; // @AD: signed operation
Q2 += D;
R += (B & D);
// now R<6*B
// quick divide
// 4*B
B2 = B + B;
B4 = B2 + B2;
R = R - B4;
// R<0 ?
D = R >> 63; // @AD: signed operation
// restore R if negative
R += (B4 & D);
Q2 += ((~D) & 4);
R = R - B2;
// R<0 ?
D = R >> 63; // @AD: signed operation
// restore R if negative
R += (B2 & D);
Q2 += ((~D) & 2);
R = R - B;
// R<0 ?
D = R >> 63; // @AD: signed operation
// restore R if negative
R += (B & D);
Q2 += ((~D) & 1);
Q += Q2;
}
if (R == 0) {
// eliminate trailing zeros
// check whether CX, CY are short
if ((UnsignedLong.isLessOrEqual(coefficient_x, 1024)) &&
(UnsignedLong.isLessOrEqual(coefficient_y, 1024))) {
i = (int) coefficient_y - 1;
j = (int) coefficient_x - 1;
// difference in powers of 2 bid_factors for Y and X
nzeros = ed2 - bid_factors_flat[i << 1] + bid_factors_flat[j << 1];
// difference in powers of 5 bid_factors
d5 = ed2 - bid_factors_flat[(i << 1) + 1] + bid_factors_flat[(j << 1) + 1];
if (d5 < nzeros)
nzeros = d5;
//__mul_64x64_to_128 (CT, Q, bid_reciprocals10_64[nzeros]);
CT_w1 = Mul64Impl.unsignedMultiplyHigh(Q, bid_reciprocals10_64[nzeros]);
// CT_w0 = Q * bid_reciprocals10_64[nzeros]; // @optimization
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[nzeros];
Q = CT_w1 >>> amount;
diff_expon += nzeros;
} else {
tdigit0 = (int) (Q & 0x3ffffff);
tdigit1 = 0;
QX = Q >>> 26;
QX32 = (int) QX;
nzeros = 0;
for (j = 0; QX32 != 0; j++, QX32 >>>= 7) {
k = (QX32 & 127);
tdigit0 += bid_convert_table[j][k][0];
tdigit1 += bid_convert_table[j][k][1];
if (/*UnsignedInteger.compare*/(tdigit0) + Integer.MIN_VALUE >= (100000000) + Integer.MIN_VALUE) {
tdigit0 -= 100000000;
tdigit1++;
}
}
digit = tdigit0;
if (digit == 0 && tdigit1 == 0)
nzeros += 16;
else {
if (digit == 0) {
nzeros += 8;
digit = tdigit1;
}
// decompose digit
PD = (long) (digit & LONG_LOW_PART) * 0x068DB8BBL;
digit_h = (int) (PD >>> 40);
digit_low = digit - digit_h * 10000;
if (digit_low == 0)
nzeros += 4;
else
digit_h = digit_low;
if ((digit_h & 1) == 0)
nzeros += 3 & (int) (bid_packed_10000_zeros[digit_h >>> 3] >>> (digit_h & 7));
}
if (nzeros != 0) {
//__mul_64x64_to_128 (CT, Q, bid_reciprocals10_64[nzeros]);
CT_w1 = Mul64Impl.unsignedMultiplyHigh(Q, bid_reciprocals10_64[nzeros]);
// CT_w0 = Q * bid_reciprocals10_64[nzeros]; // @optimization
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
amount = bid_short_recip_scale[nzeros];
Q = CT_w1 >>> amount;
}
diff_expon += nzeros;
}
if (diff_expon >= 0) {
return fast_get_BID64_check_OF(sign_x ^ sign_y, diff_expon, Q);
}
}
if (diff_expon >= 0) {
// R*10
R += R;
R = (R << 2) + R;
B5 = B4 + B;
// compare 10*R to 5*B
R = B5 - R;
// correction for (R==0 && (Q&1))
R -= Q & 1;
// R<0 ?
D = R >>> 63;
Q += D;
return fast_get_BID64_check_OF(sign_x ^ sign_y, diff_expon, Q);
} else {
// UF occurs
return get_BID64_UF(sign_x ^ sign_y, diff_expon, Q, R);
}
}