*
****************************************************************************/
-#include "bid_conf.h"
-#include "bid_functions.h"
#include "bid_internal.h"
static void
-rounding_correction (
- unsigned int rnd_mode,
- unsigned int is_inexact_lt_midpoint,
- unsigned int is_inexact_gt_midpoint,
- unsigned int is_midpoint_lt_even,
- unsigned int is_midpoint_gt_even,
- int unbexp,
- UINT128 * ptrres,
- _IDEC_flags * ptrfpsf) {
- // unbiased true exponent unbexp may be larger than emax
+rounding_correction (unsigned int rnd_mode,
+ unsigned int is_inexact_lt_midpoint,
+ unsigned int is_inexact_gt_midpoint,
+ unsigned int is_midpoint_lt_even,
+ unsigned int is_midpoint_gt_even,
+ int unbexp,
+ UINT128 * ptrres, _IDEC_flags * ptrfpsf) {
+ // unbiased true exponent unbexp may be larger than emax
UINT128 res = *ptrres; // expected to have the correct sign and coefficient
- // (the exponent field is ignored, as unbexp is used instead)
+ // (the exponent field is ignored, as unbexp is used instead)
UINT64 sign, exp;
UINT64 C_hi, C_lo;
exp = (UINT64) (unbexp + 6176) << 49; // valid only if expmin<=unbexp<=expmax
C_hi = res.w[1] & MASK_COEFF;
C_lo = res.w[0];
- if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
+ if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) &&
- is_midpoint_gt_even))) ||
- (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
+ is_midpoint_gt_even))) ||
+ (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) &&
is_midpoint_gt_even)))) {
// C = C + 1
exp = (UINT64) (unbexp + 6176) << 49;
}
} else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) &&
- ((sign && (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) ||
+ ((sign && (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) ||
(!sign && (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TO_ZERO)))) {
// C = C - 1
C_lo = C_lo - 1;
if (C_lo == 0xffffffffffffffffull)
C_hi--;
// check if we crossed into the lower decade
- if (exp > 0 && C_hi == 0x0000314dc6448d93ull &&
- C_lo == 0x38c15b09ffffffffull) { // 10^33 - 1
- C_hi = 0x0001ed09bead87c0ull; // 10^34 - 1
- C_lo = 0x378d8e63ffffffffull;
- // exp = exp - EXP_P1;
- unbexp = unbexp - 1;
- exp = (UINT64) (unbexp + 6176) << 49;
+ if (C_hi == 0x0000314dc6448d93ull && C_lo == 0x38c15b09ffffffffull) {
+ // C = 10^33 - 1
+ if (exp > 0) {
+ C_hi = 0x0001ed09bead87c0ull; // 10^34 - 1
+ C_lo = 0x378d8e63ffffffffull;
+ // exp = exp - EXP_P1;
+ unbexp = unbexp - 1;
+ exp = (UINT64) (unbexp + 6176) << 49;
+ } else { // if exp = 0
+ if (is_midpoint_lt_even || is_midpoint_lt_even ||
+ is_inexact_gt_midpoint || is_inexact_gt_midpoint) // tiny & inexact
+ *ptrfpsf |= UNDERFLOW_EXCEPTION;
+ }
}
} else {
; // the result is already correct
static int
-__bid_nr_digits256 (UINT256 R256) {
+nr_digits256 (UINT256 R256) {
int ind;
// determine the number of decimal digits in R256
if (R256.w[3] == 0x0 && R256.w[2] == 0x0 && R256.w[1] == 0x0) {
// between 1 and 19 digits
for (ind = 1; ind <= 19; ind++) {
- if (R256.w[0] < __bid_ten2k64[ind]) {
+ if (R256.w[0] < ten2k64[ind]) {
break;
}
}
// ind digits
} else if (R256.w[3] == 0x0 && R256.w[2] == 0x0 &&
- (R256.w[1] < __bid_ten2k128[0].w[1]
- || (R256.w[1] == __bid_ten2k128[0].w[1]
- && R256.w[0] < __bid_ten2k128[0].w[0]))) {
+ (R256.w[1] < ten2k128[0].w[1] ||
+ (R256.w[1] == ten2k128[0].w[1]
+ && R256.w[0] < ten2k128[0].w[0]))) {
// 20 digits
ind = 20;
} else if (R256.w[3] == 0x0 && R256.w[2] == 0x0) {
// between 21 and 38 digits
for (ind = 1; ind <= 18; ind++) {
- if (R256.w[1] < __bid_ten2k128[ind].w[1] ||
- (R256.w[1] == __bid_ten2k128[ind].w[1] &&
- R256.w[0] < __bid_ten2k128[ind].w[0])) {
+ if (R256.w[1] < ten2k128[ind].w[1] ||
+ (R256.w[1] == ten2k128[ind].w[1] &&
+ R256.w[0] < ten2k128[ind].w[0])) {
break;
}
}
// ind + 20 digits
ind = ind + 20;
- } else if (((R256.w[3] == 0x0 &&
- (R256.w[2] < __bid_ten2k256[0].w[2])) ||
- (R256.w[2] == __bid_ten2k256[0].w[2] &&
- R256.w[1] < __bid_ten2k256[0].w[1]) ||
- (R256.w[2] == __bid_ten2k256[0].w[2] &&
- R256.w[1] == __bid_ten2k256[0].w[1] &&
- R256.w[0] < __bid_ten2k256[0].w[0]))) {
+ } else if (R256.w[3] == 0x0 &&
+ (R256.w[2] < ten2k256[0].w[2] ||
+ (R256.w[2] == ten2k256[0].w[2] &&
+ R256.w[1] < ten2k256[0].w[1]) ||
+ (R256.w[2] == ten2k256[0].w[2] &&
+ R256.w[1] == ten2k256[0].w[1] &&
+ R256.w[0] < ten2k256[0].w[0]))) {
// 39 digits
ind = 39;
} else {
// between 40 and 68 digits
for (ind = 1; ind <= 29; ind++) {
- if (R256.w[3] < __bid_ten2k256[ind].w[3] ||
- (R256.w[3] == __bid_ten2k256[ind].w[3] &&
- R256.w[2] < __bid_ten2k256[ind].w[2]) ||
- (R256.w[3] == __bid_ten2k256[ind].w[3] &&
- R256.w[2] == __bid_ten2k256[ind].w[2] &&
- R256.w[1] < __bid_ten2k256[ind].w[1]) ||
- (R256.w[3] == __bid_ten2k256[ind].w[3] &&
- R256.w[2] == __bid_ten2k256[ind].w[2] &&
- R256.w[1] == __bid_ten2k256[ind].w[1] &&
- R256.w[0] < __bid_ten2k256[ind].w[0])) {
+ if (R256.w[3] < ten2k256[ind].w[3] ||
+ (R256.w[3] == ten2k256[ind].w[3] &&
+ R256.w[2] < ten2k256[ind].w[2]) ||
+ (R256.w[3] == ten2k256[ind].w[3] &&
+ R256.w[2] == ten2k256[ind].w[2] &&
+ R256.w[1] < ten2k256[ind].w[1]) ||
+ (R256.w[3] == ten2k256[ind].w[3] &&
+ R256.w[2] == ten2k256[ind].w[2] &&
+ R256.w[1] == ten2k256[ind].w[1] &&
+ R256.w[0] < ten2k256[ind].w[0])) {
break;
}
}
R256.w[0] = C3.w[0];
} else if (scale <= 19) { // 10^scale fits in 64 bits
P128.w[1] = 0;
- P128.w[0] = __bid_ten2k64[scale];
+ P128.w[0] = ten2k64[scale];
__mul_128x128_to_256 (R256, P128, C3);
} else if (scale <= 38) { // 10^scale fits in 128 bits
- __mul_128x128_to_256 (R256, __bid_ten2k128[scale - 20], C3);
+ __mul_128x128_to_256 (R256, ten2k128[scale - 20], C3);
} else if (scale <= 57) { // 39 <= scale <= 57
// 10^scale fits in 192 bits but C3 * 10^scale fits in 223 or 230 bits
// (10^67 has 223 bits; 10^69 has 230 bits);
// 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127
// bits and so 10^(scale-38) * C3 fits in 128 bits with certainty
// Note that 1 <= scale - 38 <= 19 => 10^(scale-38) fits in 64 bits
- __mul_64x128_to_128 (R128, __bid_ten2k64[scale - 38], C3);
+ __mul_64x128_to_128 (R128, ten2k64[scale - 38], C3);
// now multiply R128 by 10^38
- __mul_128x128_to_256 (R256, R128, __bid_ten2k128[18]);
+ __mul_128x128_to_256 (R256, R128, ten2k128[18]);
} else { // 58 <= scale <= 66
// 10^scale takes between 193 and 220 bits,
// and C3 * 10^scale fits in 223 bits (10^67/10^69 has 223/230 bits)
// Note that 20 <= scale - 38 <= 30 => 10^(scale-38) fits in 128 bits
// Calculate first 10^(scale-38) * C3, which fits in 128 bits; because
// 10^(scale-38) takes more than 64 bits, C3 will take less than 64
- __mul_64x128_to_128 (R128, C3.w[0], __bid_ten2k128[scale - 58]);
+ __mul_64x128_to_128 (R128, C3.w[0], ten2k128[scale - 58]);
// now calculate 10*38 * 10^(scale-38) * C3
- __mul_128x128_to_256 (R256, R128, __bid_ten2k128[18]);
+ __mul_128x128_to_256 (R256, R128, ten2k128[18]);
}
// C3 * 10^scale is now in R256
// but may require rounding
// compare first R256 = C3 * 10^scale and C4
- if (R256.w[3] > C4.w[3] || (R256.w[3] == C4.w[3] && R256.w[2] > C4.w[2]) ||
+ if (R256.w[3] > C4.w[3] || (R256.w[3] == C4.w[3] && R256.w[2] > C4.w[2]) ||
(R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] > C4.w[1]) ||
(R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] == C4.w[1] &&
R256.w[0] >= C4.w[0])) { // C3 * 10^scale >= C4
}
// determine the number of decimal digits in R256
- ind = __bid_nr_digits256 (R256);
+ ind = nr_digits256 (R256);
// the exact result is (-1)^p_sign * R256 * 10^e4 where q (R256) = ind;
// round to the destination precision, with unbounded exponent
if (ind <= 38) {
P128.w[1] = R256.w[1];
P128.w[0] = R256.w[0];
- __bid_round128_19_38 (ind, x0, P128, &R128, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round128_19_38 (ind, x0, P128, &R128, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
} else if (ind <= 57) {
P192.w[2] = R256.w[2];
P192.w[1] = R256.w[1];
P192.w[0] = R256.w[0];
- __bid_round192_39_57 (ind, x0, P192, &R192, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round192_39_57 (ind, x0, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
R128.w[1] = R192.w[1];
R128.w[0] = R192.w[0];
} else { // if (ind <= 68)
- __bid_round256_58_76 (ind, x0, R256, &R256, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round256_58_76 (ind, x0, R256, &R256, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
R128.w[1] = R256.w[1];
R128.w[0] = R256.w[0];
}
P128.w[1] = p_sign | 0x3040000000000000ull | R128.w[1];
P128.w[0] = R128.w[0];
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint, is_midpoint_lt_even,
- is_midpoint_gt_even, 0, &P128, ptrfpsf);
+ is_inexact_lt_midpoint,
+ 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) {
R128.w[1] = res.w[1] & MASK_COEFF;
R128.w[0] = res.w[0];
if (ind <= 19) {
- if (R128.w[0] < __bid_midpoint64[ind - 1]) { // < 1/2 ulp
+ if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (R128.w[0] == __bid_midpoint64[ind - 1]) { // = 1/2 ulp
+ } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
is_inexact_gt_midpoint = 1;
}
} else { // if (ind <= 38) {
- if (R128.w[1] < __bid_midpoint128[ind - 20].w[1] ||
- (R128.w[1] == __bid_midpoint128[ind - 20].w[1] &&
- R128.w[0] < __bid_midpoint128[ind - 20].w[0])) { // < 1/2 ulp
+ if (R128.w[1] < midpoint128[ind - 20].w[1] ||
+ (R128.w[1] == midpoint128[ind - 20].w[1] &&
+ R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (R128.w[1] == __bid_midpoint128[ind - 20].w[1] &&
- R128.w[0] == __bid_midpoint128[ind - 20].w[0]) { // = 1/2 ulp
+ } else if (R128.w[1] == midpoint128[ind - 20].w[1] &&
+ R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
// round the ind-digit result to ind - x0 digits
if (ind <= 18) { // 2 <= ind <= 18
- __bid_round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
res.w[1] = 0x0;
res.w[0] = R64;
} else if (ind <= 38) {
P128.w[1] = res.w[1] & MASK_COEFF;
P128.w[0] = res.w[0];
- __bid_round128_19_38 (ind, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (ind, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
}
e4 = e4 + x0; // expmin
// we want the exponent to be expmin, so if incr_exp = 1 then
// 64 x 128 -> 128
P128.w[1] = res.w[1] & MASK_COEFF;
P128.w[0] = res.w[0];
- __mul_64x128_to_128 (res, __bid_ten2k64[1], P128);
+ __mul_64x128_to_128 (res, ten2k64[1], P128);
}
res.w[1] =
p_sign | ((UINT64) (e4 + 6176) << 49) | (res.w[1] & MASK_COEFF);
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
// apply correction if not rounding to nearest
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint, is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, ptrfpsf);
+ is_inexact_lt_midpoint, is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, ptrfpsf);
}
if (is_midpoint_lt_even || is_midpoint_gt_even ||
is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
#if DECIMAL_CALL_BY_REFERENCE
-void
-__bid128_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
- _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
+static void
+bid128_ext_fma (int *ptr_is_midpoint_lt_even,
+ int *ptr_is_midpoint_gt_even,
+ int *ptr_is_inexact_lt_midpoint,
+ int *ptr_is_inexact_gt_midpoint, UINT128 * pres,
+ UINT128 * px, UINT128 * py,
+ UINT128 *
+ pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
UINT128 x = *px, y = *py, z = *pz;
#if !DECIMAL_GLOBAL_ROUNDING
unsigned int rnd_mode = *prnd_mode;
#endif
#else
-UINT128
-__bid128_fma (UINT128 x, UINT128 y, UINT128 z
- _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
+static UINT128
+bid128_ext_fma (int *ptr_is_midpoint_lt_even,
+ int *ptr_is_midpoint_gt_even,
+ int *ptr_is_inexact_lt_midpoint,
+ int *ptr_is_inexact_gt_midpoint, UINT128 x, UINT128 y,
+ UINT128 z _RND_MODE_PARAM _EXC_FLAGS_PARAM
+ _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
- UINT128 res = {{ 0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull }};
- UINT64 x_sign, y_sign, z_sign, p_sign;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT64 x_sign, y_sign, z_sign, p_sign, tmp_sign;
UINT64 x_exp = 0, y_exp = 0, z_exp = 0, p_exp;
int true_p_exp;
UINT128 C1, C2, C3;
UINT256 R256;
// the following are based on the table of special cases for fma; the NaN
- // behavior is similar to that of the IPF fma
+ // behavior is similar to that of the IA-64 Architecture fma
// identify cases where at least one operand is NaN
+ BID_SWAP128 (x);
+ BID_SWAP128 (y);
+ BID_SWAP128 (z);
if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
// if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
// check first for non-canonical NaN payload
*pfpsf |= INVALID_EXCEPTION;
}
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((z.w[1] & MASK_NAN) == MASK_NAN) { // z is NAN
// if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
*pfpsf |= INVALID_EXCEPTION;
}
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
// if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
res.w[0] = x.w[0];
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
// x, y, z are 0, f, or inf but not NaN => unpack the arguments and check
} else { // G0_G1 != 11
x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
+ (C1.w[1] == 0x0001ed09bead87c0ull &&
+ C1.w[0] > 0x378d8e63ffffffffull)) {
// x is non-canonical if coefficient is larger than 10^34 -1
C1.w[1] = 0;
C1.w[0] = 0;
} else { // G0_G1 != 11
y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bits
if (C2.w[1] > 0x0001ed09bead87c0ull ||
- (C2.w[1] == 0x0001ed09bead87c0ull
- && C2.w[0] > 0x378d8e63ffffffffull)) {
+ (C2.w[1] == 0x0001ed09bead87c0ull &&
+ C2.w[0] > 0x378d8e63ffffffffull)) {
// y is non-canonical if coefficient is larger than 10^34 -1
C2.w[1] = 0;
C2.w[0] = 0;
} else { // G0_G1 != 11
z_exp = z.w[1] & MASK_EXP; // biased and shifted left 49 bits
if (C3.w[1] > 0x0001ed09bead87c0ull ||
- (C3.w[1] == 0x0001ed09bead87c0ull
- && C3.w[0] > 0x378d8e63ffffffffull)) {
+ (C3.w[1] == 0x0001ed09bead87c0ull &&
+ C3.w[0] > 0x378d8e63ffffffffull)) {
// z is non-canonical if coefficient is larger than 10^34 -1
C3.w[1] = 0;
C3.w[0] = 0;
// set invalid flag
*pfpsf |= INVALID_EXCEPTION;
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
// set invalid flag
*pfpsf |= INVALID_EXCEPTION;
} else {
- res.w[1] = z.w[1];
- res.w[0] = z.w[0];
+ res.w[1] = z_sign | MASK_INF;
+ res.w[0] = 0x0;
}
} else if (C1.w[1] == 0x0 && C1.w[0] == 0x0) { // x = 0
// z = 0, f, inf
res.w[1] = p_sign | MASK_INF;
res.w[0] = 0x0;
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
// x = 0, f and y = 0, f, necessarily
res.w[1] = z_sign | MASK_INF;
res.w[0] = 0x0;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
if (true_p_exp < -6176)
p_exp = 0; // cannot be less than EXP_MIN
- // MACH DEBUG else if (true_p_exp > 6111)
- // p_exp = (UINT64) (6111 + 6176) << 49; // cannot be more than EXP_MAX
else
p_exp = (UINT64) (true_p_exp + 6176) << 49;
- if (((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
- (C2.w[1] == 0x0 && C2.w[0] == 0x0)) &&
- C3.w[1] == 0x0 && C3.w[0] == 0x0) { // (x = 0 or y = 0) and z = 0
+ if (((C1.w[1] == 0x0 && C1.w[0] == 0x0) || (C2.w[1] == 0x0 && C2.w[0] == 0x0)) && C3.w[1] == 0x0 && C3.w[0] == 0x0) { // (x = 0 or y = 0) and z = 0
// the result is 0
if (p_exp < z_exp)
res.w[1] = p_exp; // preferred exponent
res.w[0] = 0x0;
}
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
// from this point on, we may need to know the number of decimal digits
x_nr_bits =
65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
- q1 = __bid_nr_digits[x_nr_bits - 1].digits;
+ q1 = nr_digits[x_nr_bits - 1].digits;
if (q1 == 0) {
- q1 = __bid_nr_digits[x_nr_bits - 1].digits1;
- if (C1.w[1] > __bid_nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == __bid_nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= __bid_nr_digits[x_nr_bits - 1].threshold_lo))
+ 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))
q1++;
}
}
64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
- q2 = __bid_nr_digits[y_nr_bits].digits;
+ q2 = nr_digits[y_nr_bits].digits;
if (q2 == 0) {
- q2 = __bid_nr_digits[y_nr_bits].digits1;
- if (C2.w[1] > __bid_nr_digits[y_nr_bits].threshold_hi ||
- (C2.w[1] == __bid_nr_digits[y_nr_bits].threshold_hi &&
- C2.w[0] >= __bid_nr_digits[y_nr_bits].threshold_lo))
+ 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++;
}
}
64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
}
- q3 = __bid_nr_digits[z_nr_bits].digits;
+ q3 = nr_digits[z_nr_bits].digits;
if (q3 == 0) {
- q3 = __bid_nr_digits[z_nr_bits].digits1;
- if (C3.w[1] > __bid_nr_digits[z_nr_bits].threshold_hi ||
- (C3.w[1] == __bid_nr_digits[z_nr_bits].threshold_hi &&
- C3.w[0] >= __bid_nr_digits[z_nr_bits].threshold_lo))
+ 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++;
}
}
- if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
+ if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
(C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
- // x = 0 or y = 0
+ // x = 0 or y = 0
// z = f, necessarily; for 0 + z return z, with the preferred exponent
// the result is z, but need to get the preferred exponent
if (z_exp <= p_exp) { // the preferred exponent is z_exp
res.w[0] = z.w[0];
} else if (q3 <= 19) { // z fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
}
} else { // z fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
}
// subtract scale from the exponent
z_exp = z_exp - ((UINT64) scale << 49);
res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else {
; // continue with x = f, y = f, z = 0 or x = f, y = f, z = f
if (q1 + q2 <= 19) { // if 2 <= q1 + q2 <= 19, C4 = C1 * C2 fits in 64 bits
C4.w[0] = C1.w[0] * C2.w[0];
// if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
- if (C4.w[0] < __bid_ten2k64[q1 + q2 - 1])
+ if (C4.w[0] < ten2k64[q1 + q2 - 1])
q4 = q1 + q2 - 1; // q4 in [1, 18]
else
q4 = q1 + q2; // q4 in [2, 19]
// q1 <= 19 and q2 <= 19 so both C1 and C2 fit in 64 bits
__mul_64x64_to_128MACH (C4, C1.w[0], C2.w[0]);
// if C4 < 10^(q1+q2-1) = 10^19 then q4 = q1+q2-1 = 19 else q4 = q1+q2 = 20
- if (C4.w[1] == 0 && C4.w[0] < __bid_ten2k64[19]) { // 19 = q1+q2-1
+ if (C4.w[1] == 0 && C4.w[0] < ten2k64[19]) { // 19 = q1+q2-1
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
q4 = 19; // 19 = q1 + q2 - 1
} else {
__mul_128x64_to_128 (C4, C2.w[0], C1);
}
// if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
- if (C4.w[1] < __bid_ten2k128[q1 + q2 - 21].w[1] ||
- (C4.w[1] == __bid_ten2k128[q1 + q2 - 21].w[1] &&
- C4.w[0] < __bid_ten2k128[q1 + q2 - 21].w[0])) {
+ if (C4.w[1] < ten2k128[q1 + q2 - 21].w[1] ||
+ (C4.w[1] == ten2k128[q1 + q2 - 21].w[1] &&
+ C4.w[0] < ten2k128[q1 + q2 - 21].w[0])) {
// if (C4.w[1] == 0) // q4 = 20, necessarily
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
// else
// may replace this by 128x128_to192
__mul_128x128_to_256 (C4, C1, C2); // C4.w[3] is 0
// if C4 < 10^(q1+q2-1) = 10^38 then q4 = q1+q2-1 = 38 else q4 = q1+q2 = 39
- if (C4.w[2] == 0 && (C4.w[1] < __bid_ten2k128[18].w[1] ||
- (C4.w[1] == __bid_ten2k128[18].w[1] && C4.w[0] < __bid_ten2k128[18].w[0]))) {
+ if (C4.w[2] == 0 && (C4.w[1] < ten2k128[18].w[1] ||
+ (C4.w[1] == ten2k128[18].w[1]
+ && C4.w[0] < ten2k128[18].w[0]))) {
// 18 = 38 - 20 = q1+q2-1 - 20
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
q4 = 38; // 38 = q1 + q2 - 1
__mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
}
// if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
- if (C4.w[2] < __bid_ten2k256[q1 + q2 - 40].w[2] ||
- (C4.w[2] == __bid_ten2k256[q1 + q2 - 40].w[2] &&
- (C4.w[1] < __bid_ten2k256[q1 + q2 - 40].w[1] ||
- (C4.w[1] == __bid_ten2k256[q1 + q2 - 40].w[1] &&
- C4.w[0] < __bid_ten2k256[q1 + q2 - 40].w[0])))) {
+ if (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
+ (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
+ (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
+ (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
+ C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))) {
// if (C4.w[2] == 0) // q4 = 39, necessarily
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
// else
__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] < __bid_ten2k256[18].w[2] ||
- (C4.w[2] == __bid_ten2k256[18].w[2] && (C4.w[1] < __bid_ten2k256[18].w[1] ||
- (C4.w[1] == __bid_ten2k256[18].w[1] && C4.w[0] < __bid_ten2k256[18].w[0]))))) {
- // 18 = 57 - 39 = q1+q2-1 - 39
+ 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]))))) {
+ // 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
} else {
// may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1);
__mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
// if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
- if (C4.w[3] < __bid_ten2k256[q1 + q2 - 40].w[3] ||
- (C4.w[3] == __bid_ten2k256[q1 + q2 - 40].w[3] &&
- (C4.w[2] < __bid_ten2k256[q1 + q2 - 40].w[2] ||
- (C4.w[2] == __bid_ten2k256[q1 + q2 - 40].w[2] &&
- (C4.w[1] < __bid_ten2k256[q1 + q2 - 40].w[1] ||
- (C4.w[1] == __bid_ten2k256[q1 + q2 - 40].w[1] &&
- C4.w[0] < __bid_ten2k256[q1 + q2 - 40].w[0])))))) {
+ if (C4.w[3] < ten2k256[q1 + q2 - 40].w[3] ||
+ (C4.w[3] == ten2k256[q1 + q2 - 40].w[3] &&
+ (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
+ (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
+ (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
+ (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
+ C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))))) {
// if (C4.w[3] == 0) // q4 = 58, necessarily
// length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
// else
if (q4 <= 38) {
P128.w[1] = C4.w[1];
P128.w[0] = C4.w[0];
- __bid_round128_19_38 (q4, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q4, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
} else if (q4 <= 57) { // 35 <= q4 <= 57
P192.w[2] = C4.w[2];
P192.w[1] = C4.w[1];
P192.w[0] = C4.w[0];
- __bid_round192_39_57 (q4, x0, P192, &R192, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
res.w[0] = R192.w[0];
res.w[1] = R192.w[1];
} else { // if (q4 <= 68)
- __bid_round256_58_76 (q4, x0, C4, &R256, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
res.w[0] = R256.w[0];
res.w[1] = R256.w[1];
}
// res = (C4 * 10^scale) * 10^expmax
if (q4 <= 19) { // C4 fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C4.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C4.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C4.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C4.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C4.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C4.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C4.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C4.w[0], ten2k128[scale - 20]);
}
} else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * CC43
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C4);
+ // 64 x 128 ten2k64[scale] * CC43
+ __mul_128x64_to_128 (res, ten2k64[scale], C4);
}
e4 = e4 - scale; // expmax
q4 = q4 + scale;
} else {
res.w[1] = p_sign | res.w[1];
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, pfpsf);
}
*pfpsf |= save_fpsf;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
// check for underflow
// the number of decimal digits in res is q4
if (x0 < q4) { // 1 <= x0 <= q4-1 => round res to q4 - x0 digits
if (q4 <= 18) { // 2 <= q4 <= 18, 1 <= x0 <= 17
- __bid_round64_2_18 (q4, x0, res.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round64_2_18 (q4, x0, res.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
if (incr_exp) {
// R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
- R64 = __bid_ten2k64[q4 - x0];
+ R64 = ten2k64[q4 - x0];
}
// res.w[1] = 0; (from above)
res.w[0] = R64;
// 19 <= q4 <= 38
P128.w[1] = res.w[1];
P128.w[0] = res.w[0];
- __bid_round128_19_38 (q4, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q4, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
if (incr_exp) {
// increase coefficient by a factor of 10; this will be <= 10^33
// R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
- // res.w[1] = 0;
- res.w[0] = __bid_ten2k64[q4 - x0];
+ // res.w[1] = 0;
+ res.w[0] = ten2k64[q4 - x0];
} else { // 20 <= q4 - x0 <= 37
- res.w[0] = __bid_ten2k128[q4 - x0 - 20].w[0];
- res.w[1] = __bid_ten2k128[q4 - x0 - 20].w[1];
+ res.w[0] = ten2k128[q4 - x0 - 20].w[0];
+ res.w[1] = ten2k128[q4 - x0 - 20].w[1];
}
}
}
// the second rounding is for 0.d(0)d(1)...d(q4-1) * 10^emin
// determine relationship with 1/2 ulp
if (q4 <= 19) {
- if (res.w[0] < __bid_midpoint64[q4 - 1]) { // < 1/2 ulp
+ if (res.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (res.w[0] == __bid_midpoint64[q4 - 1]) { // = 1/2 ulp
+ } else if (res.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
is_inexact_gt_midpoint = 1;
}
} else { // if (q4 <= 34)
- if (res.w[1] < __bid_midpoint128[q4 - 20].w[1] ||
- (res.w[1] == __bid_midpoint128[q4 - 20].w[1] &&
- res.w[0] < __bid_midpoint128[q4 - 20].w[0])) { // < 1/2 ulp
+ if (res.w[1] < midpoint128[q4 - 20].w[1] ||
+ (res.w[1] == midpoint128[q4 - 20].w[1] &&
+ res.w[0] < midpoint128[q4 - 20].w[0])) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (res.w[1] == __bid_midpoint128[q4 - 20].w[1] &&
- res.w[0] == __bid_midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
+ } else if (res.w[1] == midpoint128[q4 - 20].w[1] &&
+ res.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
}
// avoid a double rounding error
if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
- is_midpoint_lt_even) { // double rounding error upward
+ is_midpoint_lt_even) { // double rounding error upward
// res = res - 1
res.w[0]--;
if (res.w[0] == 0xffffffffffffffffull)
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
; // res and e4 are unchanged
} else if (q4 <= 19) { // C4 fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 res.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, res.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 res.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, res.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 res.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, res.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 res.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, res.w[0], ten2k128[scale - 20]);
}
} else { // res fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], res);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], res);
}
// subtract scale from the exponent
e4 = e4 - scale;
}
}
+
// check for inexact result
if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
is_midpoint_lt_even || is_midpoint_gt_even) {
res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, pfpsf);
}
*pfpsf |= save_fpsf;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
+ // no overflow, and no underflow for rounding to nearest
+ res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
+
+ if (rnd_mode != ROUNDING_TO_NEAREST) {
+ rounding_correction (rnd_mode,
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, pfpsf);
+ // if e4 = expmin && significand < 10^33 => result is tiny (for RD, RZ)
+ if (e4 == expmin) {
+ if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
+ ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
+ res.w[0] < 0x38c15b0a00000000ull)) {
+ is_tiny = 1;
+ }
+ }
+ }
- // no overflow and no underflow
if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
is_midpoint_lt_even || is_midpoint_gt_even) {
// set the inexact flag
*pfpsf |= INEXACT_EXCEPTION;
- }
- res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
-
- if (rnd_mode != ROUNDING_TO_NEAREST) {
- rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, pfpsf);
+ if (is_tiny)
+ *pfpsf |= UNDERFLOW_EXCEPTION;
}
if ((*pfpsf & INEXACT_EXCEPTION) == 0) { // x * y is exact
; // leave res unchanged
} else if (q4 <= 19) { // x * y fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
}
res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
} else { // x * y fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
}
} // else leave the result as it is, because p_exp <= z_exp
}
*pfpsf |= save_fpsf;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} // else we have f * f + f
delta_ge_zero:
if (delta >= 0) {
- if (p34 <= delta - 1 || // Case (1')
+ 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')
if ((q3 + e3) > (p34 + expmax) && p34 <= delta - 1) {
} else {
if (q3 <= 19) { // C3 fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0],
- __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0],
+ ten2k128[scale - 20]);
}
} else { // C3 fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
}
// the coefficient in res has q3 + scale = p34 digits
}
e3 = e3 - scale;
res.w[1] = z_sign | res.w[1];
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e3, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e3, &res, pfpsf);
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
// res = z
res.w[0] = C3.w[0];
} else if (q3 <= 19) { // z fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
}
} else { // z fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
}
// the coefficient in res has q3 + scale digits
// subtract scale from the exponent
if ((p_sign != z_sign) && (delta == (q3 + scale + 1))) {
// there is a gap of exactly one digit between the scaled C3 and C4
// C3 * 10^ scale = 10^(q3+scale-1) <=> C3 = 10^(q3-1) is special case
- if ((q3 <= 19 && C3.w[0] != __bid_ten2k64[q3 - 1]) ||
- (q3 == 20 && (C3.w[1] != 0 || C3.w[0] != __bid_ten2k64[19])) ||
- (q3 >= 21 && (C3.w[1] != __bid_ten2k128[q3 - 21].w[1] ||
- C3.w[0] != __bid_ten2k128[q3 - 21].w[0]))) {
- // C3 * 10^ scale != 10^(q3-1)
+ if ((q3 <= 19 && C3.w[0] != ten2k64[q3 - 1]) ||
+ (q3 == 20 && (C3.w[1] != 0 || C3.w[0] != ten2k64[19])) ||
+ (q3 >= 21 && (C3.w[1] != ten2k128[q3 - 21].w[1] ||
+ C3.w[0] != ten2k128[q3 - 21].w[0]))) {
+ // C3 * 10^ scale != 10^(q3-1)
// if ((res.w[1] & MASK_COEFF) != 0x0000314dc6448d93ull ||
// res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33
is_inexact_gt_midpoint = 1; // if (z_sign), set as if for abs. value
// if q4 > 1 then truncate C4 from q4 digits to 1 digit;
// x = q4-1, 1 <= x <= 67 and check if this operation is exact
if (q4 <= 18) { // 2 <= q4 <= 18
- __bid_round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
} else if (q4 <= 38) {
P128.w[1] = C4.w[1];
P128.w[0] = C4.w[0];
- __bid_round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
R64 = R128.w[0]; // one decimal digit
} else if (q4 <= 57) {
P192.w[2] = C4.w[2];
P192.w[1] = C4.w[1];
P192.w[0] = C4.w[0];
- __bid_round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
R64 = R192.w[0]; // one decimal digit
} else { // if (q4 <= 68)
- __bid_round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
R64 = R256.w[0]; // one decimal digit
}
if (incr_exp) {
is_midpoint_lt_even = 1;
is_midpoint_gt_even = 0;
} else if ((e3 == expmin) ||
- R64 < 5 || (R64 == 5 && is_inexact_gt_midpoint)) {
+ R64 < 5 || (R64 == 5 && is_inexact_gt_midpoint)) {
// result does not change
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
// result decremented is 10^(q3+scale) - 1
if ((q3 + scale) <= 19) {
res.w[1] = 0;
- res.w[0] = __bid_ten2k64[q3 + scale];
+ res.w[0] = ten2k64[q3 + scale];
} else { // if ((q3 + scale + 1) <= 35)
- res.w[1] = __bid_ten2k128[q3 + scale - 20].w[1];
- res.w[0] = __bid_ten2k128[q3 + scale - 20].w[0];
+ res.w[1] = ten2k128[q3 + scale - 20].w[1];
+ res.w[0] = ten2k128[q3 + scale - 20].w[0];
}
res.w[0] = res.w[0] - 1; // borrow never occurs
z_exp = z_exp - EXP_P1;
*pfpsf |= UNDERFLOW_EXCEPTION;
}
}
- } // end 10^(q3+scale-1)
+ } // end 10^(q3+scale-1)
// set the inexact flag
*pfpsf |= INEXACT_EXCEPTION;
} else {
// endif
// endif
// 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) ||
+ 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)))) {
*pfpsf |= UNDERFLOW_EXCEPTION;
}
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e3, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e3, &res, pfpsf);
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if (p34 == delta) { // Case (1''B)
res.w[0] = C3.w[0];
} else if (q3 <= 19) { // z fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
}
} else { // z fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
}
// subtract scale from the exponent
z_exp = z_exp - ((UINT64) scale << 49);
// determine whether x * y is less than, equal to, or greater than
// 1/2 ulp (z)
if (q4 <= 19) {
- if (C4.w[0] < __bid_midpoint64[q4 - 1]) { // < 1/2 ulp
+ if (C4.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
lt_half_ulp = 1;
- } else if (C4.w[0] == __bid_midpoint64[q4 - 1]) { // = 1/2 ulp
+ } else if (C4.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
eq_half_ulp = 1;
} else { // > 1/2 ulp
gt_half_ulp = 1;
}
} else if (q4 <= 38) {
- if (C4.w[2] == 0 && (C4.w[1] < __bid_midpoint128[q4 - 20].w[1] ||
- (C4.w[1] == __bid_midpoint128[q4 - 20].w[1] &&
- C4.w[0] < __bid_midpoint128[q4 - 20].w[0]))) { // < 1/2 ulp
+ if (C4.w[2] == 0 && (C4.w[1] < midpoint128[q4 - 20].w[1] ||
+ (C4.w[1] == midpoint128[q4 - 20].w[1] &&
+ C4.w[0] < midpoint128[q4 - 20].w[0]))) { // < 1/2 ulp
lt_half_ulp = 1;
- } else if (C4.w[2] == 0 && C4.w[1] == __bid_midpoint128[q4 - 20].w[1] &&
- C4.w[0] == __bid_midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
+ } else if (C4.w[2] == 0 && C4.w[1] == midpoint128[q4 - 20].w[1] &&
+ C4.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
} else { // > 1/2 ulp
gt_half_ulp = 1;
}
} else if (q4 <= 58) {
- if (C4.w[3] == 0 && (C4.w[2] < __bid_midpoint192[q4 - 39].w[2] ||
- (C4.w[2] == __bid_midpoint192[q4 - 39].w[2] &&
- C4.w[1] < __bid_midpoint192[q4 - 39].w[1]) ||
- (C4.w[2] == __bid_midpoint192[q4 - 39].w[2] &&
- C4.w[1] == __bid_midpoint192[q4 - 39].w[1] &&
- C4.w[0] < __bid_midpoint192[q4 - 39].w[0]))) { // < 1/2 ulp
+ if (C4.w[3] == 0 && (C4.w[2] < midpoint192[q4 - 39].w[2] ||
+ (C4.w[2] == midpoint192[q4 - 39].w[2] &&
+ C4.w[1] < midpoint192[q4 - 39].w[1]) ||
+ (C4.w[2] == midpoint192[q4 - 39].w[2] &&
+ C4.w[1] == midpoint192[q4 - 39].w[1] &&
+ C4.w[0] < midpoint192[q4 - 39].w[0]))) { // < 1/2 ulp
lt_half_ulp = 1;
- } else if (C4.w[3] == 0 && C4.w[2] == __bid_midpoint192[q4 - 39].w[2] &&
- C4.w[1] == __bid_midpoint192[q4 - 39].w[1] &&
- C4.w[0] == __bid_midpoint192[q4 - 39].w[0]) { // = 1/2 ulp
+ } else if (C4.w[3] == 0 && C4.w[2] == midpoint192[q4 - 39].w[2] &&
+ C4.w[1] == midpoint192[q4 - 39].w[1] &&
+ C4.w[0] == midpoint192[q4 - 39].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
} else { // > 1/2 ulp
gt_half_ulp = 1;
}
} else {
- if (C4.w[3] < __bid_midpoint256[q4 - 59].w[3] ||
- (C4.w[3] == __bid_midpoint256[q4 - 59].w[3] &&
- C4.w[2] < __bid_midpoint256[q4 - 59].w[2]) ||
- (C4.w[3] == __bid_midpoint256[q4 - 59].w[3] &&
- C4.w[2] == __bid_midpoint256[q4 - 59].w[2] &&
- C4.w[1] < __bid_midpoint256[q4 - 59].w[1]) ||
- (C4.w[3] == __bid_midpoint256[q4 - 59].w[3] &&
- C4.w[2] == __bid_midpoint256[q4 - 59].w[2] &&
- C4.w[1] == __bid_midpoint256[q4 - 59].w[1] &&
- C4.w[0] < __bid_midpoint256[q4 - 59].w[0])) { // < 1/2 ulp
+ if (C4.w[3] < midpoint256[q4 - 59].w[3] ||
+ (C4.w[3] == midpoint256[q4 - 59].w[3] &&
+ C4.w[2] < midpoint256[q4 - 59].w[2]) ||
+ (C4.w[3] == midpoint256[q4 - 59].w[3] &&
+ C4.w[2] == midpoint256[q4 - 59].w[2] &&
+ C4.w[1] < midpoint256[q4 - 59].w[1]) ||
+ (C4.w[3] == midpoint256[q4 - 59].w[3] &&
+ C4.w[2] == midpoint256[q4 - 59].w[2] &&
+ C4.w[1] == midpoint256[q4 - 59].w[1] &&
+ C4.w[0] < midpoint256[q4 - 59].w[0])) { // < 1/2 ulp
lt_half_ulp = 1;
- } else if (C4.w[3] == __bid_midpoint256[q4 - 59].w[3] &&
- C4.w[2] == __bid_midpoint256[q4 - 59].w[2] &&
- C4.w[1] == __bid_midpoint256[q4 - 59].w[1] &&
- C4.w[0] == __bid_midpoint256[q4 - 59].w[0]) { // = 1/2 ulp
+ } else if (C4.w[3] == midpoint256[q4 - 59].w[3] &&
+ C4.w[2] == midpoint256[q4 - 59].w[2] &&
+ C4.w[1] == midpoint256[q4 - 59].w[1] &&
+ C4.w[0] == midpoint256[q4 - 59].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
} else { // > 1/2 ulp
gt_half_ulp = 1;
res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
res.w[0] = 0x0000000000000000ull;
*pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e3, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e3, &res, pfpsf);
+ z_exp = res.w[1] & MASK_EXP;
}
} else { // if (p_sign != z_sign)
// consider two cases, because C3 * 10^scale = 10^33 is a special case
*pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
} else {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even,
- is_midpoint_gt_even, e3, &res,
- pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even,
+ is_midpoint_gt_even, e3, &res,
+ pfpsf);
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
// set the inexact flag
*pfpsf |= INEXACT_EXCEPTION;
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even,
- is_midpoint_gt_even, e3, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even,
+ is_midpoint_gt_even, e3, &res, pfpsf);
}
+ z_exp = res.w[1] & MASK_EXP;
} else { // if C3 * 10^scale = 10^33
e3 = (z_exp >> 49) - 6176;
if (e3 > expmin) {
// if q4 > 1 then truncate C4 from q4 digits to 1 digit;
// x = q4-1, 1 <= x <= 67 and check if this operation is exact
if (q4 <= 18) { // 2 <= q4 <= 18
- __bid_round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
} else if (q4 <= 38) {
- P128.w[1] = C4.w[1];
- P128.w[0] = C4.w[0];
- __bid_round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
- R64 = R128.w[0]; // one decimal digit
+ P128.w[1] = C4.w[1];
+ P128.w[0] = C4.w[0];
+ round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
+ R64 = R128.w[0]; // one decimal digit
} else if (q4 <= 57) {
- P192.w[2] = C4.w[2];
- P192.w[1] = C4.w[1];
- P192.w[0] = C4.w[0];
- __bid_round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
- R64 = R192.w[0]; // one decimal digit
+ P192.w[2] = C4.w[2];
+ P192.w[1] = C4.w[1];
+ P192.w[0] = C4.w[0];
+ round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
+ R64 = R192.w[0]; // one decimal digit
} else { // if (q4 <= 68)
- __bid_round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
- &is_midpoint_lt_even,
- &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
- R64 = R256.w[0]; // one decimal digit
+ round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
+ &is_midpoint_lt_even,
+ &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
+ R64 = R256.w[0]; // one decimal digit
}
if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
- // the result is exact: 10^34 - R64
- // incr_exp = 0 with certainty
- z_exp = z_exp - EXP_P1;
- e3 = e3 - 1;
- res.w[1] =
- z_sign | (z_exp & MASK_EXP) | 0x0001ed09bead87c0ull;
- res.w[0] = 0x378d8e6400000000ull - R64;
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ // the result is exact: 10^34 - R64
+ // incr_exp = 0 with certainty
+ z_exp = z_exp - EXP_P1;
+ e3 = e3 - 1;
+ res.w[1] =
+ z_sign | (z_exp & MASK_EXP) | 0x0001ed09bead87c0ull;
+ res.w[0] = 0x378d8e6400000000ull - R64;
} else {
- // We want R64 to be the top digit of C4, but we actually
- // obtained (C4 * 10^(-q4+1))RN; a correction may be needed,
- // because the top digit is (C4 * 10^(-q4+1))RZ
- // however, if incr_exp = 1 then R64 = 10 with certainty
- if (incr_exp) {
- R64 = 10;
- }
- // the result is inexact as C4 has more than 1 significant digit
- // and C3 * 10^scale = 10^33
- // example of case that is treated here:
- // 100...0 * 10^e3 - 0.41 * 10^e3 =
- // 0999...9.59 * 10^e3 -> rounds to 99...96*10^(e3-1)
- // note that (e3 > expmin}
- // in order to round, subtract R64 from 10^34 and then compare
- // C4 - R64 * 10^(q4-1) with 1/2 ulp
- // calculate 10^34 - R64
- res.w[1] = 0x0001ed09bead87c0ull;
- res.w[0] = 0x378d8e6400000000ull - R64;
- z_exp = z_exp - EXP_P1; // will be OR-ed with sign & significand
- // calculate C4 - R64 * 10^(q4-1); this is a rare case and
- // R64 is small, 1 <= R64 <= 9
- e3 = e3 - 1;
- if (is_inexact_lt_midpoint) {
- is_inexact_lt_midpoint = 0;
- is_inexact_gt_midpoint = 1;
- } else if (is_inexact_gt_midpoint) {
- is_inexact_gt_midpoint = 0;
- is_inexact_lt_midpoint = 1;
- } else if (is_midpoint_lt_even) {
- is_midpoint_lt_even = 0;
- is_midpoint_gt_even = 1;
- } else if (is_midpoint_gt_even) {
- is_midpoint_gt_even = 0;
- is_midpoint_lt_even = 1;
- } else {
- ;
- }
- // the result is always inexact, and never tiny
- // check for overflow for RN
- if (e3 > expmax) {
- if (rnd_mode == ROUNDING_TO_NEAREST) {
- res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
- } else {
- rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even,
- is_midpoint_gt_even, e3, &res,
- pfpsf);
- }
- BID_RETURN (res)
- }
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- res.w[1] =
- z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
- if (rnd_mode != ROUNDING_TO_NEAREST) {
- rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even,
- is_midpoint_gt_even, e3, &res,
- pfpsf);
- }
- z_exp = res.w[1] & MASK_EXP;
- e3 = (z_exp >> 49) - 6176;
- } // end result is inexact
- } // end q4 > 1
+ // We want R64 to be the top digit of C4, but we actually
+ // obtained (C4 * 10^(-q4+1))RN; a correction may be needed,
+ // because the top digit is (C4 * 10^(-q4+1))RZ
+ // however, if incr_exp = 1 then R64 = 10 with certainty
+ if (incr_exp) {
+ R64 = 10;
+ }
+ // the result is inexact as C4 has more than 1 significant digit
+ // and C3 * 10^scale = 10^33
+ // example of case that is treated here:
+ // 100...0 * 10^e3 - 0.41 * 10^e3 =
+ // 0999...9.59 * 10^e3 -> rounds to 99...96*10^(e3-1)
+ // note that (e3 > expmin}
+ // in order to round, subtract R64 from 10^34 and then compare
+ // C4 - R64 * 10^(q4-1) with 1/2 ulp
+ // calculate 10^34 - R64
+ res.w[1] = 0x0001ed09bead87c0ull;
+ res.w[0] = 0x378d8e6400000000ull - R64;
+ z_exp = z_exp - EXP_P1; // will be OR-ed with sign & significand
+ // calculate C4 - R64 * 10^(q4-1); this is a rare case and
+ // R64 is small, 1 <= R64 <= 9
+ e3 = e3 - 1;
+ if (is_inexact_lt_midpoint) {
+ is_inexact_lt_midpoint = 0;
+ is_inexact_gt_midpoint = 1;
+ } else if (is_inexact_gt_midpoint) {
+ is_inexact_gt_midpoint = 0;
+ is_inexact_lt_midpoint = 1;
+ } else if (is_midpoint_lt_even) {
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 1;
+ } else if (is_midpoint_gt_even) {
+ is_midpoint_gt_even = 0;
+ is_midpoint_lt_even = 1;
+ } else {
+ ;
+ }
+ // the result is always inexact, and never tiny
+ // check for overflow for RN
+ if (e3 > expmax) {
+ if (rnd_mode == ROUNDING_TO_NEAREST) {
+ res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+ res.w[0] = 0x0000000000000000ull;
+ *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+ } else {
+ rounding_correction (rnd_mode,
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even,
+ is_midpoint_gt_even, e3, &res,
+ pfpsf);
+ }
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
+ BID_RETURN (res)
+ }
+ // set the inexact flag
+ *pfpsf |= INEXACT_EXCEPTION;
+ res.w[1] =
+ z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
+ if (rnd_mode != ROUNDING_TO_NEAREST) {
+ rounding_correction (rnd_mode,
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even,
+ is_midpoint_gt_even, e3, &res,
+ pfpsf);
+ }
+ z_exp = res.w[1] & MASK_EXP;
+ } // end result is inexact
+ } // end q4 > 1
} else { // if (e3 = emin)
// if e3 = expmin the result is also tiny (the condition for
// tininess is C4 > 050...0 [q4 digits] which is met because
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even,
- is_midpoint_gt_even, e3, &res,
- pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even,
+ is_midpoint_gt_even, e3, &res,
+ pfpsf);
+ z_exp = res.w[1] & MASK_EXP;
}
- } // end e3 = emin
+ } // end e3 = emin
// set the inexact flag (if the result was not exact)
if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
is_midpoint_lt_even || is_midpoint_gt_even)
*pfpsf |= INEXACT_EXCEPTION;
- } // end 10^33
- } // end if (p_sign != z_sign)
+ } // end 10^33
+ } // end if (p_sign != z_sign)
res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if (((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2)
(q3 <= delta && delta + q4 <= p34) || // Case (3)
(delta < q3 && p34 < delta + q4) || // Case (4)
(delta < q3 && q3 <= delta + q4 && delta + q4 <= p34) || // Case (5)
- (delta + q4 < q3)) &&// Case (6)
+ (delta + q4 < q3)) && // Case (6)
!(delta <= 1 && p_sign != z_sign)) { // Case (2), (3), (4), (5) or (6)
// the result has the sign of z
scale = q3 - delta - q4; // 1 <= scale <= 33
if (q4 <= 19) { // 1 <= scale <= 19; C4 fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C4.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (P128, C4.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C4.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (P128, C4.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C4.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (P128, C4.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C4.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (P128, C4.w[0], ten2k128[scale - 20]);
}
} else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C4
- __mul_128x64_to_128 (P128, __bid_ten2k64[scale], C4);
+ // 64 x 128 ten2k64[scale] * C4
+ __mul_128x64_to_128 (P128, ten2k64[scale], C4);
}
C4.w[0] = P128.w[0];
C4.w[1] = P128.w[1];
res.w[0] = C3.w[0];
} else if (q3 <= 19) { // 1 <= scale <= 19; z fits in 64 bits
if (scale <= 19) { // 10^scale fits in 64 bits
- // 64 x 64 C3.w[0] * __bid_ten2k64[scale]
- __mul_64x64_to_128MACH (res, C3.w[0], __bid_ten2k64[scale]);
+ // 64 x 64 C3.w[0] * ten2k64[scale]
+ __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
} else { // 10^scale fits in 128 bits
- // 64 x 128 C3.w[0] * __bid_ten2k128[scale - 20]
- __mul_128x64_to_128 (res, C3.w[0], __bid_ten2k128[scale - 20]);
+ // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+ __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
}
} else { // z fits in 128 bits, but 10^scale must fit in 64 bits
- // 64 x 128 __bid_ten2k64[scale] * C3
- __mul_128x64_to_128 (res, __bid_ten2k64[scale], C3);
+ // 64 x 128 ten2k64[scale] * C3
+ __mul_128x64_to_128 (res, ten2k64[scale], C3);
}
// e3 is already calculated
e3 = e3 - scale;
// x0 = delta + q4 - p34 (calculated before reaching case2_repeat)
// because q3 + q4 - x0 <= P => x0 >= q3 + q4 - p34
if (x0 == 0) { // this could happen only if we return to case2_repeat, or
- // for Case (6)
+ // for Case (3) or Case (6)
R128.w[1] = C4.w[1];
R128.w[0] = C4.w[0];
} else if (q4 <= 18) {
// 2 <= q4 <= 18, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 17
- __bid_round64_2_18 (q4, x0, C4.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round64_2_18 (q4, x0, C4.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
if (incr_exp) {
// R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
- R64 = __bid_ten2k64[q4 - x0];
+ R64 = ten2k64[q4 - x0];
}
R128.w[1] = 0;
R128.w[0] = R64;
// 19 <= q4 <= 38, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 37
P128.w[1] = C4.w[1];
P128.w[0] = C4.w[0];
- __bid_round128_19_38 (q4, x0, P128, &R128, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q4, x0, P128, &R128, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
if (incr_exp) {
// R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
- R128.w[0] = __bid_ten2k64[q4 - x0];
+ R128.w[0] = ten2k64[q4 - x0];
// R128.w[1] stays 0
} else { // 20 <= q4 - x0 <= 37
- R128.w[0] = __bid_ten2k128[q4 - x0 - 20].w[0];
- R128.w[1] = __bid_ten2k128[q4 - x0 - 20].w[1];
+ R128.w[0] = ten2k128[q4 - x0 - 20].w[0];
+ R128.w[1] = ten2k128[q4 - x0 - 20].w[1];
}
}
} else if (q4 <= 57) {
P192.w[2] = C4.w[2];
P192.w[1] = C4.w[1];
P192.w[0] = C4.w[0];
- __bid_round192_39_57 (q4, x0, P192, &R192, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
// R192.w[2] is always 0
if (incr_exp) {
// R192 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 5, 1 <= q4 - x0 <= 52
if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
- R192.w[0] = __bid_ten2k64[q4 - x0];
+ R192.w[0] = ten2k64[q4 - x0];
// R192.w[1] stays 0
// R192.w[2] stays 0
} else { // 20 <= q4 - x0 <= 33
- R192.w[0] = __bid_ten2k128[q4 - x0 - 20].w[0];
- R192.w[1] = __bid_ten2k128[q4 - x0 - 20].w[1];
+ R192.w[0] = ten2k128[q4 - x0 - 20].w[0];
+ R192.w[1] = ten2k128[q4 - x0 - 20].w[1];
// R192.w[2] stays 0
}
}
R128.w[0] = R192.w[0];
} else {
// 58 <= q4 <= 68, max(1, q3+q4-p34) <= x0 <= q4 - 1, 25 <= x0 <= 67
- __bid_round256_58_76 (q4, x0, C4, &R256, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
// R256.w[3] and R256.w[2] are always 0
if (incr_exp) {
// R256 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 25, 1 <= q4 - x0 <= 43
if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
- R256.w[0] = __bid_ten2k64[q4 - x0];
+ R256.w[0] = ten2k64[q4 - x0];
// R256.w[1] stays 0
// R256.w[2] stays 0
// R256.w[3] stays 0
} else { // 20 <= q4 - x0 <= 33
- R256.w[0] = __bid_ten2k128[q4 - x0 - 20].w[0];
- R256.w[1] = __bid_ten2k128[q4 - x0 - 20].w[1];
+ R256.w[0] = ten2k128[q4 - x0 - 20].w[0];
+ R256.w[1] = ten2k128[q4 - x0 - 20].w[1];
// R256.w[2] stays 0
// R256.w[3] stays 0
}
// if res > 10^34 - 1 need to increase x0 and decrease scale by 1
if (res.w[1] > 0x0001ed09bead87c0ull ||
(res.w[1] == 0x0001ed09bead87c0ull &&
- res.w[0] > 0x378d8e63ffffffffull)) {
+ res.w[0] > 0x378d8e63ffffffffull)) {
// avoid double rounding error
is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
is_midpoint_gt_even = 0;
P128.w[1] = res.w[1];
P128.w[0] = res.w[0];
- __bid_round128_19_38 (35, 1, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (35, 1, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
// incr_exp is 0 with certainty in this case
// avoid a double rounding error
if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint
- && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint
+ && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0
+ || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0
+ || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
!is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
if (is_midpoint_lt_even0 || is_midpoint_gt_even0 ||
- is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
+ is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
is_inexact_lt_midpoint = 1;
}
}
is_midpoint_lt_even = 1;
res.w[0]++;
if (res.w[0] == 0x0)
- res.w[1]++;
+ res.w[1]++;
// check for rounding overflow
if (res.w[1] == 0x0001ed09bead87c0ull &&
- res.w[0] == 0x378d8e6400000000ull) {
- // res = 10^34 => rounding overflow
- res.w[1] = 0x0000314dc6448d93ull;
- res.w[0] = 0x38c15b0a00000000ull; // 10^33
- e3++;
+ res.w[0] == 0x378d8e6400000000ull) {
+ // res = 10^34 => rounding overflow
+ res.w[1] = 0x0000314dc6448d93ull;
+ res.w[0] = 0x38c15b0a00000000ull; // 10^33
+ e3++;
}
} else if (is_midpoint_lt_even) {
// res = res - 1
is_midpoint_gt_even = 1;
res.w[0]--;
if (res.w[0] == 0xffffffffffffffffull)
- res.w[1]--;
+ res.w[1]--;
// if the result is pure zero, the sign depends on the rounding
// mode (x*y and z had opposite signs)
if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
- if (rnd_mode != ROUNDING_DOWN)
- z_sign = 0x0000000000000000ull;
- else
- z_sign = 0x8000000000000000ull;
- // the exponent is max (e3, expmin)
- res.w[1] = 0x0;
- res.w[0] = 0x0;
- BID_RETURN (res)
+ if (rnd_mode != ROUNDING_DOWN)
+ z_sign = 0x0000000000000000ull;
+ else
+ z_sign = 0x8000000000000000ull;
+ // the exponent is max (e3, expmin)
+ res.w[1] = 0x0;
+ res.w[0] = 0x0;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
+ BID_RETURN (res)
}
} else {
;
res.w[1]--; // borrow
// if res < 10^33 and exp > expmin need to decrease x0 and
// increase scale by 1
- if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull ||
- (res.w[1] == 0x0000314dc6448d93ull &&
- res.w[0] < 0x38c15b0a00000000ull)) ||
- (is_inexact_lt_midpoint && res.w[1] == 0x0000314dc6448d93ull &&
- res.w[0] == 0x38c15b0a00000000ull)) && x0 >= 1) {
+ if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull ||
+ (res.w[1] == 0x0000314dc6448d93ull &&
+ res.w[0] < 0x38c15b0a00000000ull)) ||
+ (is_inexact_lt_midpoint
+ && res.w[1] == 0x0000314dc6448d93ull
+ && res.w[0] == 0x38c15b0a00000000ull))
+ && x0 >= 1) {
x0 = x0 - 1;
// first restore e3, otherwise it will be too small
e3 = e3 + scale;
incr_exp = 0;
goto case2_repeat;
}
- // else this is the result rounded with unbounded exponent
+ // else this is the result rounded with unbounded exponent;
// because the result has opposite sign to that of C4 which was
// rounded, need to change the rounding indicators
if (is_inexact_lt_midpoint) {
res.w[1]++;
// check for rounding overflow
if (res.w[1] == 0x0001ed09bead87c0ull &&
- res.w[0] == 0x378d8e6400000000ull) {
+ res.w[0] == 0x378d8e6400000000ull) {
// res = 10^34 => rounding overflow
res.w[1] = 0x0000314dc6448d93ull;
res.w[0] = 0x38c15b0a00000000ull; // 10^33
// mode (x*y and z had opposite signs)
if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
if (rnd_mode != ROUNDING_DOWN)
- z_sign = 0x0000000000000000ull;
+ z_sign = 0x0000000000000000ull;
else
- z_sign = 0x8000000000000000ull;
+ z_sign = 0x8000000000000000ull;
// the exponent is max (e3, expmin)
res.w[1] = 0x0;
res.w[0] = 0x0;
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
} else {
}
}
// check for underflow
- if (e3 < expmin) {
+ if (e3 == expmin) { // and if significand < 10^33 => result is tiny
+ if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
+ ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
+ res.w[0] < 0x38c15b0a00000000ull)) {
+ is_tiny = 1;
+ }
+ } else if (e3 < expmin) {
// the result is tiny, so we must truncate more of res
is_tiny = 1;
x0 = expmin - e3;
if (res.w[1] == 0x0) {
// between 1 and 19 digits
for (ind = 1; ind <= 19; ind++) {
- if (res.w[0] < __bid_ten2k64[ind]) {
+ if (res.w[0] < ten2k64[ind]) {
break;
}
}
// ind digits
- } else if ((res.w[1] < __bid_ten2k128[0].w[1] ||
- (res.w[1] == __bid_ten2k128[0].w[1]
- && res.w[0] < __bid_ten2k128[0].w[0]))) {
+ } else if (res.w[1] < ten2k128[0].w[1] ||
+ (res.w[1] == ten2k128[0].w[1]
+ && res.w[0] < ten2k128[0].w[0])) {
// 20 digits
ind = 20;
} else { // between 21 and 38 digits
for (ind = 1; ind <= 18; ind++) {
- if (res.w[1] < __bid_ten2k128[ind].w[1] ||
- (res.w[1] == __bid_ten2k128[ind].w[1] &&
- res.w[0] < __bid_ten2k128[ind].w[0])) {
+ if (res.w[1] < ten2k128[ind].w[1] ||
+ (res.w[1] == ten2k128[ind].w[1] &&
+ res.w[0] < ten2k128[ind].w[0])) {
break;
}
}
}
// at this point ind >= x0; because delta >= 2 on this path, the case
- // ind = x0 cab occur only in Case (2) or case (3), when C3 has one
+ // ind = x0 can occur only in Case (2) or case (3), when C3 has one
// digit (q3 = 1) equal to 1 (C3 = 1), e3 is expmin (e3 = expmin),
// the signs of x * y and z are opposite, and through cancellation
// the most significant decimal digit in res has the weight
is_inexact_gt_midpoint = 1;
} else if (ind <= 18) { // check that 2 <= ind
// 2 <= ind <= 18, 1 <= x0 <= 17
- __bid_round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
if (incr_exp) {
- // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
- R64 = __bid_ten2k64[q4 - x0];
+ // R64 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 17
+ R64 = ten2k64[ind - x0];
}
res.w[1] = 0;
res.w[0] = R64;
// 19 <= ind <= 38
P128.w[1] = res.w[1];
P128.w[0] = res.w[0];
- __bid_round128_19_38 (ind, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (ind, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
if (incr_exp) {
- // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
+ // R128 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 37
if (ind - x0 <= 19) { // 1 <= ind - x0 <= 19
- res.w[0] = __bid_ten2k64[ind - x0];
+ res.w[0] = ten2k64[ind - x0];
// res.w[1] stays 0
} else { // 20 <= ind - x0 <= 37
- res.w[0] = __bid_ten2k128[ind - x0 - 20].w[0];
- res.w[1] = __bid_ten2k128[ind - x0 - 20].w[1];
+ res.w[0] = ten2k128[ind - x0 - 20].w[0];
+ res.w[1] = ten2k128[ind - x0 - 20].w[1];
}
}
}
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
is_inexact_lt_midpoint = 1;
}
}
+ } else {
+ ; // not underflow
}
// check for inexact result
if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
}
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e3, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e3, &res, pfpsf);
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else {
ind = e3;
e3 = e4;
e4 = ind;
- tmp64 = z_sign;
+ tmp_sign = z_sign;
z_sign = p_sign;
- p_sign = tmp64;
+ p_sign = tmp_sign;
} else { // from Cases (2), (3), (4), (5)
// In Cases (2), (3), (4), (5) with 0 <= delta <= 1 C3 has to be
// scaled up by q4 + delta - q3; this is the same as in Cases (15),
delta = -delta;
}
add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
- rnd_mode, &is_midpoint_lt_even,
- &is_midpoint_gt_even, &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint, pfpsf, &res);
+ rnd_mode, &is_midpoint_lt_even,
+ &is_midpoint_gt_even, &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, pfpsf, &res);
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
if (q4 <= 38) {
P128.w[1] = C4.w[1];
P128.w[0] = C4.w[0];
- __bid_round128_19_38 (q4, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q4, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
} else if (q4 <= 57) { // 35 <= q4 <= 57
P192.w[2] = C4.w[2];
P192.w[1] = C4.w[1];
P192.w[0] = C4.w[0];
- __bid_round192_39_57 (q4, x0, P192, &R192, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
res.w[0] = R192.w[0];
res.w[1] = R192.w[1];
} else { // if (q4 <= 68)
- __bid_round256_58_76 (q4, x0, C4, &R256, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
res.w[0] = R256.w[0];
res.w[1] = R256.w[1];
}
is_inexact_gt_midpoint = 1;
} else { // if (delta == p34 + 1)
if (q3 <= 19) {
- if (C3.w[0] < __bid_midpoint64[q3 - 1]) { // C3 < 1/2 ulp
- // res = 10^33, unchanged
- is_inexact_gt_midpoint = 1;
- } else if (C3.w[0] == __bid_midpoint64[q3 - 1]) { // C3 = 1/2 ulp
- // res = 10^33, unchanged
- is_midpoint_lt_even = 1;
- } else { // if (C3.w[0] > __bid_midpoint64[q3-1]), C3 > 1/2 ulp
- res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
- res.w[0] = 0x378d8e63ffffffffull;
- e4 = e4 - 1;
- is_inexact_lt_midpoint = 1;
- }
+ if (C3.w[0] < midpoint64[q3 - 1]) { // C3 < 1/2 ulp
+ // res = 10^33, unchanged
+ is_inexact_gt_midpoint = 1;
+ } else if (C3.w[0] == midpoint64[q3 - 1]) { // C3 = 1/2 ulp
+ // res = 10^33, unchanged
+ is_midpoint_lt_even = 1;
+ } else { // if (C3.w[0] > midpoint64[q3-1]), C3 > 1/2 ulp
+ res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+ res.w[0] = 0x378d8e63ffffffffull;
+ e4 = e4 - 1;
+ is_inexact_lt_midpoint = 1;
+ }
} else { // if (20 <= q3 <=34)
- if (C3.w[1] < __bid_midpoint128[q3 - 20].w[1] ||
- (C3.w[1] == __bid_midpoint128[q3 - 20].w[1] &&
- C3.w[0] < __bid_midpoint128[q3 - 20].w[0])) { // C3 < 1/2 ulp
- // res = 10^33, unchanged
- is_inexact_gt_midpoint = 1;
- } else if (C3.w[1] == __bid_midpoint128[q3 - 20].w[1] &&
- C3.w[0] == __bid_midpoint128[q3 - 20].w[0]) { // C3 = 1/2 ulp
- // res = 10^33, unchanged
- is_midpoint_lt_even = 1;
- } else { // if (C3 > __bid_midpoint128[q3-20]), C3 > 1/2 ulp
- res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
- res.w[0] = 0x378d8e63ffffffffull;
- e4 = e4 - 1;
- is_inexact_lt_midpoint = 1;
- }
+ if (C3.w[1] < midpoint128[q3 - 20].w[1] ||
+ (C3.w[1] == midpoint128[q3 - 20].w[1] &&
+ C3.w[0] < midpoint128[q3 - 20].w[0])) { // C3 < 1/2 ulp
+ // res = 10^33, unchanged
+ is_inexact_gt_midpoint = 1;
+ } else if (C3.w[1] == midpoint128[q3 - 20].w[1] &&
+ C3.w[0] == midpoint128[q3 - 20].w[0]) { // C3 = 1/2 ulp
+ // res = 10^33, unchanged
+ is_midpoint_lt_even = 1;
+ } else { // if (C3 > midpoint128[q3-20]), C3 > 1/2 ulp
+ res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+ res.w[0] = 0x378d8e63ffffffffull;
+ e4 = e4 - 1;
+ is_inexact_lt_midpoint = 1;
+ }
}
}
}
}
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, pfpsf);
}
if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
is_midpoint_lt_even || is_midpoint_gt_even) {
// set the inexact flag
*pfpsf |= INEXACT_EXCEPTION;
}
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((q4 <= p34 && p34 <= delta) || // Case (8)
ind = e3;
e3 = e4;
e4 = ind;
- tmp64 = z_sign;
+ tmp_sign = z_sign;
z_sign = p_sign;
- p_sign = tmp64;
- tmp64 = z_exp;
+ p_sign = tmp_sign;
+ tmp.ui64 = z_exp;
z_exp = p_exp;
- p_exp = tmp64;
+ p_exp = tmp.ui64;
goto delta_ge_zero;
} else if ((p34 <= delta && delta < q4 && q4 < delta + q3) || // Case (11)
// 1 <= x0 <= q3 - 1 <= p34 - 1
x0 = e4 - e3; // or x0 = delta + q3 - q4
if (q3 <= 18) { // 2 <= q3 <= 18
- __bid_round64_2_18 (q3, x0, C3.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ round64_2_18 (q3, x0, C3.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
// C3.w[1] = 0;
C3.w[0] = R64;
} else if (q3 <= 38) {
- __bid_round128_19_38 (q3, x0, C3, &R128, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (q3, x0, C3, &R128, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
C3.w[1] = R128.w[1];
C3.w[0] = R128.w[0];
}
// 64 x 128 -> 128
P128.w[1] = C3.w[1];
P128.w[0] = C3.w[0];
- __mul_64x128_to_128 (C3, __bid_ten2k64[1], P128);
+ __mul_64x128_to_128 (C3, ten2k64[1], P128);
}
e3 = e3 + x0; // this is e4
// now add/subtract the 256-bit C4 and the new (and shorter) 128-bit C3;
if (R256.w[0] == 0x0) {
R256.w[1]++;
if (R256.w[1] == 0x0) {
- R256.w[2]++;
- if (R256.w[2] == 0x0) {
- R256.w[3]++;
- }
+ R256.w[2]++;
+ if (R256.w[2] == 0x0) {
+ R256.w[3]++;
+ }
}
}
// no check for rounding overflow - R256 was a difference
if (R256.w[0] == 0xffffffffffffffffull) {
R256.w[1]--;
if (R256.w[1] == 0xffffffffffffffffull) {
- R256.w[2]--;
- if (R256.w[2] == 0xffffffffffffffffull) {
- R256.w[3]--;
- }
+ R256.w[2]--;
+ if (R256.w[2] == 0xffffffffffffffffull) {
+ R256.w[3]--;
+ }
}
}
} else {
}
}
// determine the number of decimal digits in R256
- ind = __bid_nr_digits256 (R256); // ind >= p34
+ ind = nr_digits256 (R256); // ind >= p34
// if R256 is sum, then ind > p34; if R256 is a difference, then
// ind >= p34; this means that we can calculate the result rounded to
// the destination precision, with unbounded exponent, starting from R256
if (ind <= 38) {
P128.w[1] = R256.w[1];
P128.w[0] = R256.w[0];
- __bid_round128_19_38 (ind, x0, P128, &R128, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (ind, x0, P128, &R128, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
} else if (ind <= 57) {
P192.w[2] = R256.w[2];
P192.w[1] = R256.w[1];
P192.w[0] = R256.w[0];
- __bid_round192_39_57 (ind, x0, P192, &R192, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round192_39_57 (ind, x0, P192, &R192, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
R128.w[1] = R192.w[1];
R128.w[0] = R192.w[0];
} else { // if (ind <= 68)
- __bid_round256_58_76 (ind, x0, R256, &R256, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round256_58_76 (ind, x0, R256, &R256, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
R128.w[1] = R256.w[1];
R128.w[0] = R256.w[0];
}
// not possible in Cases (2)-(6) or (15)-(17) which may get here
// if this is 10^33 - 1 make it 10^34 - 1 and decrement exponent
if (res.w[1] == 0x0000314dc6448d93ull &&
- res.w[0] == 0x38c15b09ffffffffull) { // 10^33 - 1
+ res.w[0] == 0x38c15b09ffffffffull) { // 10^33 - 1
res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
res.w[0] = 0x378d8e63ffffffffull;
e4--;
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
P128.w[1] = p_sign | 0x3040000000000000ull | res.w[1];
P128.w[0] = res.w[0];
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- 0, &P128, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ 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) {
res.w[1] = p_sign | 0x7800000000000000ull;
res.w[0] = 0x0000000000000000ull;
*pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
- } // else not overflow or not RN, so continue
+ } // else not overflow or not RN, so continue
// from this point on this is similar to the last part of the computation
// for Cases (15), (16), (17)
R128.w[1] = res.w[1] & MASK_COEFF;
R128.w[0] = res.w[0];
if (ind <= 19) {
- if (R128.w[0] < __bid_midpoint64[ind - 1]) { // < 1/2 ulp
+ if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (R128.w[0] == __bid_midpoint64[ind - 1]) { // = 1/2 ulp
+ } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
is_inexact_gt_midpoint = 1;
}
} else { // if (ind <= 38)
- if (R128.w[1] < __bid_midpoint128[ind - 20].w[1] ||
- (R128.w[1] == __bid_midpoint128[ind - 20].w[1] &&
- R128.w[0] < __bid_midpoint128[ind - 20].w[0])) { // < 1/2 ulp
+ if (R128.w[1] < midpoint128[ind - 20].w[1] ||
+ (R128.w[1] == midpoint128[ind - 20].w[1] &&
+ R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
lt_half_ulp = 1;
is_inexact_lt_midpoint = 1;
- } else if (R128.w[1] == __bid_midpoint128[ind - 20].w[1] &&
- R128.w[0] == __bid_midpoint128[ind - 20].w[0]) { // = 1/2 ulp
+ } else if (R128.w[1] == midpoint128[ind - 20].w[1] &&
+ R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
eq_half_ulp = 1;
is_midpoint_gt_even = 1;
} else { // > 1/2 ulp
// round the ind-digit result to ind - x0 digits
if (ind <= 18) { // 2 <= ind <= 18
- __bid_round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
res.w[1] = 0x0;
res.w[0] = R64;
} else if (ind <= 38) {
P128.w[1] = res.w[1] & MASK_COEFF;
P128.w[0] = res.w[0];
- __bid_round128_19_38 (ind, x0, P128, &res, &incr_exp,
- &is_midpoint_lt_even, &is_midpoint_gt_even,
- &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint);
+ round128_19_38 (ind, x0, P128, &res, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint);
}
e4 = e4 + x0; // expmin
// we want the exponent to be expmin, so if incr_exp = 1 then
// 64 x 128 -> 128
P128.w[1] = res.w[1] & MASK_COEFF;
P128.w[0] = res.w[0];
- __mul_64x128_to_128 (res, __bid_ten2k64[1], P128);
+ __mul_64x128_to_128 (res, ten2k64[1], P128);
}
res.w[1] =
- p_sign | ((UINT64) (e4 + 6176) << 49) | (res.w[1] & MASK_COEFF);
+ p_sign | ((UINT64) (e4 + 6176) << 49) | (res.
+ w[1] & MASK_COEFF);
// avoid a double rounding error
if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
- is_midpoint_lt_even) { // double rounding error upward
+ is_midpoint_lt_even) { // double rounding error upward
// res = res - 1
res.w[0]--;
if (res.w[0] == 0xffffffffffffffffull)
is_midpoint_lt_even = 0;
is_inexact_lt_midpoint = 1;
} else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
- is_midpoint_gt_even) { // double rounding error downward
+ is_midpoint_gt_even) { // double rounding error downward
// res = res + 1
res.w[0]++;
if (res.w[0] == 0)
is_midpoint_gt_even = 0;
is_inexact_gt_midpoint = 1;
} else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
- !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ !is_inexact_lt_midpoint
+ && !is_inexact_gt_midpoint) {
// if this second rounding was exact the result may still be
// inexact because of the first rounding
if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
is_inexact_lt_midpoint = 1;
}
} else if (is_midpoint_gt_even &&
- (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ (is_inexact_gt_midpoint0
+ || is_midpoint_lt_even0)) {
// pulled up to a midpoint
is_inexact_lt_midpoint = 1;
is_inexact_gt_midpoint = 0;
is_midpoint_lt_even = 0;
is_midpoint_gt_even = 0;
} else if (is_midpoint_lt_even &&
- (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ (is_inexact_lt_midpoint0
+ || is_midpoint_gt_even0)) {
// pulled down to a midpoint
is_inexact_lt_midpoint = 0;
is_inexact_gt_midpoint = 1;
// apply correction if not rounding to nearest
if (rnd_mode != ROUNDING_TO_NEAREST) {
rounding_correction (rnd_mode,
- is_inexact_lt_midpoint,
- is_inexact_gt_midpoint,
- is_midpoint_lt_even, is_midpoint_gt_even,
- e4, &res, pfpsf);
+ is_inexact_lt_midpoint,
+ is_inexact_gt_midpoint,
+ is_midpoint_lt_even, is_midpoint_gt_even,
+ e4, &res, pfpsf);
}
if (is_midpoint_lt_even || is_midpoint_gt_even ||
is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
if (is_tiny)
*pfpsf |= UNDERFLOW_EXCEPTION;
}
-
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else if ((p34 <= delta && delta + q3 <= q4) || // Case (15)
// unbounded exponent
add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
- rnd_mode, &is_midpoint_lt_even,
- &is_midpoint_gt_even, &is_inexact_lt_midpoint,
- &is_inexact_gt_midpoint, pfpsf, &res);
-
+ rnd_mode, &is_midpoint_lt_even,
+ &is_midpoint_gt_even, &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, pfpsf, &res);
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
} else {
;
}
- } // end if delta < 0
+ } // end if delta < 0
+ *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+ *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+ *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+ *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+ BID_SWAP128 (res);
BID_RETURN (res)
}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT128 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128_fma (UINT128 x, UINT128 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even, is_midpoint_gt_even,
+ is_inexact_lt_midpoint, is_inexact_gt_midpoint;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, &x, &y, &z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x, y,
+ z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128ddd_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128ddd_fma (UINT64 x, UINT64 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 x1, y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, &x1, &y1, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x1, y1,
+ z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128ddq_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px, y = *py;
+ UINT128 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128ddq_fma (UINT64 x, UINT64 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 x1, y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, &x1, &y1, &z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x1, y1,
+ z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128dqd_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128dqd_fma (UINT64 x, UINT128 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 x1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, &x1, py, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x1, y,
+ z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128dqq_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128dqq_fma (UINT64 x, UINT128 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 x1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, &x1, py, pz
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x1, y,
+ z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qdd_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qdd_fma (UINT128 x, UINT64 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, px, &y1, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x, y1,
+ z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qdq_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qdq_fma (UINT128 x, UINT64 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, px, &y1, pz
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x, y1,
+ z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qqd_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qqd_fma (UINT128 x, UINT128 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+ &res, px, py, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint,
+ &is_inexact_gt_midpoint, x, y,
+ z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res);
+}
+
+// Note: bid128qqq_fma is represented by bid128_fma
+
+// Note: bid64ddd_fma is represented by bid64_fma
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64ddq_fma (UINT64 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px, y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64ddq_fma (UINT64 x, UINT64 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 x1, y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, &x1, &y1, pz
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x1, y1, z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64dqd_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64dqd_fma (UINT64 x, UINT128 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 x1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, &x1, py, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x1, y, z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64dqq_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x = *px;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64dqq_fma (UINT64 x, UINT128 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 x1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, &x1, py, pz
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x1, y, z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qdd_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qdd_fma (UINT128 x, UINT64 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, px, &y1, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x, y1, z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qdq_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qdq_fma (UINT128 x, UINT64 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, px, &y1, pz
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x, y1, z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qqd_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qqd_fma (UINT128 x, UINT128 y, UINT64 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ UINT128 z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ bid64qqq_fma (&res1, px, py, &z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+ res1 = bid64qqq_fma (x, y, z1
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+ BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qqq_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT128 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+ unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z
+ _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+#endif
+ int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0,
+ is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0;
+ int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+ is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+ int incr_exp;
+ UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT128 res128 = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+ UINT64 res1 = 0xbaddbaddbaddbaddull;
+ unsigned int save_fpsf; // needed because of the call to bid128_ext_fma
+ UINT64 sign;
+ UINT64 exp;
+ int unbexp;
+ UINT128 C;
+ BID_UI64DOUBLE tmp;
+ int nr_bits;
+ int q, x0;
+ int scale;
+ int lt_half_ulp = 0, eq_half_ulp = 0;
+
+ // Note: for rounding modes other than RN or RA, the result can be obtained
+ // by rounding first to BID128 and then to BID64
+
+ save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
+ *pfpsf = 0;
+
+#if DECIMAL_CALL_BY_REFERENCE
+ bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
+ &is_inexact_lt_midpoint0, &is_inexact_gt_midpoint0,
+ &res, &x, &y, &z
+ _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#else
+ res = bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
+ &is_inexact_lt_midpoint0,
+ &is_inexact_gt_midpoint0, x, y,
+ z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+ _EXC_INFO_ARG);
+#endif
+
+ if ((rnd_mode == ROUNDING_DOWN) || (rnd_mode == ROUNDING_UP) ||
+ (rnd_mode == ROUNDING_TO_ZERO) || // no double rounding error is possible
+ ((res.w[HIGH_128W] & MASK_NAN) == MASK_NAN) || //res=QNaN (cannot be SNaN)
+ ((res.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF)) { // result is infinity
+#if DECIMAL_CALL_BY_REFERENCE
+ bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#else
+ res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#endif
+ // determine the unbiased exponent of the result
+ unbexp = ((res1 >> 53) & 0x3ff) - 398;
+
+ // 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);
+ } 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;
+ }
+
+ *pfpsf |= save_fpsf;
+ BID_RETURN (res1);
+ } // else continue, and use rounding to nearest to round to 16 digits
+
+ // at this point the result is rounded to nearest (even or away) to 34 digits
+ // (or less if exact), and it is zero or finite non-zero canonical [sub]normal
+ sign = res.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
+ exp = res.w[HIGH_128W] & MASK_EXP; // biased and shifted left 49 bits
+ unbexp = (exp >> 49) - 6176;
+ C.w[1] = res.w[HIGH_128W] & MASK_COEFF;
+ C.w[0] = res.w[LOW_128W];
+
+ if ((C.w[1] == 0x0 && C.w[0] == 0x0) || // result is zero
+ (unbexp <= (-398 - 35)) || (unbexp >= (369 + 16))) {
+ // clear under/overflow
+#if DECIMAL_CALL_BY_REFERENCE
+ bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#else
+ res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#endif
+ *pfpsf |= save_fpsf;
+ BID_RETURN (res1);
+ } // else continue
+
+ // -398 - 34 <= unbexp <= 369 + 15
+ if (rnd_mode == ROUNDING_TIES_AWAY) {
+ // apply correction, if needed, to make the result rounded to nearest-even
+ if (is_midpoint_gt_even) {
+ // res = res - 1
+ res1--; // res1 is now even
+ } // else the result is already correctly rounded to nearest-even
+ }
+ // at this point the result is finite, non-zero canonical normal or subnormal,
+ // and in most cases overflow or underflow will not occur
+
+ // determine the number of digits q in the result
+ // q = nr. of decimal digits in x
+ // determine first the nr. of bits in x
+ if (C.w[1] == 0) {
+ if (C.w[0] >= 0x0020000000000000ull) { // x >= 2^53
+ // split the 64-bit value in two 32-bit halves to avoid rounding errors
+ if (C.w[0] >= 0x0000000100000000ull) { // x >= 2^32
+ tmp.d = (double) (C.w[0] >> 32); // exact conversion
+ nr_bits =
+ 33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ } else { // x < 2^32
+ tmp.d = (double) (C.w[0]); // exact conversion
+ nr_bits =
+ 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ }
+ } else { // if x < 2^53
+ tmp.d = (double) C.w[0]; // exact conversion
+ nr_bits =
+ 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ }
+ } else { // C.w[1] != 0 => nr. bits = 64 + nr_bits (C.w[1])
+ tmp.d = (double) C.w[1]; // exact conversion
+ nr_bits =
+ 65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+ }
+ q = nr_digits[nr_bits - 1].digits;
+ if (q == 0) {
+ q = nr_digits[nr_bits - 1].digits1;
+ if (C.w[1] > nr_digits[nr_bits - 1].threshold_hi ||
+ (C.w[1] == nr_digits[nr_bits - 1].threshold_hi &&
+ C.w[0] >= nr_digits[nr_bits - 1].threshold_lo))
+ q++;
+ }
+ // if q > 16, round to nearest even to 16 digits (but for underflow it may
+ // have to be truncated even more)
+ if (q > 16) {
+ x0 = q - 16;
+ if (q <= 18) {
+ round64_2_18 (q, x0, C.w[0], &res1, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ } else { // 19 <= q <= 34
+ round128_19_38 (q, x0, C, &res128, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ res1 = res128.w[0]; // the result fits in 64 bits
+ }
+ unbexp = unbexp + x0;
+ if (incr_exp)
+ unbexp++;
+ q = 16; // need to set in case denormalization is necessary
+ } else {
+ // the result does not require a second rounding (and it must have
+ // been exact in the first rounding, since q <= 16)
+ res1 = C.w[0];
+ }
+
+ // avoid a double rounding error
+ if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
+ is_midpoint_lt_even) { // double rounding error upward
+ // res = res - 1
+ res1--; // res1 becomes odd
+ is_midpoint_lt_even = 0;
+ is_inexact_lt_midpoint = 1;
+ if (res1 == 0x00038d7ea4c67fffull) { // 10^15 - 1
+ res1 = 0x002386f26fc0ffffull; // 10^16 - 1
+ unbexp--;
+ }
+ } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
+ is_midpoint_gt_even) { // double rounding error downward
+ // res = res + 1
+ res1++; // res1 becomes odd (so it cannot be 10^16)
+ is_midpoint_gt_even = 0;
+ is_inexact_gt_midpoint = 1;
+ } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ // if this second rounding was exact the result may still be
+ // inexact because of the first rounding
+ if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+ is_inexact_gt_midpoint = 1;
+ }
+ if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+ is_inexact_lt_midpoint = 1;
+ }
+ } else if (is_midpoint_gt_even &&
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ // pulled up to a midpoint
+ is_inexact_lt_midpoint = 1;
+ is_inexact_gt_midpoint = 0;
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 0;
+ } else if (is_midpoint_lt_even &&
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ // pulled down to a midpoint
+ is_inexact_lt_midpoint = 0;
+ is_inexact_gt_midpoint = 1;
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 0;
+ } else {
+ ;
+ }
+ // this is the result rounded correctly to nearest even, with unbounded exp.
+
+ // check for overflow
+ if (q + unbexp > P16 + expmax16) {
+ res1 = sign | 0x7800000000000000ull;
+ *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+ *pfpsf |= save_fpsf;
+ BID_RETURN (res1)
+ } else if (unbexp > expmax16) { // q + unbexp <= P16 + expmax16
+ // not overflow; the result must be exact, and we can multiply res1 by
+ // 10^(unbexp - expmax16) and the product will fit in 16 decimal digits
+ scale = unbexp - expmax16;
+ res1 = res1 * ten2k64[scale]; // res1 * 10^scale
+ unbexp = expmax16; // unbexp - scale
+ } else {
+ ; // continue
+ }
+
+ // check for underflow
+ if (q + unbexp < P16 + expmin16) {
+ if (unbexp < expmin16) {
+ // we must truncate more of res
+ x0 = expmin16 - unbexp; // x0 >= 1
+ is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+ is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+ is_midpoint_lt_even0 = is_midpoint_lt_even;
+ is_midpoint_gt_even0 = is_midpoint_gt_even;
+ is_inexact_lt_midpoint = 0;
+ is_inexact_gt_midpoint = 0;
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 0;
+ // the number of decimal digits in res1 is q
+ if (x0 < q) { // 1 <= x0 <= q-1 => round res to q - x0 digits
+ // 2 <= q <= 16, 1 <= x0 <= 15
+ round64_2_18 (q, x0, res1, &res1, &incr_exp,
+ &is_midpoint_lt_even, &is_midpoint_gt_even,
+ &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+ if (incr_exp) {
+ // res1 = 10^(q-x0), 1 <= q - x0 <= q - 1, 1 <= q - x0 <= 15
+ res1 = ten2k64[q - x0];
+ }
+ unbexp = unbexp + x0; // expmin16
+ } else if (x0 == q) {
+ // the second rounding is for 0.d(0)d(1)...d(q-1) * 10^emin
+ // determine relationship with 1/2 ulp
+ // q <= 16
+ if (res1 < midpoint64[q - 1]) { // < 1/2 ulp
+ lt_half_ulp = 1;
+ is_inexact_lt_midpoint = 1;
+ } else if (res1 == midpoint64[q - 1]) { // = 1/2 ulp
+ eq_half_ulp = 1;
+ is_midpoint_gt_even = 1;
+ } else { // > 1/2 ulp
+ // gt_half_ulp = 1;
+ is_inexact_gt_midpoint = 1;
+ }
+ if (lt_half_ulp || eq_half_ulp) {
+ // res = +0.0 * 10^expmin16
+ res1 = 0x0000000000000000ull;
+ } else { // if (gt_half_ulp)
+ // res = +1 * 10^expmin16
+ res1 = 0x0000000000000001ull;
+ }
+ unbexp = expmin16;
+ } else { // if (x0 > q)
+ // the second rounding is for 0.0...d(0)d(1)...d(q-1) * 10^emin
+ res1 = 0x0000000000000000ull;
+ unbexp = expmin16;
+ is_inexact_lt_midpoint = 1;
+ }
+ // avoid a double rounding error
+ if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
+ is_midpoint_lt_even) { // double rounding error upward
+ // res = res - 1
+ res1--; // res1 becomes odd
+ is_midpoint_lt_even = 0;
+ is_inexact_lt_midpoint = 1;
+ } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
+ is_midpoint_gt_even) { // double rounding error downward
+ // res = res + 1
+ res1++; // res1 becomes odd
+ is_midpoint_gt_even = 0;
+ is_inexact_gt_midpoint = 1;
+ } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+ !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+ // if this rounding was exact the result may still be
+ // inexact because of the previous roundings
+ if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+ is_inexact_gt_midpoint = 1;
+ }
+ if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+ is_inexact_lt_midpoint = 1;
+ }
+ } else if (is_midpoint_gt_even &&
+ (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+ // pulled up to a midpoint
+ is_inexact_lt_midpoint = 1;
+ is_inexact_gt_midpoint = 0;
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 0;
+ } else if (is_midpoint_lt_even &&
+ (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+ // pulled down to a midpoint
+ is_inexact_lt_midpoint = 0;
+ is_inexact_gt_midpoint = 1;
+ is_midpoint_lt_even = 0;
+ is_midpoint_gt_even = 0;
+ } else {
+ ;
+ }
+ }
+ // else if unbexp >= emin then q < P (because q + unbexp < P16 + expmin16)
+ // and the result is tiny and exact
+
+ // check for inexact result
+ if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+ is_midpoint_lt_even || is_midpoint_gt_even ||
+ 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_midpoint || is_inexact_gt_midpoint ||
+ is_midpoint_lt_even || is_midpoint_gt_even) {
+ *pfpsf |= INEXACT_EXCEPTION;
+ }
+ // this is the result rounded correctly to nearest, with bounded exponent
+
+ if (rnd_mode == ROUNDING_TIES_AWAY && is_midpoint_gt_even) { // correction
+ // res = res + 1
+ res1++; // res1 is now odd
+ } // else the result is already correct
+
+ // assemble the result
+ if (res1 < 0x0020000000000000ull) { // res < 2^53
+ res1 = sign | ((UINT64) (unbexp + 398) << 53) | res1;
+ } else { // res1 >= 2^53
+ res1 = sign | MASK_STEERING_BITS |
+ ((UINT64) (unbexp + 398) << 51) | (res1 & MASK_BINARY_SIG2);
+ }
+ *pfpsf |= save_fpsf;
+ BID_RETURN (res1);
+}
OSDN Git Service
