R128.w[1] = R256.w[1];
R128.w[0] = R256.w[0];
}
+ if (e4 + x0 < expmin) { // for all rounding modes
+ is_tiny = 1;
+ }
// the rounded result has p34 = 34 digits
e4 = e4 + x0 + incr_exp;
- if (rnd_mode == ROUNDING_TO_NEAREST) {
- if (e4 < expmin) {
- is_tiny = 1; // for other rounding modes apply correction
- }
- } else {
+ if (rnd_mode != ROUNDING_TO_NEAREST) {
// for RM, RP, RZ, RA apply correction in order to determine tininess
// but do not save the result; apply the correction to
// (-1)^p_sign * significand * 10^0
is_inexact_gt_midpoint, is_midpoint_lt_even,
is_midpoint_gt_even, 0, &P128, ptrfpsf);
scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
- // the number of digits in the significand is p34 = 34
- if (e4 + scale < expmin) {
- is_tiny = 1;
- }
}
ind = p34; // the number of decimal digits in the signifcand of res
res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R128.w[1]; // RN
}
}
}
-
p_sign = x_sign ^ y_sign; // sign of the product
// identify cases where at least one operand is infinity
if (C1.w[1] == 0) {
if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
// split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
+
tmp.d = (double) (C1.w[0] >> 32); // exact conversion
x_nr_bits =
33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
} else { // if x < 2^53
tmp.d = (double) C1.w[0]; // exact conversion
x_nr_bits =
if (q1 == 0) {
q1 = nr_digits[x_nr_bits - 1].digits1;
if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
+ (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
+ C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
q1++;
}
}
-
+ // q2 = nr. of decimal digits in y
+ // determine first the nr. of bits in y
if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f (non-zero finite)
if (C2.w[1] == 0) {
if (C2.w[0] >= 0x0020000000000000ull) { // y >= 2^53
// split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C2.w[0] >= 0x0000000100000000ull) { // y >= 2^32
tmp.d = (double) (C2.w[0] >> 32); // exact conversion
y_nr_bits =
- 32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // y < 2^32
- tmp.d = (double) C2.w[0]; // exact conversion
- y_nr_bits =
- ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
+ 33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
} else { // if y < 2^53
tmp.d = (double) C2.w[0]; // exact conversion
y_nr_bits =
- ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
} else { // C2.w[1] != 0 => nr. bits = 64 + nr_bits (C2.w[1])
tmp.d = (double) C2.w[1]; // exact conversion
y_nr_bits =
- 64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ 65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
-
- q2 = nr_digits[y_nr_bits].digits;
+ q2 = nr_digits[y_nr_bits - 1].digits;
if (q2 == 0) {
- q2 = nr_digits[y_nr_bits].digits1;
- if (C2.w[1] > nr_digits[y_nr_bits].threshold_hi ||
- (C2.w[1] == nr_digits[y_nr_bits].threshold_hi &&
- C2.w[0] >= nr_digits[y_nr_bits].threshold_lo))
+ q2 = nr_digits[y_nr_bits - 1].digits1;
+ if (C2.w[1] > nr_digits[y_nr_bits - 1].threshold_hi ||
+ (C2.w[1] == nr_digits[y_nr_bits - 1].threshold_hi &&
+ C2.w[0] >= nr_digits[y_nr_bits - 1].threshold_lo))
q2++;
}
}
if (C3.w[1] == 0) {
if (C3.w[0] >= 0x0020000000000000ull) { // z >= 2^53
// split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C3.w[0] >= 0x0000000100000000ull) { // z >= 2^32
tmp.d = (double) (C3.w[0] >> 32); // exact conversion
z_nr_bits =
- 32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // z < 2^32
- tmp.d = (double) C3.w[0]; // exact conversion
- z_nr_bits =
- ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
+ 33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
} else { // if z < 2^53
tmp.d = (double) C3.w[0]; // exact conversion
z_nr_bits =
- ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
} else { // C3.w[1] != 0 => nr. bits = 64 + nr_bits (C3.w[1])
tmp.d = (double) C3.w[1]; // exact conversion
z_nr_bits =
- 64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ 65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
-
- q3 = nr_digits[z_nr_bits].digits;
+ q3 = nr_digits[z_nr_bits - 1].digits;
if (q3 == 0) {
- q3 = nr_digits[z_nr_bits].digits1;
- if (C3.w[1] > nr_digits[z_nr_bits].threshold_hi ||
- (C3.w[1] == nr_digits[z_nr_bits].threshold_hi &&
- C3.w[0] >= nr_digits[z_nr_bits].threshold_lo))
+ q3 = nr_digits[z_nr_bits - 1].digits1;
+ if (C3.w[1] > nr_digits[z_nr_bits - 1].threshold_hi ||
+ (C3.w[1] == nr_digits[z_nr_bits - 1].threshold_hi &&
+ C3.w[0] >= nr_digits[z_nr_bits - 1].threshold_lo))
q3++;
}
}
} else {
; // continue with x = f, y = f, z = 0 or x = f, y = f, z = f
}
-
e1 = (x_exp >> 49) - 6176; // unbiased exponent of x
e2 = (y_exp >> 49) - 6176; // unbiased exponent of y
e3 = (z_exp >> 49) - 6176; // unbiased exponent of z
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
q4 = q1 + q2; // q4 in [40, 57]
}
- } else if (q1 + q2 == 58) { // C4 = C1 * C2 fits in 192 or 256 bits
- // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
- // may fit in 64 bits
- if (C1.w[1] == 0) { // C1 * C2 will fit in 192 bits
- __mul_64x128_full (C4.w[2], C4, C1.w[0], C2); // may use 64x128_to_192
- } else if (C2.w[1] == 0) { // C1 * C2 will fit in 192 bits
- __mul_64x128_full (C4.w[2], C4, C2.w[0], C1); // may use 64x128_to_192
- } else { // C1 * C2 will fit in 192 bits or in 256 bits
- __mul_128x128_to_256 (C4, C1, C2);
- }
+ } else if (q1 + q2 == 58) { // C4 = C1 * C2 fits in 192 or 256 bits;
+ // both C1 and C2 fit in 128 bits (actually in 113 bits); none can
+ // fit in 64 bits, because each number must have at least 24 decimal
+ // digits for the sum to have 58 (as the max. nr. of digits is 34) =>
+ // C1.w[1] != 0 and C2.w[1] != 0
+ __mul_128x128_to_256 (C4, C1, C2);
// if C4 < 10^(q1+q2-1) = 10^57 then q4 = q1+q2-1 = 57 else q4 = q1+q2 = 58
if (C4.w[3] == 0 && (C4.w[2] < ten2k256[18].w[2] ||
- (C4.w[2] == ten2k256[18].w[2]
- && (C4.w[1] < ten2k256[18].w[1]
- || (C4.w[1] == ten2k256[18].w[1]
- && C4.w[0] < ten2k256[18].w[0]))))) {
+ (C4.w[2] == ten2k256[18].w[2]
+ && (C4.w[1] < ten2k256[18].w[1]
+ || (C4.w[1] == ten2k256[18].w[1]
+ && C4.w[0] < ten2k256[18].w[0]))))) {
// 18 = 57 - 39 = q1+q2-1 - 39
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
q4 = 57; // 57 = q1 + q2 - 1
q4 = q1 + q2; // q4 in [59, 68]
}
}
-
if (C3.w[1] == 0x0 && C3.w[0] == 0x0) { // x = f, y = f, z = 0
save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
*pfpsf = 0;
res.w[1] = R256.w[1];
}
e4 = e4 + x0;
+ q4 = p34;
if (incr_exp) {
e4 = e4 + 1;
+ if (q4 + e4 == expmin + p34) *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
}
- q4 = p34;
// res is now the coefficient of the result rounded to the destination
// precision, with unbounded exponent; the exponent is e4; q4=digits(res)
} else { // if (q4 <= p34)
delta = q3 + e3 - q4 - e4;
delta_ge_zero:
if (delta >= 0) {
-
if (p34 <= delta - 1 || // Case (1')
(p34 == delta && e3 + 6176 < p34 - q3)) { // Case (1''A)
// check for overflow, which can occur only in Case (1')
res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | C3.w[1];
res.w[0] = C3.w[0];
}
-
+
// use the following to avoid double rounding errors when operating on
// mixed formats in rounding to nearest, and for correcting the result
// if not rounding to nearest
R64 = 10;
}
}
- if (q4 == 1 && C4.w[0] == 5) {
+
+ if (R64 == 5 && !is_inexact_lt_midpoint && !is_inexact_gt_midpoint &&
+ !is_midpoint_lt_even && !is_midpoint_gt_even) {
+ //if (q4 == 1 && C4.w[0] == 5) {
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 1;
res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
}
if (e3 == expmin) {
- if (R64 < 5 || (R64 == 5 && !is_inexact_lt_midpoint)) {
- ; // result not tiny (in round-to-nearest mode)
- } else {
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
+ *pfpsf |= UNDERFLOW_EXCEPTION; // tiny if detected before rounding
}
} // end 10^(q3+scale-1)
// set the inexact flag
// endif
if ((e3 == expmin && (q3 + scale) < p34) ||
(e3 == expmin && (q3 + scale) == p34 &&
- (res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull && // 10^33_high
- res.w[0] == 0x38c15b0a00000000ull && // 10^33_low
- z_sign != p_sign && ((!z_sign && rnd_mode != ROUNDING_UP) ||
- (z_sign && rnd_mode != ROUNDING_DOWN)))) {
+ (res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull && // 10^33_high
+ res.w[0] == 0x38c15b0a00000000ull && // 10^33_low
+ z_sign != p_sign)) {
*pfpsf |= UNDERFLOW_EXCEPTION;
}
if (rnd_mode != ROUNDING_TO_NEAREST) {
if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull ||
(res.w[1] == 0x0000314dc6448d93ull &&
res.w[0] < 0x38c15b0a00000000ull)) ||
- (is_inexact_lt_midpoint
+ ((is_inexact_lt_midpoint | is_midpoint_gt_even)
&& res.w[1] == 0x0000314dc6448d93ull
&& res.w[0] == 0x38c15b0a00000000ull))
&& x0 >= 1) {
res.w[0] < 0x38c15b0a00000000ull)) {
is_tiny = 1;
}
+ if (((res.w[1] & 0x7fffffffffffffffull) == 0x0000314dc6448d93ull) &&
+ (res.w[0] == 0x38c15b0a00000000ull) && // 10^33*10^-6176
+ (z_sign != p_sign)) is_tiny = 1;
} else if (e3 < expmin) {
// the result is tiny, so we must truncate more of res
is_tiny = 1;
0, &P128, pfpsf);
scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
// the number of digits in the significand is p34 = 34
- if (e4 + scale < expmin) {
- is_tiny = 1;
- }
}
// the result rounded to the destination precision with unbounded exponent
is_midpoint_lt_even, is_midpoint_gt_even,
e4, &res, pfpsf);
}
+ // correction needed for tininess detection before rounding
+ if ((((res.w[1] & 0x7fffffffffffffffull) == 0x0000314dc6448d93ull) &&
+ // 10^33*10^-6176_high
+ (res.w[0] == 0x38c15b0a00000000ull)) && // 10^33*10^-6176_low
+ (((rnd_mode == ROUNDING_TO_NEAREST ||
+ rnd_mode == ROUNDING_TIES_AWAY) &&
+ (is_midpoint_lt_even || is_inexact_gt_midpoint)) ||
+ ((((rnd_mode == ROUNDING_UP) && !(res.w[1] & MASK_SIGN)) ||
+ ((rnd_mode == ROUNDING_DOWN) && (res.w[1] & MASK_SIGN)))
+ && (is_midpoint_lt_even || is_midpoint_gt_even ||
+ is_inexact_lt_midpoint || is_inexact_gt_midpoint)))) {
+ is_tiny = 1;
+ }
if (is_midpoint_lt_even || is_midpoint_gt_even ||
is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
// set the inexact flag
// determine the unbiased exponent of the result
unbexp = ((res1 >> 53) & 0x3ff) - 398;
+ if (!((res1 & MASK_NAN) == MASK_NAN)) { // res1 not NaN
// if subnormal, res1 must have exp = -398
// if tiny and inexact set underflow and inexact status flags
- if (!((res1 & MASK_NAN) == MASK_NAN) && // res1 not NaN
- (unbexp == -398)
- && ((res1 & MASK_BINARY_SIG1) < 1000000000000000ull)
- && (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0
- || is_midpoint_lt_even0 || is_midpoint_gt_even0)) {
- // set the inexact flag and the underflow flag
- *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
+ if ((unbexp == -398)
+ && ((res1 & MASK_BINARY_SIG1) < 1000000000000000ull)
+ && (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0
+ || is_midpoint_lt_even0 || is_midpoint_gt_even0)) {
+ // set the inexact flag and the underflow flag
+ *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
} else if (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 ||
is_midpoint_lt_even0 || is_midpoint_gt_even0) {
// set the inexact flag and the underflow flag
*pfpsf |= INEXACT_EXCEPTION;
- }
+ }
+ // correction needed for tininess detection before rounding
+ if (((res1 & 0x7fffffffffffffffull) == 1000000000000000ull) &&
+ // 10^15*10^-398
+ (((rnd_mode == ROUNDING_TO_NEAREST ||
+ rnd_mode == ROUNDING_TIES_AWAY) &&
+ (is_midpoint_lt_even || is_inexact_gt_midpoint)) ||
+ ((((rnd_mode == ROUNDING_UP) && !(res1 & MASK_SIGN)) ||
+ ((rnd_mode == ROUNDING_DOWN) && (res1 & MASK_SIGN)))
+ && (is_midpoint_lt_even || is_midpoint_gt_even ||
+ is_inexact_lt_midpoint || is_inexact_gt_midpoint)))) {
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+ }
*pfpsf |= save_fpsf;
BID_RETURN (res1);
} // else continue, and use rounding to nearest to round to 16 digits
res1 = sign | MASK_STEERING_BITS |
((UINT64) (unbexp + 398) << 51) | (res1 & MASK_BINARY_SIG2);
}
+
+ // correction needed for tininess detection before rounding
+ if (((res1 & 0x7fffffffffffffffull) == 1000000000000000ull) &&
+ // 10^15*10^-398
+ (((rnd_mode == ROUNDING_TO_NEAREST ||
+ rnd_mode == ROUNDING_TIES_AWAY) &&
+ (is_midpoint_lt_even || is_inexact_gt_midpoint)) ||
+ ((((rnd_mode == ROUNDING_UP) && !(res1 & MASK_SIGN)) ||
+ ((rnd_mode == ROUNDING_DOWN) && (res1 & MASK_SIGN)))
+ && (is_midpoint_lt_even || is_midpoint_gt_even ||
+ is_inexact_lt_midpoint || is_inexact_gt_midpoint)))) {
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
*pfpsf |= save_fpsf;
BID_RETURN (res1);
}
bid128_class (int *pres, UINT128 * px _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
UINT128 x = *px;
#else
-int
+class_t
bid128_class (UINT128 x _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
int res;
BID_RETURN (res);
}
break;
+default: break; // default added to avoid compiler warning
}
// q = nr. of decimal digits in x
BID_RETURN (res);
}
break;
+default: break; // default added to avoid compiler warning
}
BID_RETURN (res);
UINT128 C1;
unsigned int k = 0; // pointer in the string
unsigned int d0, d123;
+ unsigned int zero_digit = (unsigned int) '0';
UINT64 HI_18Dig, LO_18Dig, Tmp;
UINT32 MiDi[12], *ptr;
char *c_ptr_start, *c_ptr;
d123 = exp - 1000 * d0;
if (d0) { // 1000 <= exp <= 6144 => 4 digits to return
- str[k++] = d0 + 0x30;// ASCII for decimal digit d0
+ str[k++] = d0 + zero_digit; // ASCII for decimal digit d0
ind = 3 * d123;
str[k++] = char_table3[ind];
str[k++] = char_table3[ind + 1];
str[k++] = char_table3[ind + 2];
} else { // 0 <= exp <= 999 => d0 = 0
if (d123 < 10) { // 0 <= exp <= 9 => 1 digit to return
- str[k++] = d123 + 0x30;// ASCII
+ str[k++] = d123 + zero_digit; // ASCII
} else if (d123 < 100) { // 10 <= exp <= 99 => 2 digits to return
ind = 2 * (d123 - 10);
str[k++] = char_table2[ind];
}
break;
-
+ default: break; // default added to avoid compiler warning
}
// now form the coefficient as coeff_high*10^17+coeff_low+carry
scale_high = 100000000000000000ull;
T128 = round_const_table_128[rmode][extra_digits];
__add_carry_out (CX1.w[0], carry, T128.w[0], CX.w[0]);
CX1.w[1] = CX.w[1] + T128.w[1] + carry;
- if (__unsigned_compare_ge_128
- (CX1, power10_table_128[extra_digits + 7]))
- uf_check = 0;
}
extra_digits =
extra_digits + DECIMAL_EXPONENT_BIAS_128 -
UINT128 Q;
UINT64 sign_x, coefficient_x, remainder_h, carry, Stemp;
UINT32 res;
+ UINT64 t64;
int_float tempx;
int exponent_x, bin_expon_cx, extra_digits, rmode = 0, amount;
unsigned status = 0;
// unpack arguments, check for NaN or Infinity, 0
if (!unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x)) {
if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
- res = (coefficient_x & 0x0003ffffffffffffull);
- res /= 1000000000ull;
+ t64 = (coefficient_x & 0x0003ffffffffffffull);
+ res = t64/1000000000ull;
+ //res = (coefficient_x & 0x0003ffffffffffffull);
+ //res /= 1000000000ull;
res |= ((coefficient_x >> 32) & 0xfc000000);
#ifdef SET_STATUS_FLAGS
if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
exponent_x += extra_digits;
if ((exponent_x < 0) && (exponent_x + MAX_FORMAT_DIGITS_32 >= 0)) {
status = UNDERFLOW_EXCEPTION;
- if (exponent_x == -1)
- if (coefficient_x + round_const_table[rmode][extra_digits] >=
- power10_table_128[extra_digits + 7].w[0])
- status = 0;
extra_digits -= exponent_x;
exponent_x = 0;
}
bid64_class (int *pres, UINT64 * px _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
UINT64 x = *px;
#else
-int
+class_t
bid64_class (UINT64 x _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
int res;
BID_RETURN (res);
}
break;
+ default: break; // default added to avoid compiler warning
} // end switch ()
// q = nr. of decimal digits in x (1 <= q <= 54)
BID_RETURN (res);
}
break;
+ default: break; // default added to avoid compiler warning
} // end switch ()
BID_RETURN (res);
}
#endif
UINT64 sign_x, coefficient_x = 0, rounded = 0, res;
int expon_x = 0, sgn_expon, ndigits, add_expon = 0, midpoint =
- 0, rounded_up = 0;
+ 0, rounded_up = 0, dround = 0;
int dec_expon_scale = 0, right_radix_leading_zeros = 0, rdx_pt_enc =
0;
unsigned fpsc;
break;
case ROUNDING_DOWN:
- if(sign_x) { coefficient_x++; rounded_up=1; }
+ if(sign_x) { if(c>'0') {coefficient_x++; rounded_up=1;} else dround=1; }
break;
case ROUNDING_UP:
- if(!sign_x) { coefficient_x++; rounded_up=1; }
+ if(!sign_x) { if(c>'0') {coefficient_x++; rounded_up=1;} else dround=1; }
break;
case ROUNDING_TIES_AWAY:
if(c>='5') { coefficient_x++; rounded_up=1; }
midpoint = 0;
rounded_up = 1;
}
- if (c > '0')
+ if (c > '0') {
rounded = 1;
+
+ if(dround)
+ {
+ dround = 0;
+ coefficient_x ++;
+ rounded_up = 1;
+
+ if (coefficient_x == 10000000000000000ull) {
+ coefficient_x = 1000000000000000ull;
+ add_expon = 1;
+ }
+ }
+ }
}
ps++;
c = *ps;
T128 = round_const_table_128[rmode][extra_digits];
__add_carry_out (CX1.w[0], carry, T128.w[0], CX.w[0]);
CX1.w[1] = CX.w[1] + T128.w[1] + carry;
- if (__unsigned_compare_ge_128
- (CX1, power10_table_128[extra_digits + 16]))
- uf_check = 0;
}
extra_digits =
extra_digits + DECIMAL_EXPONENT_BIAS_128 -
{ if ((x & (0xFull<<27)) == (0xFull<<27)) \
{ if ((x & (0x1Full<<26)) != (0x1Full<<26)) inf; \
if ((x & (1ul<<25))!=0) *pfpsf |= INVALID_EXCEPTION; \
- nan(s,((((x) & 0xFFFFul) > 999999ul) ? 0 : \
+ nan(s,((((x) & 0xFFFFFul) > 999999ul) ? 0 : \
(((unsigned long long) x) << 44)),0ull); \
} \
e = ((x >> 21) & ((1ull<<8)-1)) - 101; \
c = (1ull<<23) + (x & ((1ull<<21)-1)); \
- if ((unsigned long)(c) > 9999999ul) c = 0; \
+ if ((unsigned long)(c) > 9999999ul) zero; \
k = 0; \
} \
else \
{ e = ((x >> 23) & ((1ull<<8)-1)) - 101; \
c = x & ((1ull<<23)-1); \
if (c == 0) zero; \
- k = clz32(c) - 8; \
+ k = clz32_nz(c) - 8; \
c = c << k; \
} \
}
} \
e = ((x >> 51) & ((1ull<<10)-1)) - 398; \
c = (1ull<<53) + (x & ((1ull<<51)-1)); \
- if ((unsigned long long)(c) > 9999999999999999ull) c = 0; \
+ if ((unsigned long long)(c) > 9999999999999999ull) zero; \
k = 0; \
} \
else \
{ e = ((x >> 53) & ((1ull<<10)-1)) - 398; \
c = x & ((1ull<<53)-1); \
if (c == 0) zero; \
- k = clz64(c) - 10; \
+ k = clz64_nz(c) - 10; \
c = c << k; \
} \
}
// We actually check if e >= ceil((sci_emax + 1) * log_10(2))
// which in this case is e >= ceil(1024 * log_10(2)) = ceil(308.25) = 309
- if (e >= 309) {
- *pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
- return_binary64_ovf (s);
- }
-// 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
m_min = (breakpoints_binary64 + 358)[e];
// Choose provisional exponent and reciprocal multiplier based on breakpoint
- if (le128 (c.w[1], c.w[0], m_min.w[1], m_min.w[0])) {
+ if (c.w[1] < m_min.w[1]) {
r = (multipliers1_binary64 + 358)[e];
} else {
r = (multipliers2_binary64 + 358)[e];
// Do the reciprocal multiplication
- __mul_128x256_to_384 (z, c, r)
+ __mul_64x256_to_320(z, c.w[1], r);
+ z.w[5]=z.w[4]; z.w[4]=z.w[3]; z.w[3]=z.w[2]; z.w[2]=z.w[1]; z.w[1]=z.w[0]; z.w[0]=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 (z.w[5], z.w[4], z.w[3], z.w[2], d);
- }
+
c_prov = z.w[5];
// Round using round-sticky words
w[1],
roundbound_128[(rnd_mode << 2) + ((s & 1) << 1) +
(c_prov & 1)].w[0], z.w[4], z.w[3])) {
- c_prov = c_prov + 1;
- if (c_prov == (1ull << 53)) {
- c_prov = 1ull << 52;
- e_out = e_out + 1;
- }
- }
-// Check for overflow
-
- if (e_out >= 2047) {
- *pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
- return_binary64_ovf (s);
+ c_prov = c_prov + 1;
}
-// Modify exponent for a tiny result, otherwise lop the implicit bit
-
- if (c_prov < (1ull << 52))
- e_out = 0;
- else
- c_prov = c_prov & ((1ull << 52) - 1);
+ c_prov = c_prov & ((1ull << 52) - 1);
// Set the inexact and underflow flag as appropriate
if ((z.w[4] != 0) || (z.w[3] != 0)) {
*pfpsf |= INEXACT_EXCEPTION;
- if (e_out == 0)
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
// Package up the result as a binary floating-point number
__mul_128x256_to_384 (z, c, r)
c_prov = z.w[5];
+// Test inexactness and underflow (when testing tininess before rounding)
+
+ if ((z.w[4] != 0) || (z.w[3] != 0)) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ if (c_prov < 1000000ull)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
// Round using round-sticky words
// If we spill over into the next decade, correct
// Flag underflow where it may be needed even for |result| = SNN
if (c_prov == 10000000ull) {
c_prov = 1000000ull;
e_out = e_out + 1;
- } else if ((c_prov == 1000000ull) && (e_out == 0)) {
- if ((((rnd_mode & 3) == 0) && (z.w[4] <= 17524406870024074035ull))
- || ((rnd_mode + (s & 1) == 2)
- && (z.w[4] <= 16602069666338596454ull)))
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
}
+
// Check for overflow
if (e_out > 90 + 101) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_bid32_ovf (s);
}
-// Set the inexact flag as appropriate and check underflow
-// It's no doubt superfluous to check inexactness, but anyway...
- if ((z.w[4] != 0) || (z.w[3] != 0)) {
- *pfpsf |= INEXACT_EXCEPTION;
- if (c_prov < 1000000ull)
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
// Package up the result
return_bid32 (s, e_out, c_prov);
__mul_128x256_to_384 (z, c, r)
c_prov = z.w[5];
+// Test inexactness and underflow (when testing tininess before rounding)
+
+ if ((z.w[4] != 0) || (z.w[3] != 0)) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ if (c_prov < 1000000ull)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
// Round using round-sticky words
// If we spill over into the next decade, correct
// Flag underflow where it may be needed even for |result| = SNN
if (c_prov == 10000000ull) {
c_prov = 1000000ull;
e_out = e_out + 1;
- } else if ((c_prov == 1000000ull) && (e_out == 0)) {
- if ((((rnd_mode & 3) == 0) && (z.w[4] <= 17524406870024074035ull))
- || ((rnd_mode + (s & 1) == 2)
- && (z.w[4] <= 16602069666338596454ull)))
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
}
+
// Check for overflow
if (e_out > 90 + 101) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_bid32_ovf (s);
}
-// Set the inexact flag as appropriate and check underflow
-// It's no doubt superfluous to check inexactness, but anyway...
- if ((z.w[4] != 0) || (z.w[3] != 0)) {
- *pfpsf |= INEXACT_EXCEPTION;
- if (c_prov < 1000000ull)
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
// Package up the result
return_bid32 (s, e_out, c_prov);
__mul_128x256_to_384 (z, c, r)
c_prov = z.w[5];
+// Test inexactness and underflow (when testing tininess before rounding)
+ if ((z.w[4] != 0) || (z.w[3] != 0)) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ if (c_prov < 1000000ull)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
// Round using round-sticky words
// If we spill over into the next decade, correct
// Flag underflow where it may be needed even for |result| = SNN
if (c_prov == 10000000ull) {
c_prov = 1000000ull;
e_out = e_out + 1;
- } else if ((c_prov == 1000000ull) && (e_out == 0)) {
- if ((((rnd_mode & 3) == 0) &&
- le128 (z.w[4], z.w[3],
- 17524406870024074035ull, 3689348814741910323ull)) ||
- ((rnd_mode + (s & 1) == 2) &&
- le128 (z.w[4], z.w[3],
- 16602069666338596454ull, 7378697629483820646ull)))
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
}
+
// Check for overflow
if (e_out > 90 + 101) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_bid32_ovf (s);
}
-// Set the inexact flag as appropriate and check underflow
-// It's no doubt superfluous to check inexactness, but anyway...
- if ((z.w[4] != 0) || (z.w[3] != 0)) {
- *pfpsf |= INEXACT_EXCEPTION;
- if (c_prov < 1000000ull)
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
// Package up the result
return_bid32 (s, e_out, c_prov);
__mul_128x256_to_384 (z, c, r)
c_prov = z.w[5];
+// Test inexactness and underflow (when testing tininess before rounding)
+
+ if ((z.w[4] != 0) || (z.w[3] != 0)) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ if (c_prov < 1000000000000000ull)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
// Round using round-sticky words
// If we spill over into the next decade, correct
// Flag underflow where it may be needed even for |result| = SNN
if (c_prov == 10000000000000000ull) {
c_prov = 1000000000000000ull;
e_out = e_out + 1;
- } else if ((c_prov == 1000000000000000ull) && (e_out == 0)) {
- if ((((rnd_mode & 3) == 0) && (z.w[4] <= 17524406870024074035ull))
- || ((rnd_mode + (s & 1) == 2)
- && (z.w[4] <= 16602069666338596454ull)))
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
}
+
// Check for overflow
if (e_out > 369 + 398) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_bid64_ovf (s);
}
-// Set the inexact flag as appropriate and check underflow
-// It's no doubt superfluous to check inexactness, but anyway...
- if ((z.w[4] != 0) || (z.w[3] != 0)) {
- *pfpsf |= INEXACT_EXCEPTION;
- if (c_prov < 1000000000000000ull)
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
// Package up the result
return_bid64 (s, e_out, c_prov);
__mul_128x256_to_384 (z, c, r)
c_prov = z.w[5];
+// Test inexactness and underflow (when testing tininess before rounding)
+
+ if ((z.w[4] != 0) || (z.w[3] != 0)) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ if (c_prov < 1000000000000000ull)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
+ }
+
// Round using round-sticky words
// If we spill over into the next decade, correct
// Flag underflow where it may be needed even for |result| = SNN
if (c_prov == 10000000000000000ull) {
c_prov = 1000000000000000ull;
e_out = e_out + 1;
- } else if ((c_prov == 1000000000000000ull) && (e_out == 0)) {
- if ((((rnd_mode & 3) == 0) && (z.w[4] <= 17524406870024074035ull))
- || ((rnd_mode + (s & 1) == 2)
- && (z.w[4] <= 16602069666338596454ull)))
- *pfpsf |= UNDERFLOW_EXCEPTION;
}
}
+
// Check for overflow
if (e_out > 369 + 398) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_bid64_ovf (s);
}
-// Set the inexact flag as appropriate and check underflow
-// It's no doubt superfluous to check inexactness, but anyway...
- if ((z.w[4] != 0) || (z.w[3] != 0)) {
- *pfpsf |= INEXACT_EXCEPTION;
- if (c_prov < 1000000000000000ull)
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
// Package up the result
return_bid64 (s, e_out, c_prov);
#define BID_BIG_ENDIAN __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
#endif
+#if BID_BIG_ENDIAN
+#define BID_HIGH_128W 0
+#define BID_LOW_128W 1
+#else
+#define BID_HIGH_128W 1
+#define BID_LOW_128W 0
+#endif
+
#ifndef BID_THREAD
#if defined (HAVE_CC_TLS) && defined (USE_TLS)
#define BID_THREAD __thread
#endif
+#if defined __NO_BINARY80__
+#define __ENABLE_BINARY80__ 0
+#else
#if !defined _MSC_VER || defined __INTEL_COMPILER
#define __ENABLE_BINARY80__ 1
#endif
+#endif
#ifndef HPUX_OS
#define BINARY80 long double
} UINT256;
typedef unsigned int FPSC; // floating-point status and control
+ typedef enum class_types {
+ signalingNaN,
+ quietNaN,
+ negativeInfinity,
+ negativeNormal,
+ negativeSubnormal,
+ negativeZero,
+ positiveZero,
+ positiveSubnormal,
+ positiveNormal,
+ positiveInfinity
+ } class_t;
+
// TYPE parameters
#define BID128_MAXDIGITS 34
#define BID64_MAXDIGITS 16
extern UINT64 bid64_copySign (UINT64 x,
UINT64 y _EXC_MASKS_PARAM
_EXC_INFO_PARAM);
- extern int bid64_class (UINT64 x _EXC_MASKS_PARAM _EXC_INFO_PARAM);
+ extern class_t bid64_class (UINT64 x _EXC_MASKS_PARAM _EXC_INFO_PARAM);
extern int bid64_sameQuantum (UINT64 x, UINT64 y
_EXC_MASKS_PARAM _EXC_INFO_PARAM);
extern int bid64_totalOrder (UINT64 x, UINT64 y
extern UINT128 bid128_copySign (UINT128 x,
UINT128 y _EXC_MASKS_PARAM
_EXC_INFO_PARAM);
- extern int bid128_class (UINT128 x _EXC_MASKS_PARAM
- _EXC_INFO_PARAM);
+ extern class_t bid128_class (UINT128 x _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM);
extern int bid128_sameQuantum (UINT128 x,
UINT128 y _EXC_MASKS_PARAM
_EXC_INFO_PARAM);
coefficient_x += D;
}
break;
+ default: break; // default added to avoid compiler warning
}
if (coefficient_x < 1000000000000000ull) {
coefficient_x -= D;
} else if (FS.w[1] | FS.w[0])
CYh++;
break;
+ default: break; // default added to avoid compiler warning
}
#endif
#endif
// round up
if (sgn)
r = SMALLEST_BID64;
+ default:
+ break;
}
return r;
}
// round up
if (sgn)
r = SMALLEST_BID64;
+ default:
+ break;
}
return r;
}
A=((tempx.i >>23) & EXPONENT_MASK32) - 0x7f;\
}
- enum class_types {
- signalingNaN,
- quietNaN,
- negativeInfinity,
- negativeNormal,
- negativeSubnormal,
- negativeZero,
- positiveZero,
- positiveSubnormal,
- positiveNormal,
- positiveInfinity
- };
-
typedef union {
UINT64 ui64;
double d;