diff options
Diffstat (limited to 'gcc-4.4.3/libgcc/config/libbid/bid_inline_add.h')
| -rw-r--r-- | gcc-4.4.3/libgcc/config/libbid/bid_inline_add.h | 1253 |
1 files changed, 1253 insertions, 0 deletions
diff --git a/gcc-4.4.3/libgcc/config/libbid/bid_inline_add.h b/gcc-4.4.3/libgcc/config/libbid/bid_inline_add.h new file mode 100644 index 000000000..76246fab6 --- /dev/null +++ b/gcc-4.4.3/libgcc/config/libbid/bid_inline_add.h @@ -0,0 +1,1253 @@ +/* Copyright (C) 2007, 2009 Free Software Foundation, Inc. + +This file is part of GCC. + +GCC is free software; you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free +Software Foundation; either version 3, or (at your option) any later +version. + +GCC is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +<http://www.gnu.org/licenses/>. */ + +/***************************************************************************** + * + * Helper add functions (for fma) + * + * __BID_INLINE__ UINT64 get_add64( + * UINT64 sign_x, int exponent_x, UINT64 coefficient_x, + * UINT64 sign_y, int exponent_y, UINT64 coefficient_y, + * int rounding_mode) + * + * __BID_INLINE__ UINT64 get_add128( + * UINT64 sign_x, int exponent_x, UINT64 coefficient_x, + * UINT64 sign_y, int final_exponent_y, UINT128 CY, + * int extra_digits, int rounding_mode) + * + ***************************************************************************** + * + * Algorithm description: + * + * get_add64: same as BID64 add, but arguments are unpacked and there + * are no special case checks + * + * get_add128: add 64-bit coefficient to 128-bit product (which contains + * 16+extra_digits decimal digits), + * return BID64 result + * - the exponents are compared and the two coefficients are + * properly aligned for addition/subtraction + * - multiple paths are needed + * - final result exponent is calculated and the lower term is + * rounded first if necessary, to avoid manipulating + * coefficients longer than 128 bits + * + ****************************************************************************/ + +#ifndef _INLINE_BID_ADD_H_ +#define _INLINE_BID_ADD_H_ + +#include "bid_internal.h" + +#define MAX_FORMAT_DIGITS 16 +#define DECIMAL_EXPONENT_BIAS 398 +#define MASK_BINARY_EXPONENT 0x7ff0000000000000ull +#define BINARY_EXPONENT_BIAS 0x3ff +#define UPPER_EXPON_LIMIT 51 + +/////////////////////////////////////////////////////////////////////// +// +// get_add64() is essentially the same as bid_add(), except that +// the arguments are unpacked +// +////////////////////////////////////////////////////////////////////// +__BID_INLINE__ UINT64 +get_add64 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x, + UINT64 sign_y, int exponent_y, UINT64 coefficient_y, + int rounding_mode, unsigned *fpsc) { + UINT128 CA, CT, CT_new; + UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab, + rem_a; + UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp, + C64_new; + int_double tempx; + int exponent_a, exponent_b, diff_dec_expon; + int bin_expon_ca, extra_digits, amount, scale_k, scale_ca; + unsigned rmode, status; + + // sort arguments by exponent + if (exponent_x <= exponent_y) { + sign_a = sign_y; + exponent_a = exponent_y; + coefficient_a = coefficient_y; + sign_b = sign_x; + exponent_b = exponent_x; + coefficient_b = coefficient_x; + } else { + sign_a = sign_x; + exponent_a = exponent_x; + coefficient_a = coefficient_x; + sign_b = sign_y; + exponent_b = exponent_y; + coefficient_b = coefficient_y; + } + + // exponent difference + diff_dec_expon = exponent_a - exponent_b; + + /* get binary coefficients of x and y */ + + //--- get number of bits in the coefficients of x and y --- + + tempx.d = (double) coefficient_a; + bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + + if (!coefficient_a) { + return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode, + fpsc); + } + if (diff_dec_expon > MAX_FORMAT_DIGITS) { + // normalize a to a 16-digit coefficient + + scale_ca = estimate_decimal_digits[bin_expon_ca]; + if (coefficient_a >= power10_table_128[scale_ca].w[0]) + scale_ca++; + + scale_k = 16 - scale_ca; + + coefficient_a *= power10_table_128[scale_k].w[0]; + + diff_dec_expon -= scale_k; + exponent_a -= scale_k; + + /* get binary coefficients of x and y */ + + //--- get number of bits in the coefficients of x and y --- + tempx.d = (double) coefficient_a; + bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + + if (diff_dec_expon > MAX_FORMAT_DIGITS) { +#ifdef SET_STATUS_FLAGS + if (coefficient_b) { + __set_status_flags (fpsc, INEXACT_EXCEPTION); + } +#endif + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (((rounding_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST + { + switch (rounding_mode) { + case ROUNDING_DOWN: + if (sign_b) { + coefficient_a -= ((((SINT64) sign_a) >> 63) | 1); + if (coefficient_a < 1000000000000000ull) { + exponent_a--; + coefficient_a = 9999999999999999ull; + } else if (coefficient_a >= 10000000000000000ull) { + exponent_a++; + coefficient_a = 1000000000000000ull; + } + } + break; + case ROUNDING_UP: + if (!sign_b) { + coefficient_a += ((((SINT64) sign_a) >> 63) | 1); + if (coefficient_a < 1000000000000000ull) { + exponent_a--; + coefficient_a = 9999999999999999ull; + } else if (coefficient_a >= 10000000000000000ull) { + exponent_a++; + coefficient_a = 1000000000000000ull; + } + } + break; + default: // RZ + if (sign_a != sign_b) { + coefficient_a--; + if (coefficient_a < 1000000000000000ull) { + exponent_a--; + coefficient_a = 9999999999999999ull; + } + } + break; + } + } else +#endif +#endif + // check special case here + if ((coefficient_a == 1000000000000000ull) + && (diff_dec_expon == MAX_FORMAT_DIGITS + 1) + && (sign_a ^ sign_b) + && (coefficient_b > 5000000000000000ull)) { + coefficient_a = 9999999999999999ull; + exponent_a--; + } + + return get_BID64 (sign_a, exponent_a, coefficient_a, + rounding_mode, fpsc); + } + } + // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62 + if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) { + // coefficient_a*10^(exponent_a-exponent_b)<2^63 + + // multiply by 10^(exponent_a-exponent_b) + coefficient_a *= power10_table_128[diff_dec_expon].w[0]; + + // sign mask + sign_b = ((SINT64) sign_b) >> 63; + // apply sign to coeff. of b + coefficient_b = (coefficient_b + sign_b) ^ sign_b; + + // apply sign to coefficient a + sign_a = ((SINT64) sign_a) >> 63; + coefficient_a = (coefficient_a + sign_a) ^ sign_a; + + coefficient_a += coefficient_b; + // get sign + sign_s = ((SINT64) coefficient_a) >> 63; + coefficient_a = (coefficient_a + sign_s) ^ sign_s; + sign_s &= 0x8000000000000000ull; + + // coefficient_a < 10^16 ? + if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) { +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (rounding_mode == ROUNDING_DOWN && (!coefficient_a) + && sign_a != sign_b) + sign_s = 0x8000000000000000ull; +#endif +#endif + return get_BID64 (sign_s, exponent_b, coefficient_a, + rounding_mode, fpsc); + } + // otherwise rounding is necessary + + // already know coefficient_a<10^19 + // coefficient_a < 10^17 ? + if (coefficient_a < power10_table_128[17].w[0]) + extra_digits = 1; + else if (coefficient_a < power10_table_128[18].w[0]) + extra_digits = 2; + else + extra_digits = 3; + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + rmode = rounding_mode; + if (sign_s && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; +#else + rmode = 0; +#endif +#else + rmode = 0; +#endif + coefficient_a += round_const_table[rmode][extra_digits]; + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_64x64_to_128 (CT, coefficient_a, + reciprocals10_64[extra_digits]); + + // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 + amount = short_recip_scale[extra_digits]; + C64 = CT.w[1] >> amount; + + } else { + // coefficient_a*10^(exponent_a-exponent_b) is large + sign_s = sign_a; + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + rmode = rounding_mode; + if (sign_s && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; +#else + rmode = 0; +#endif +#else + rmode = 0; +#endif + + // check whether we can take faster path + scale_ca = estimate_decimal_digits[bin_expon_ca]; + + sign_ab = sign_a ^ sign_b; + sign_ab = ((SINT64) sign_ab) >> 63; + + // T1 = 10^(16-diff_dec_expon) + T1 = power10_table_128[16 - diff_dec_expon].w[0]; + + // get number of digits in coefficient_a + //P_ca = power10_table_128[scale_ca].w[0]; + //P_ca_m1 = power10_table_128[scale_ca-1].w[0]; + if (coefficient_a >= power10_table_128[scale_ca].w[0]) { + scale_ca++; + //P_ca_m1 = P_ca; + //P_ca = power10_table_128[scale_ca].w[0]; + } + + scale_k = 16 - scale_ca; + + // apply sign + //Ts = (T1 + sign_ab) ^ sign_ab; + + // test range of ca + //X = coefficient_a + Ts - P_ca_m1; + + // addition + saved_ca = coefficient_a - T1; + coefficient_a = + (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0]; + extra_digits = diff_dec_expon - scale_k; + + // apply sign + saved_cb = (coefficient_b + sign_ab) ^ sign_ab; + // add 10^16 and rounding constant + coefficient_b = + saved_cb + 10000000000000000ull + + round_const_table[rmode][extra_digits]; + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_64x64_to_128 (CT, coefficient_b, + reciprocals10_64[extra_digits]); + + // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 + amount = short_recip_scale[extra_digits]; + C0_64 = CT.w[1] >> amount; + + // result coefficient + C64 = C0_64 + coefficient_a; + // filter out difficult (corner) cases + // the following test is equivalent to + // ( (initial_coefficient_a + Ts) < P_ca && + // (initial_coefficient_a + Ts) > P_ca_m1 ), + // which ensures the number of digits in coefficient_a does not change + // after adding (the appropriately scaled and rounded) coefficient_b + if ((UINT64) (C64 - 1000000000000000ull - 1) > + 9000000000000000ull - 2) { + if (C64 >= 10000000000000000ull) { + // result has more than 16 digits + if (!scale_k) { + // must divide coeff_a by 10 + saved_ca = saved_ca + T1; + __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull); + //reciprocals10_64[1]); + coefficient_a = CA.w[1] >> 1; + rem_a = + saved_ca - (coefficient_a << 3) - (coefficient_a << 1); + coefficient_a = coefficient_a - T1; + + saved_cb += + /*90000000000000000 */ +rem_a * + power10_table_128[diff_dec_expon].w[0]; + } else + coefficient_a = + (SINT64) (saved_ca - T1 - + (T1 << 3)) * (SINT64) power10_table_128[scale_k - + 1].w[0]; + + extra_digits++; + coefficient_b = + saved_cb + 100000000000000000ull + + round_const_table[rmode][extra_digits]; + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_64x64_to_128 (CT, coefficient_b, + reciprocals10_64[extra_digits]); + + // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 + amount = short_recip_scale[extra_digits]; + C0_64 = CT.w[1] >> amount; + + // result coefficient + C64 = C0_64 + coefficient_a; + } else if (C64 <= 1000000000000000ull) { + // less than 16 digits in result + coefficient_a = + (SINT64) saved_ca *(SINT64) power10_table_128[scale_k + + 1].w[0]; + //extra_digits --; + exponent_b--; + coefficient_b = + (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull + + round_const_table[rmode][extra_digits]; + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_64x64_to_128 (CT_new, coefficient_b, + reciprocals10_64[extra_digits]); + + // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 + amount = short_recip_scale[extra_digits]; + C0_64 = CT_new.w[1] >> amount; + + // result coefficient + C64_new = C0_64 + coefficient_a; + if (C64_new < 10000000000000000ull) { + C64 = C64_new; +#ifdef SET_STATUS_FLAGS + CT = CT_new; +#endif + } else + exponent_b++; + } + + } + + } + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (rmode == 0) //ROUNDING_TO_NEAREST +#endif + if (C64 & 1) { + // check whether fractional part of initial_P/10^extra_digits + // is exactly .5 + // this is the same as fractional part of + // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero + + // get remainder + remainder_h = CT.w[1] << (64 - amount); + + // test whether fractional part is 0 + if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) { + C64--; + } + } +#endif + +#ifdef SET_STATUS_FLAGS + status = INEXACT_EXCEPTION; + + // get remainder + remainder_h = CT.w[1] << (64 - amount); + + switch (rmode) { + case ROUNDING_TO_NEAREST: + case ROUNDING_TIES_AWAY: + // test whether fractional part is 0 + if ((remainder_h == 0x8000000000000000ull) + && (CT.w[0] < reciprocals10_64[extra_digits])) + status = EXACT_STATUS; + break; + case ROUNDING_DOWN: + case ROUNDING_TO_ZERO: + if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) + status = EXACT_STATUS; + break; + default: + // round up + __add_carry_out (tmp, carry, CT.w[0], + reciprocals10_64[extra_digits]); + if ((remainder_h >> (64 - amount)) + carry >= + (((UINT64) 1) << amount)) + status = EXACT_STATUS; + break; + } + __set_status_flags (fpsc, status); + +#endif + + return get_BID64 (sign_s, exponent_b + extra_digits, C64, + rounding_mode, fpsc); +} + + +/////////////////////////////////////////////////////////////////// +// round 128-bit coefficient and return result in BID64 format +// do not worry about midpoint cases +////////////////////////////////////////////////////////////////// +static UINT64 +__bid_simple_round64_sticky (UINT64 sign, int exponent, UINT128 P, + int extra_digits, int rounding_mode, + unsigned *fpsc) { + UINT128 Q_high, Q_low, C128; + UINT64 C64; + int amount, rmode; + + rmode = rounding_mode; +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (sign && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; +#endif +#endif + __add_128_64 (P, P, round_const_table[rmode][extra_digits]); + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_128x128_full (Q_high, Q_low, P, + reciprocals10_128[extra_digits]); + + // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 + amount = recip_scale[extra_digits]; + __shr_128 (C128, Q_high, amount); + + C64 = __low_64 (C128); + +#ifdef SET_STATUS_FLAGS + + __set_status_flags (fpsc, INEXACT_EXCEPTION); + +#endif + + return get_BID64 (sign, exponent, C64, rounding_mode, fpsc); +} + +/////////////////////////////////////////////////////////////////// +// round 128-bit coefficient and return result in BID64 format +/////////////////////////////////////////////////////////////////// +static UINT64 +__bid_full_round64 (UINT64 sign, int exponent, UINT128 P, + int extra_digits, int rounding_mode, + unsigned *fpsc) { + UINT128 Q_high, Q_low, C128, Stemp, PU; + UINT64 remainder_h, C64, carry, CY; + int amount, amount2, rmode, status = 0; + + if (exponent < 0) { + if (exponent >= -16 && (extra_digits + exponent < 0)) { + extra_digits = -exponent; +#ifdef SET_STATUS_FLAGS + if (extra_digits > 0) { + rmode = rounding_mode; + if (sign && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; + __add_128_128 (PU, P, + round_const_table_128[rmode][extra_digits]); + if (__unsigned_compare_gt_128 + (power10_table_128[extra_digits + 15], PU)) + status = UNDERFLOW_EXCEPTION; + } +#endif + } + } + + if (extra_digits > 0) { + exponent += extra_digits; + rmode = rounding_mode; + if (sign && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; + __add_128_128 (P, P, round_const_table_128[rmode][extra_digits]); + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_128x128_full (Q_high, Q_low, P, + reciprocals10_128[extra_digits]); + + // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 + amount = recip_scale[extra_digits]; + __shr_128_long (C128, Q_high, amount); + + C64 = __low_64 (C128); + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (rmode == 0) //ROUNDING_TO_NEAREST +#endif + if (C64 & 1) { + // check whether fractional part of initial_P/10^extra_digits + // is exactly .5 + + // get remainder + amount2 = 64 - amount; + remainder_h = 0; + remainder_h--; + remainder_h >>= amount2; + remainder_h = remainder_h & Q_high.w[0]; + + if (!remainder_h + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) { + C64--; + } + } +#endif + +#ifdef SET_STATUS_FLAGS + status |= INEXACT_EXCEPTION; + + // get remainder + remainder_h = Q_high.w[0] << (64 - amount); + + switch (rmode) { + case ROUNDING_TO_NEAREST: + case ROUNDING_TIES_AWAY: + // test whether fractional part is 0 + if (remainder_h == 0x8000000000000000ull + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) + status = EXACT_STATUS; + break; + case ROUNDING_DOWN: + case ROUNDING_TO_ZERO: + if (!remainder_h + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) + status = EXACT_STATUS; + break; + default: + // round up + __add_carry_out (Stemp.w[0], CY, Q_low.w[0], + reciprocals10_128[extra_digits].w[0]); + __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], + reciprocals10_128[extra_digits].w[1], CY); + if ((remainder_h >> (64 - amount)) + carry >= + (((UINT64) 1) << amount)) + status = EXACT_STATUS; + } + + __set_status_flags (fpsc, status); + +#endif + } else { + C64 = P.w[0]; + if (!C64) { + sign = 0; + if (rounding_mode == ROUNDING_DOWN) + sign = 0x8000000000000000ull; + } + } + return get_BID64 (sign, exponent, C64, rounding_mode, fpsc); +} + +///////////////////////////////////////////////////////////////////////////////// +// round 192-bit coefficient (P, remainder_P) and return result in BID64 format +// the lowest 64 bits (remainder_P) are used for midpoint checking only +//////////////////////////////////////////////////////////////////////////////// +static UINT64 +__bid_full_round64_remainder (UINT64 sign, int exponent, UINT128 P, + int extra_digits, UINT64 remainder_P, + int rounding_mode, unsigned *fpsc, + unsigned uf_status) { + UINT128 Q_high, Q_low, C128, Stemp; + UINT64 remainder_h, C64, carry, CY; + int amount, amount2, rmode, status = uf_status; + + rmode = rounding_mode; + if (sign && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; + if (rmode == ROUNDING_UP && remainder_P) { + P.w[0]++; + if (!P.w[0]) + P.w[1]++; + } + + if (extra_digits) { + __add_128_64 (P, P, round_const_table[rmode][extra_digits]); + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_128x128_full (Q_high, Q_low, P, + reciprocals10_128[extra_digits]); + + // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 + amount = recip_scale[extra_digits]; + __shr_128 (C128, Q_high, amount); + + C64 = __low_64 (C128); + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (rmode == 0) //ROUNDING_TO_NEAREST +#endif + if (!remainder_P && (C64 & 1)) { + // check whether fractional part of initial_P/10^extra_digits + // is exactly .5 + + // get remainder + amount2 = 64 - amount; + remainder_h = 0; + remainder_h--; + remainder_h >>= amount2; + remainder_h = remainder_h & Q_high.w[0]; + + if (!remainder_h + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) { + C64--; + } + } +#endif + +#ifdef SET_STATUS_FLAGS + status |= INEXACT_EXCEPTION; + + if (!remainder_P) { + // get remainder + remainder_h = Q_high.w[0] << (64 - amount); + + switch (rmode) { + case ROUNDING_TO_NEAREST: + case ROUNDING_TIES_AWAY: + // test whether fractional part is 0 + if (remainder_h == 0x8000000000000000ull + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) + status = EXACT_STATUS; + break; + case ROUNDING_DOWN: + case ROUNDING_TO_ZERO: + if (!remainder_h + && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] + || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] + && Q_low.w[0] < + reciprocals10_128[extra_digits].w[0]))) + status = EXACT_STATUS; + break; + default: + // round up + __add_carry_out (Stemp.w[0], CY, Q_low.w[0], + reciprocals10_128[extra_digits].w[0]); + __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], + reciprocals10_128[extra_digits].w[1], CY); + if ((remainder_h >> (64 - amount)) + carry >= + (((UINT64) 1) << amount)) + status = EXACT_STATUS; + } + } + __set_status_flags (fpsc, status); + +#endif + } else { + C64 = P.w[0]; +#ifdef SET_STATUS_FLAGS + if (remainder_P) { + __set_status_flags (fpsc, uf_status | INEXACT_EXCEPTION); + } +#endif + } + + return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode, + fpsc); +} + + +/////////////////////////////////////////////////////////////////// +// get P/10^extra_digits +// result fits in 64 bits +/////////////////////////////////////////////////////////////////// +__BID_INLINE__ UINT64 +__truncate (UINT128 P, int extra_digits) +// extra_digits <= 16 +{ + UINT128 Q_high, Q_low, C128; + UINT64 C64; + int amount; + + // get P*(2^M[extra_digits])/10^extra_digits + __mul_128x128_full (Q_high, Q_low, P, + reciprocals10_128[extra_digits]); + + // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 + amount = recip_scale[extra_digits]; + __shr_128 (C128, Q_high, amount); + + C64 = __low_64 (C128); + + return C64; +} + + +/////////////////////////////////////////////////////////////////// +// return number of decimal digits in 128-bit value X +/////////////////////////////////////////////////////////////////// +__BID_INLINE__ int +__get_dec_digits64 (UINT128 X) { + int_double tempx; + int digits_x, bin_expon_cx; + + if (!X.w[1]) { + //--- get number of bits in the coefficients of x and y --- + tempx.d = (double) X.w[0]; + bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + // get number of decimal digits in the coeff_x + digits_x = estimate_decimal_digits[bin_expon_cx]; + if (X.w[0] >= power10_table_128[digits_x].w[0]) + digits_x++; + return digits_x; + } + tempx.d = (double) X.w[1]; + bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + // get number of decimal digits in the coeff_x + digits_x = estimate_decimal_digits[bin_expon_cx + 64]; + if (__unsigned_compare_ge_128 (X, power10_table_128[digits_x])) + digits_x++; + + return digits_x; +} + + +//////////////////////////////////////////////////////////////////////////////// +// +// add 64-bit coefficient to 128-bit coefficient, return result in BID64 format +// +//////////////////////////////////////////////////////////////////////////////// +__BID_INLINE__ UINT64 +get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x, + UINT64 sign_y, int final_exponent_y, UINT128 CY, + int extra_digits, int rounding_mode, unsigned *fpsc) { + UINT128 CY_L, CX, FS, F, CT, ST, T2; + UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y; + SINT64 D = 0; + int_double tempx; + int diff_dec_expon, extra_digits2, exponent_y, status; + int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode; + + // CY has more than 16 decimal digits + + exponent_y = final_exponent_y - extra_digits; + +#ifdef IEEE_ROUND_NEAREST_TIES_AWAY + rounding_mode = 0; +#endif +#ifdef IEEE_ROUND_NEAREST + rounding_mode = 0; +#endif + + if (exponent_x > exponent_y) { + // normalize x + //--- get number of bits in the coefficients of x and y --- + tempx.d = (double) coefficient_x; + bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + // get number of decimal digits in the coeff_x + digits_x = estimate_decimal_digits[bin_expon_cx]; + if (coefficient_x >= power10_table_128[digits_x].w[0]) + digits_x++; + + extra_dx = 16 - digits_x; + coefficient_x *= power10_table_128[extra_dx].w[0]; + if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) { + extra_dx++; + coefficient_x = 10000000000000000ull; + } + exponent_x -= extra_dx; + + if (exponent_x > exponent_y) { + + // exponent_x > exponent_y + diff_dec_expon = exponent_x - exponent_y; + + if (exponent_x <= final_exponent_y + 1) { + __mul_64x64_to_128 (CX, coefficient_x, + power10_table_128[diff_dec_expon].w[0]); + + if (sign_x == sign_y) { + __add_128_128 (CT, CY, CX); + if ((exponent_x > + final_exponent_y) /*&& (final_exponent_y>0) */ ) + extra_digits++; + if (__unsigned_compare_ge_128 + (CT, power10_table_128[16 + extra_digits])) + extra_digits++; + } else { + __sub_128_128 (CT, CY, CX); + if (((SINT64) CT.w[1]) < 0) { + CT.w[0] = 0 - CT.w[0]; + CT.w[1] = 0 - CT.w[1]; + if (CT.w[0]) + CT.w[1]--; + sign_y = sign_x; + } else if (!(CT.w[1] | CT.w[0])) { + sign_y = + (rounding_mode != + ROUNDING_DOWN) ? 0 : 0x8000000000000000ull; + } + if ((exponent_x + 1 >= + final_exponent_y) /*&& (final_exponent_y>=0) */ ) { + extra_digits = __get_dec_digits64 (CT) - 16; + if (extra_digits <= 0) { + if (!CT.w[0] && rounding_mode == ROUNDING_DOWN) + sign_y = 0x8000000000000000ull; + return get_BID64 (sign_y, exponent_y, CT.w[0], + rounding_mode, fpsc); + } + } else + if (__unsigned_compare_gt_128 + (power10_table_128[15 + extra_digits], CT)) + extra_digits--; + } + + return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits, + rounding_mode, fpsc); + } + // diff_dec2+extra_digits is the number of digits to eliminate from + // argument CY + diff_dec2 = exponent_x - final_exponent_y; + + if (diff_dec2 >= 17) { +#ifndef IEEE_ROUND_NEAREST +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY + if ((rounding_mode) & 3) { + switch (rounding_mode) { + case ROUNDING_UP: + if (!sign_y) { + D = ((SINT64) (sign_x ^ sign_y)) >> 63; + D = D + D + 1; + coefficient_x += D; + } + break; + case ROUNDING_DOWN: + if (sign_y) { + D = ((SINT64) (sign_x ^ sign_y)) >> 63; + D = D + D + 1; + coefficient_x += D; + } + break; + case ROUNDING_TO_ZERO: + if (sign_y != sign_x) { + D = 0 - 1; + coefficient_x += D; + } + break; + } + if (coefficient_x < 1000000000000000ull) { + coefficient_x -= D; + coefficient_x = + D + (coefficient_x << 1) + (coefficient_x << 3); + exponent_x--; + } + } +#endif +#endif +#ifdef SET_STATUS_FLAGS + if (CY.w[1] | CY.w[0]) + __set_status_flags (fpsc, INEXACT_EXCEPTION); +#endif + return get_BID64 (sign_x, exponent_x, coefficient_x, + rounding_mode, fpsc); + } + // here exponent_x <= 16+final_exponent_y + + // truncate CY to 16 dec. digits + CYh = __truncate (CY, extra_digits); + + // get remainder + T = power10_table_128[extra_digits].w[0]; + __mul_64x64_to_64 (CY0L, CYh, T); + + remainder_y = CY.w[0] - CY0L; + + // align coeff_x, CYh + __mul_64x64_to_128 (CX, coefficient_x, + power10_table_128[diff_dec2].w[0]); + + if (sign_x == sign_y) { + __add_128_64 (CT, CX, CYh); + if (__unsigned_compare_ge_128 + (CT, power10_table_128[16 + diff_dec2])) + diff_dec2++; + } else { + if (remainder_y) + CYh++; + __sub_128_64 (CT, CX, CYh); + if (__unsigned_compare_gt_128 + (power10_table_128[15 + diff_dec2], CT)) + diff_dec2--; + } + + return __bid_full_round64_remainder (sign_x, final_exponent_y, CT, + diff_dec2, remainder_y, + rounding_mode, fpsc, 0); + } + } + // Here (exponent_x <= exponent_y) + { + diff_dec_expon = exponent_y - exponent_x; + + if (diff_dec_expon > MAX_FORMAT_DIGITS) { + rmode = rounding_mode; + + if ((sign_x ^ sign_y)) { + if (!CY.w[0]) + CY.w[1]--; + CY.w[0]--; + if (__unsigned_compare_gt_128 + (power10_table_128[15 + extra_digits], CY)) { + if (rmode & 3) { + extra_digits--; + final_exponent_y--; + } else { + CY.w[0] = 1000000000000000ull; + CY.w[1] = 0; + extra_digits = 0; + } + } + } + __scale128_10 (CY, CY); + extra_digits++; + CY.w[0] |= 1; + + return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY, + extra_digits, rmode, fpsc); + } + // apply sign to coeff_x + sign_x ^= sign_y; + sign_x = ((SINT64) sign_x) >> 63; + CX.w[0] = (coefficient_x + sign_x) ^ sign_x; + CX.w[1] = sign_x; + + // check whether CY (rounded to 16 digits) and CX have + // any digits in the same position + diff_dec2 = final_exponent_y - exponent_x; + + if (diff_dec2 <= 17) { + // align CY to 10^ex + S = power10_table_128[diff_dec_expon].w[0]; + __mul_64x128_short (CY_L, S, CY); + + __add_128_128 (ST, CY_L, CX); + extra_digits2 = __get_dec_digits64 (ST) - 16; + return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2, + rounding_mode, fpsc); + } + // truncate CY to 16 dec. digits + CYh = __truncate (CY, extra_digits); + + // get remainder + T = power10_table_128[extra_digits].w[0]; + __mul_64x64_to_64 (CY0L, CYh, T); + + coefficient_y = CY.w[0] - CY0L; + // add rounding constant + rmode = rounding_mode; + if (sign_y && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (!(rmode & 3)) //ROUNDING_TO_NEAREST +#endif +#endif + { + coefficient_y += round_const_table[rmode][extra_digits]; + } + // align coefficient_y, coefficient_x + S = power10_table_128[diff_dec_expon].w[0]; + __mul_64x64_to_128 (F, coefficient_y, S); + + // fraction + __add_128_128 (FS, F, CX); + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (rmode == 0) //ROUNDING_TO_NEAREST +#endif + { + // rounding code, here RN_EVEN + // 10^(extra_digits+diff_dec_expon) + T2 = power10_table_128[diff_dec_expon + extra_digits]; + if (__unsigned_compare_gt_128 (FS, T2) + || ((CYh & 1) && __test_equal_128 (FS, T2))) { + CYh++; + __sub_128_128 (FS, FS, T2); + } + } +#endif +#ifndef IEEE_ROUND_NEAREST +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY + if (rmode == 4) //ROUNDING_TO_NEAREST +#endif + { + // rounding code, here RN_AWAY + // 10^(extra_digits+diff_dec_expon) + T2 = power10_table_128[diff_dec_expon + extra_digits]; + if (__unsigned_compare_ge_128 (FS, T2)) { + CYh++; + __sub_128_128 (FS, FS, T2); + } + } +#endif +#ifndef IEEE_ROUND_NEAREST +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY + switch (rmode) { + case ROUNDING_DOWN: + case ROUNDING_TO_ZERO: + if ((SINT64) FS.w[1] < 0) { + CYh--; + if (CYh < 1000000000000000ull) { + CYh = 9999999999999999ull; + final_exponent_y--; + } + } else { + T2 = power10_table_128[diff_dec_expon + extra_digits]; + if (__unsigned_compare_ge_128 (FS, T2)) { + CYh++; + __sub_128_128 (FS, FS, T2); + } + } + break; + case ROUNDING_UP: + if ((SINT64) FS.w[1] < 0) + break; + T2 = power10_table_128[diff_dec_expon + extra_digits]; + if (__unsigned_compare_gt_128 (FS, T2)) { + CYh += 2; + __sub_128_128 (FS, FS, T2); + } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) { + CYh++; + FS.w[1] = FS.w[0] = 0; + } else if (FS.w[1] | FS.w[0]) + CYh++; + break; + } +#endif +#endif + +#ifdef SET_STATUS_FLAGS + status = INEXACT_EXCEPTION; +#ifndef IEEE_ROUND_NEAREST +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY + if (!(rmode & 3)) +#endif +#endif + { + // RN modes + if ((FS.w[1] == + round_const_table_128[0][diff_dec_expon + extra_digits].w[1]) + && (FS.w[0] == + round_const_table_128[0][diff_dec_expon + + extra_digits].w[0])) + status = EXACT_STATUS; + } +#ifndef IEEE_ROUND_NEAREST +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY + else if (!FS.w[1] && !FS.w[0]) + status = EXACT_STATUS; +#endif +#endif + + __set_status_flags (fpsc, status); +#endif + + return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode, + fpsc); + } + +} + +////////////////////////////////////////////////////////////////////////// +// +// If coefficient_z is less than 16 digits long, normalize to 16 digits +// +///////////////////////////////////////////////////////////////////////// +static UINT64 +BID_normalize (UINT64 sign_z, int exponent_z, + UINT64 coefficient_z, UINT64 round_dir, int round_flag, + int rounding_mode, unsigned *fpsc) { + SINT64 D; + int_double tempx; + int digits_z, bin_expon, scale, rmode; + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + rmode = rounding_mode; + if (sign_z && (unsigned) (rmode - 1) < 2) + rmode = 3 - rmode; +#else + if (coefficient_z >= power10_table_128[15].w[0]) + return z; +#endif +#endif + + //--- get number of bits in the coefficients of x and y --- + tempx.d = (double) coefficient_z; + bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + // get number of decimal digits in the coeff_x + digits_z = estimate_decimal_digits[bin_expon]; + if (coefficient_z >= power10_table_128[digits_z].w[0]) + digits_z++; + + scale = 16 - digits_z; + exponent_z -= scale; + if (exponent_z < 0) { + scale += exponent_z; + exponent_z = 0; + } + coefficient_z *= power10_table_128[scale].w[0]; + +#ifdef SET_STATUS_FLAGS + if (round_flag) { + __set_status_flags (fpsc, INEXACT_EXCEPTION); + if (coefficient_z < 1000000000000000ull) + __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); + else if ((coefficient_z == 1000000000000000ull) && !exponent_z + && ((SINT64) (round_dir ^ sign_z) < 0) && round_flag + && (rmode == ROUNDING_DOWN || rmode == ROUNDING_TO_ZERO)) + __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); + } +#endif + +#ifndef IEEE_ROUND_NEAREST_TIES_AWAY +#ifndef IEEE_ROUND_NEAREST + if (round_flag && (rmode & 3)) { + D = round_dir ^ sign_z; + + if (rmode == ROUNDING_UP) { + if (D >= 0) + coefficient_z++; + } else { + if (D < 0) + coefficient_z--; + if (coefficient_z < 1000000000000000ull && exponent_z) { + coefficient_z = 9999999999999999ull; + exponent_z--; + } + } + } +#endif +#endif + + return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode, + fpsc); +} + + +////////////////////////////////////////////////////////////////////////// +// +// 0*10^ey + cz*10^ez, ey<ez +// +////////////////////////////////////////////////////////////////////////// + +__BID_INLINE__ UINT64 +add_zero64 (int exponent_y, UINT64 sign_z, int exponent_z, + UINT64 coefficient_z, unsigned *prounding_mode, + unsigned *fpsc) { + int_double tempx; + int bin_expon, scale_k, scale_cz; + int diff_expon; + + diff_expon = exponent_z - exponent_y; + + tempx.d = (double) coefficient_z; + bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; + scale_cz = estimate_decimal_digits[bin_expon]; + if (coefficient_z >= power10_table_128[scale_cz].w[0]) + scale_cz++; + + scale_k = 16 - scale_cz; + if (diff_expon < scale_k) + scale_k = diff_expon; + coefficient_z *= power10_table_128[scale_k].w[0]; + + return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z, + *prounding_mode, fpsc); +} +#endif |