Update libbid according to the latest Intel Decimal Floating-Point Math Library.

The Intel Decimal Floating-Point Math Library is available as open-source on Netlib[1].

[1] https://www.netlib.org/misc/intel/.

libgcc/config/libbid/ChangeLog:

	* bid128_fma.c (add_and_round): Fix bug: the result
	of (+5E+368)*(+10E-34)+(-10E+369) was returning
	-9999999999999999999999999999999999E+336 instead of expected
	result -1000000000000000000000000000000000E+337.
	(bid128_ext_fma): Ditto.
	(bid64qqq_fma): Ditto.
	* bid128_noncomp.c: Change return type of bid128_class from
	int to class_t.
	* bid128_round_integral.c: Add default case to avoid compiler
	warning.
	* bid128_string.c (bid128_to_string): Replace 0x30 with '0'
	for zero digit.
	(bid128_from_string): Ditto.
	* bid32_to_bid128.c (bid128_to_bid32): Fix Bug. In addition
	to the INEXACT flag, the UNDERFLOW flag needs to be set (and
	was not) when converting an input such as
	+6931674235302037148946035460357709E+1857 to +1000000E-101
	* bid32_to_bid64.c (bid64_to_bid32): fix Bug, In addition to
	the INEXACT flag, the UNDERFLOW flag needs to be set (and was
	not) when converting an input such as +9999999000000001E-111
	to +1000000E-101. Furthermore, significant bits of NaNs are
	set correctly now. For example,  0x7c00003b9aca0000 was
	returning 0x7c000002 instead of 0x 7c000100.
	* bid64_noncomp.c: Change return type of bid64_class from int
	to class_t.
	* bid64_round_integral.c (bid64_round_integral_exact): Add
	default case to avoid compiler warning.
	* bid64_string.c (bid64_from_string): Fix bug for rounding
	up. The input string "10000000000000000" was returning
	+1000000000000001E+1 instead of +1000000000000000E+1.
	* bid64_to_bid128.c (bid128_to_bid64): Fix bug, in addition to
	the INEXACT flag, the UNDERFLOW flag needs to be set (and was
	not) when converting an input such as
	+9999999999999999999999999999999999E-417 to
	+1000000000000000E-398.
	* bid_binarydecimal.c (bid32_to_binary64): Fix bug for
	conversion between binary and bid types. For example,
	0x7c0F4240 was returning 0x7FFFA12000000000 instead of
	expected double precision 0x7FF8000000000000.
	(binary64_to_bid32): Ditto.
	(binary80_to_bid32): Ditto.
	(binary128_to_bid32): Ditto.
	(binary80_to_bid64): Ditto.
	(binary128_to_bid64): Ditto.
	* bid_conf.h (BID_HIGH_128W): New macro.
	(BID_LOW_128W): Ditto.
	* bid_functions.h (__ENABLE_BINARY80__): Ditto.
	(ALIGN): Ditto.
	* bid_inline_add.h (get_add128): Add default case to avoid compiler
	warning.
	* bid_internal.h (get_BID64): Ditto.
	(fast_get_BID64_check_OF): Ditto.
	(ALIGN): New macro.

	Co-authored-by: Anderson, Cristina S <cristina.s.anderson@intel.com>
	Co-authored-by: Akkas, Ahmet <ahmet.akkas@intel.com>
	Co-authored-by: Cornea, Marius <marius.cornea@intel.com>
This commit is contained in:
liuhongt 2024-03-27 08:20:13 +08:00
parent ffb521f157
commit affd77d3fe
15 changed files with 220 additions and 238 deletions

View file

@ -417,13 +417,12 @@ add_and_round (int q3,
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
@ -434,10 +433,6 @@ add_and_round (int q3,
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
@ -851,7 +846,6 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
}
}
}
p_sign = x_sign ^ y_sign; // sign of the product
// identify cases where at least one operand is infinity
@ -988,15 +982,10 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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 =
@ -1011,42 +1000,36 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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++;
}
}
@ -1055,32 +1038,25 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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++;
}
}
@ -1128,7 +1104,6 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
} 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
@ -1232,22 +1207,18 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
// 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
@ -1283,7 +1254,6 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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;
@ -1319,10 +1289,11 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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)
@ -1648,7 +1619,6 @@ bid128_ext_fma (int *ptr_is_midpoint_lt_even,
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')
@ -1736,7 +1706,7 @@ delta_ge_zero:
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
@ -1795,7 +1765,10 @@ delta_ge_zero:
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;
@ -1826,11 +1799,7 @@ delta_ge_zero:
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
@ -1877,10 +1846,9 @@ delta_ge_zero:
// 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) {
@ -2594,7 +2562,7 @@ delta_ge_zero:
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) {
@ -2678,6 +2646,9 @@ delta_ge_zero:
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;
@ -3328,9 +3299,6 @@ delta_ge_zero:
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
@ -3521,6 +3489,19 @@ delta_ge_zero:
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
@ -4162,21 +4143,34 @@ bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z
// 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
@ -4453,6 +4447,20 @@ bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z
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);
}

View file

@ -443,7 +443,7 @@ void
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;

View file

@ -177,6 +177,7 @@ case ROUNDING_TO_ZERO:
BID_RETURN (res);
}
break;
default: break; // default added to avoid compiler warning
}
// q = nr. of decimal digits in x
@ -804,6 +805,7 @@ case ROUNDING_TO_ZERO:
BID_RETURN (res);
}
break;
default: break; // default added to avoid compiler warning
}
BID_RETURN (res);

View file

@ -56,6 +56,7 @@ bid128_to_string (char *str, UINT128 x
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;
@ -232,14 +233,14 @@ bid128_to_string (char *str, UINT128 x
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];
@ -643,7 +644,7 @@ bid128_from_string (char *ps _RND_MODE_PARAM _EXC_FLAGS_PARAM
}
break;
default: break; // default added to avoid compiler warning
}
// now form the coefficient as coeff_high*10^17+coeff_low+carry
scale_high = 100000000000000000ull;

View file

@ -155,9 +155,6 @@ bid128_to_bid32 (UINT128 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
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 -

View file

@ -79,6 +79,7 @@ bid64_to_bid32 (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
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;
@ -93,8 +94,10 @@ bid64_to_bid32 (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
// 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
@ -139,10 +142,6 @@ bid64_to_bid32 (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
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;
}

View file

@ -358,7 +358,7 @@ void
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;

View file

@ -142,6 +142,7 @@ bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
BID_RETURN (res);
}
break;
default: break; // default added to avoid compiler warning
} // end switch ()
// q = nr. of decimal digits in x (1 <= q <= 54)
@ -483,6 +484,7 @@ bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
BID_RETURN (res);
}
break;
default: break; // default added to avoid compiler warning
} // end switch ()
BID_RETURN (res);
}

View file

@ -251,7 +251,7 @@ bid64_from_string (char *ps
#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;
@ -419,10 +419,10 @@ bid64_from_string (char *ps
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; }
@ -443,8 +443,21 @@ bid64_from_string (char *ps
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;

View file

@ -153,9 +153,6 @@ bid128_to_bid64 (UINT128 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
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 -

View file

@ -566,19 +566,19 @@ BID_BINARY80LDOUBLE;
{ 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; \
} \
}
@ -594,14 +594,14 @@ BID_BINARY80LDOUBLE;
} \
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; \
} \
}
@ -144302,20 +144302,6 @@ bid32_to_binary64 (UINT32 x
// 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];
@ -144323,7 +144309,7 @@ bid32_to_binary64 (UINT32 x
// 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];
@ -144332,17 +144318,12 @@ bid32_to_binary64 (UINT32 x
// 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
@ -144353,31 +144334,14 @@ bid32_to_binary64 (UINT32 x
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;
}
c_prov = c_prov + 1;
}
// Check for overflow
if (e_out >= 2047) {
*pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
return_binary64_ovf (s);
}
// 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
@ -145756,6 +145720,14 @@ binary64_to_bid32 (double x
__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
@ -145769,27 +145741,16 @@ binary64_to_bid32 (double x
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);
@ -145919,6 +145880,14 @@ binary80_to_bid32 (BINARY80 x
__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
@ -145932,27 +145901,16 @@ binary80_to_bid32 (BINARY80 x
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);
@ -146071,6 +146029,13 @@ binary128_to_bid32 (BINARY128 x
__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
@ -146086,30 +146051,16 @@ binary128_to_bid32 (BINARY128 x
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);
@ -146562,6 +146513,14 @@ binary80_to_bid64 (BINARY80 x
__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
@ -146575,27 +146534,16 @@ binary80_to_bid64 (BINARY80 x
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);
@ -146723,6 +146671,14 @@ binary128_to_bid64 (BINARY128 x
__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
@ -146736,27 +146692,16 @@ binary128_to_bid64 (BINARY128 x
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);

View file

@ -519,6 +519,14 @@ see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
#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

View file

@ -67,9 +67,13 @@ ALIGN (16)
#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
@ -91,6 +95,19 @@ ALIGN (16)
} 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
@ -2948,7 +2965,7 @@ ALIGN (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
@ -2984,8 +3001,8 @@ ALIGN (16)
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);

View file

@ -918,6 +918,7 @@ get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
coefficient_x += D;
}
break;
default: break; // default added to avoid compiler warning
}
if (coefficient_x < 1000000000000000ull) {
coefficient_x -= D;
@ -1107,6 +1108,7 @@ get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
} else if (FS.w[1] | FS.w[0])
CYh++;
break;
default: break; // default added to avoid compiler warning
}
#endif
#endif

View file

@ -970,6 +970,8 @@ get_BID64 (UINT64 sgn, int expon, UINT64 coeff, int rmode,
// round up
if (sgn)
r = SMALLEST_BID64;
default:
break;
}
return r;
}
@ -1086,6 +1088,8 @@ fast_get_BID64_check_OF (UINT64 sgn, int expon, UINT64 coeff, int rmode,
// round up
if (sgn)
r = SMALLEST_BID64;
default:
break;
}
return r;
}
@ -2582,19 +2586,6 @@ ALIGN (16)
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;