blob: d2f70a1075da4232faac78e8294171eb8bb76ae0 [file] [log] [blame]
 /* Copyright (C) 2007-2019 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 . */ /***************************************************************************** * BID64 add ***************************************************************************** * * Algorithm description: * * if(exponent_a < exponent_b) * switch a, b * diff_expon = exponent_a - exponent_b * if(diff_expon > 16) * return normalize(a) * if(coefficient_a*10^diff_expon guaranteed below 2^62) * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b * if(|S|<10^16) * return get_BID64(sign(S),exponent_b,|S|) * else * determine number of extra digits in S (1, 2, or 3) * return rounded result * else // large exponent difference * if(number_digits(coefficient_a*10^diff_expon) +/- 10^16) * guaranteed the same as * number_digits(coefficient_a*10^diff_expon) ) * S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-diff_expon)) * corr = 10^16 + (sign_a^sign_b)*coefficient_b * corr*10^exponent_b is rounded so it aligns with S*10^exponent_S * return get_BID64(sign_a,exponent(S),S+rounded(corr)) * else * add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b * in 128-bit integer arithmetic, then round to 16 decimal digits * * ****************************************************************************/ #include "bid_internal.h" #if DECIMAL_CALL_BY_REFERENCE void bid64_add (UINT64 * pres, UINT64 * px, UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM); #else UINT64 bid64_add (UINT64 x, UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM); #endif #if DECIMAL_CALL_BY_REFERENCE void bid64_sub (UINT64 * pres, UINT64 * px, UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { UINT64 y = *py; #if !DECIMAL_GLOBAL_ROUNDING _IDEC_round rnd_mode = *prnd_mode; #endif // check if y is not NaN if (((y & NAN_MASK64) != NAN_MASK64)) y ^= 0x8000000000000000ull; bid64_add (pres, px, &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); } #else UINT64 bid64_sub (UINT64 x, UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { // check if y is not NaN if (((y & NAN_MASK64) != NAN_MASK64)) y ^= 0x8000000000000000ull; return bid64_add (x, y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); } #endif #if DECIMAL_CALL_BY_REFERENCE void bid64_add (UINT64 * pres, UINT64 * px, UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { UINT64 x, y; #else UINT64 bid64_add (UINT64 x, UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { #endif UINT128 CA, CT, CT_new; UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new; UINT64 valid_x, valid_y; UINT64 res; 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; int_double tempx; int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon; int bin_expon_ca, extra_digits, amount, scale_k, scale_ca; unsigned rmode, status; #if DECIMAL_CALL_BY_REFERENCE #if !DECIMAL_GLOBAL_ROUNDING _IDEC_round rnd_mode = *prnd_mode; #endif x = *px; y = *py; #endif valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); // unpack arguments, check for NaN or Infinity if (!valid_x) { // x is Inf. or NaN // test if x is NaN if ((x & NAN_MASK64) == NAN_MASK64) { #ifdef SET_STATUS_FLAGS if (((x & SNAN_MASK64) == SNAN_MASK64) // sNaN || ((y & SNAN_MASK64) == SNAN_MASK64)) __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif res = coefficient_x & QUIET_MASK64; BID_RETURN (res); } // x is Infinity? if ((x & INFINITY_MASK64) == INFINITY_MASK64) { // check if y is Inf if (((y & NAN_MASK64) == INFINITY_MASK64)) { if (sign_x == (y & 0x8000000000000000ull)) { res = coefficient_x; BID_RETURN (res); } // return NaN { #ifdef SET_STATUS_FLAGS __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif res = NAN_MASK64; BID_RETURN (res); } } // check if y is NaN if (((y & NAN_MASK64) == NAN_MASK64)) { res = coefficient_y & QUIET_MASK64; #ifdef SET_STATUS_FLAGS if (((y & SNAN_MASK64) == SNAN_MASK64)) __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif BID_RETURN (res); } // otherwise return +/-Inf { res = coefficient_x; BID_RETURN (res); } } // x is 0 { if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) { if (exponent_y <= exponent_x) { res = y; BID_RETURN (res); } } } } if (!valid_y) { // y is Inf. or NaN? if (((y & INFINITY_MASK64) == INFINITY_MASK64)) { #ifdef SET_STATUS_FLAGS if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif res = coefficient_y & QUIET_MASK64; BID_RETURN (res); } // y is 0 if (!coefficient_x) { // x==0 if (exponent_x <= exponent_y) res = ((UINT64) exponent_x) << 53; else res = ((UINT64) exponent_y) << 53; if (sign_x == sign_y) res |= sign_x; #ifndef IEEE_ROUND_NEAREST_TIES_AWAY #ifndef IEEE_ROUND_NEAREST if (rnd_mode == ROUNDING_DOWN && sign_x != sign_y) res |= 0x8000000000000000ull; #endif #endif BID_RETURN (res); } else if (exponent_y >= exponent_x) { res = x; BID_RETURN (res); } } // 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 --- // version 2 (original) tempx.d = (double) coefficient_a; bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 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 (pfpsf, INEXACT_EXCEPTION); } #endif #ifndef IEEE_ROUND_NEAREST_TIES_AWAY #ifndef IEEE_ROUND_NEAREST if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST { switch (rnd_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--; } res = fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a, rnd_mode, pfpsf); BID_RETURN (res); } } // 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 (rnd_mode == ROUNDING_DOWN && (!coefficient_a) && sign_a != sign_b) sign_s = 0x8000000000000000ull; #endif #endif res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a); BID_RETURN (res); } // 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 = rnd_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 = rnd_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 if (coefficient_a >= power10_table_128[scale_ca].w[0]) { scale_ca++; } scale_k = 16 - scale_ca; // 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 // this test 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 += 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; //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y; 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 (pfpsf, status); #endif res = fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64, rnd_mode, pfpsf); BID_RETURN (res); }