1*ec02198aSmrg /* Copyright (C) 2007-2020 Free Software Foundation, Inc.
263d1a8abSmrg
363d1a8abSmrg This file is part of GCC.
463d1a8abSmrg
563d1a8abSmrg GCC is free software; you can redistribute it and/or modify it under
663d1a8abSmrg the terms of the GNU General Public License as published by the Free
763d1a8abSmrg Software Foundation; either version 3, or (at your option) any later
863d1a8abSmrg version.
963d1a8abSmrg
1063d1a8abSmrg GCC is distributed in the hope that it will be useful, but WITHOUT ANY
1163d1a8abSmrg WARRANTY; without even the implied warranty of MERCHANTABILITY or
1263d1a8abSmrg FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
1363d1a8abSmrg for more details.
1463d1a8abSmrg
1563d1a8abSmrg Under Section 7 of GPL version 3, you are granted additional
1663d1a8abSmrg permissions described in the GCC Runtime Library Exception, version
1763d1a8abSmrg 3.1, as published by the Free Software Foundation.
1863d1a8abSmrg
1963d1a8abSmrg You should have received a copy of the GNU General Public License and
2063d1a8abSmrg a copy of the GCC Runtime Library Exception along with this program;
2163d1a8abSmrg see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
2263d1a8abSmrg <http://www.gnu.org/licenses/>. */
2363d1a8abSmrg
2463d1a8abSmrg /*****************************************************************************
2563d1a8abSmrg * BID64 add
2663d1a8abSmrg *****************************************************************************
2763d1a8abSmrg *
2863d1a8abSmrg * Algorithm description:
2963d1a8abSmrg *
3063d1a8abSmrg * if(exponent_a < exponent_b)
3163d1a8abSmrg * switch a, b
3263d1a8abSmrg * diff_expon = exponent_a - exponent_b
3363d1a8abSmrg * if(diff_expon > 16)
3463d1a8abSmrg * return normalize(a)
3563d1a8abSmrg * if(coefficient_a*10^diff_expon guaranteed below 2^62)
3663d1a8abSmrg * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
3763d1a8abSmrg * if(|S|<10^16)
3863d1a8abSmrg * return get_BID64(sign(S),exponent_b,|S|)
3963d1a8abSmrg * else
4063d1a8abSmrg * determine number of extra digits in S (1, 2, or 3)
4163d1a8abSmrg * return rounded result
4263d1a8abSmrg * else // large exponent difference
4363d1a8abSmrg * if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
4463d1a8abSmrg * guaranteed the same as
4563d1a8abSmrg * number_digits(coefficient_a*10^diff_expon) )
4663d1a8abSmrg * S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-diff_expon))
4763d1a8abSmrg * corr = 10^16 + (sign_a^sign_b)*coefficient_b
4863d1a8abSmrg * corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
4963d1a8abSmrg * return get_BID64(sign_a,exponent(S),S+rounded(corr))
5063d1a8abSmrg * else
5163d1a8abSmrg * add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
5263d1a8abSmrg * in 128-bit integer arithmetic, then round to 16 decimal digits
5363d1a8abSmrg *
5463d1a8abSmrg *
5563d1a8abSmrg ****************************************************************************/
5663d1a8abSmrg
5763d1a8abSmrg #include "bid_internal.h"
5863d1a8abSmrg
5963d1a8abSmrg #if DECIMAL_CALL_BY_REFERENCE
6063d1a8abSmrg void bid64_add (UINT64 * pres, UINT64 * px,
6163d1a8abSmrg UINT64 *
6263d1a8abSmrg py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
6363d1a8abSmrg _EXC_INFO_PARAM);
6463d1a8abSmrg #else
6563d1a8abSmrg UINT64 bid64_add (UINT64 x,
6663d1a8abSmrg UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
6763d1a8abSmrg _EXC_MASKS_PARAM _EXC_INFO_PARAM);
6863d1a8abSmrg #endif
6963d1a8abSmrg
7063d1a8abSmrg #if DECIMAL_CALL_BY_REFERENCE
7163d1a8abSmrg
7263d1a8abSmrg void
bid64_sub(UINT64 * pres,UINT64 * px,UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)7363d1a8abSmrg bid64_sub (UINT64 * pres, UINT64 * px,
7463d1a8abSmrg UINT64 *
7563d1a8abSmrg py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
7663d1a8abSmrg _EXC_INFO_PARAM) {
7763d1a8abSmrg UINT64 y = *py;
7863d1a8abSmrg #if !DECIMAL_GLOBAL_ROUNDING
7963d1a8abSmrg _IDEC_round rnd_mode = *prnd_mode;
8063d1a8abSmrg #endif
8163d1a8abSmrg // check if y is not NaN
8263d1a8abSmrg if (((y & NAN_MASK64) != NAN_MASK64))
8363d1a8abSmrg y ^= 0x8000000000000000ull;
8463d1a8abSmrg bid64_add (pres, px,
8563d1a8abSmrg &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
8663d1a8abSmrg _EXC_INFO_ARG);
8763d1a8abSmrg }
8863d1a8abSmrg #else
8963d1a8abSmrg
9063d1a8abSmrg UINT64
bid64_sub(UINT64 x,UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)9163d1a8abSmrg bid64_sub (UINT64 x,
9263d1a8abSmrg UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
9363d1a8abSmrg _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
9463d1a8abSmrg // check if y is not NaN
9563d1a8abSmrg if (((y & NAN_MASK64) != NAN_MASK64))
9663d1a8abSmrg y ^= 0x8000000000000000ull;
9763d1a8abSmrg
9863d1a8abSmrg return bid64_add (x,
9963d1a8abSmrg y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
10063d1a8abSmrg _EXC_INFO_ARG);
10163d1a8abSmrg }
10263d1a8abSmrg #endif
10363d1a8abSmrg
10463d1a8abSmrg
10563d1a8abSmrg
10663d1a8abSmrg #if DECIMAL_CALL_BY_REFERENCE
10763d1a8abSmrg
10863d1a8abSmrg void
bid64_add(UINT64 * pres,UINT64 * px,UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)10963d1a8abSmrg bid64_add (UINT64 * pres, UINT64 * px,
11063d1a8abSmrg UINT64 *
11163d1a8abSmrg py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
11263d1a8abSmrg _EXC_INFO_PARAM) {
11363d1a8abSmrg UINT64 x, y;
11463d1a8abSmrg #else
11563d1a8abSmrg
11663d1a8abSmrg UINT64
11763d1a8abSmrg bid64_add (UINT64 x,
11863d1a8abSmrg UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
11963d1a8abSmrg _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
12063d1a8abSmrg #endif
12163d1a8abSmrg
12263d1a8abSmrg UINT128 CA, CT, CT_new;
12363d1a8abSmrg UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
12463d1a8abSmrg UINT64 valid_x, valid_y;
12563d1a8abSmrg UINT64 res;
12663d1a8abSmrg UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
12763d1a8abSmrg rem_a;
12863d1a8abSmrg UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
12963d1a8abSmrg int_double tempx;
13063d1a8abSmrg int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
13163d1a8abSmrg int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
13263d1a8abSmrg unsigned rmode, status;
13363d1a8abSmrg
13463d1a8abSmrg #if DECIMAL_CALL_BY_REFERENCE
13563d1a8abSmrg #if !DECIMAL_GLOBAL_ROUNDING
13663d1a8abSmrg _IDEC_round rnd_mode = *prnd_mode;
13763d1a8abSmrg #endif
13863d1a8abSmrg x = *px;
13963d1a8abSmrg y = *py;
14063d1a8abSmrg #endif
14163d1a8abSmrg
14263d1a8abSmrg valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
14363d1a8abSmrg valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
14463d1a8abSmrg
14563d1a8abSmrg // unpack arguments, check for NaN or Infinity
14663d1a8abSmrg if (!valid_x) {
14763d1a8abSmrg // x is Inf. or NaN
14863d1a8abSmrg
14963d1a8abSmrg // test if x is NaN
15063d1a8abSmrg if ((x & NAN_MASK64) == NAN_MASK64) {
15163d1a8abSmrg #ifdef SET_STATUS_FLAGS
15263d1a8abSmrg if (((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
15363d1a8abSmrg || ((y & SNAN_MASK64) == SNAN_MASK64))
15463d1a8abSmrg __set_status_flags (pfpsf, INVALID_EXCEPTION);
15563d1a8abSmrg #endif
15663d1a8abSmrg res = coefficient_x & QUIET_MASK64;
15763d1a8abSmrg BID_RETURN (res);
15863d1a8abSmrg }
15963d1a8abSmrg // x is Infinity?
16063d1a8abSmrg if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
16163d1a8abSmrg // check if y is Inf
16263d1a8abSmrg if (((y & NAN_MASK64) == INFINITY_MASK64)) {
16363d1a8abSmrg if (sign_x == (y & 0x8000000000000000ull)) {
16463d1a8abSmrg res = coefficient_x;
16563d1a8abSmrg BID_RETURN (res);
16663d1a8abSmrg }
16763d1a8abSmrg // return NaN
16863d1a8abSmrg {
16963d1a8abSmrg #ifdef SET_STATUS_FLAGS
17063d1a8abSmrg __set_status_flags (pfpsf, INVALID_EXCEPTION);
17163d1a8abSmrg #endif
17263d1a8abSmrg res = NAN_MASK64;
17363d1a8abSmrg BID_RETURN (res);
17463d1a8abSmrg }
17563d1a8abSmrg }
17663d1a8abSmrg // check if y is NaN
17763d1a8abSmrg if (((y & NAN_MASK64) == NAN_MASK64)) {
17863d1a8abSmrg res = coefficient_y & QUIET_MASK64;
17963d1a8abSmrg #ifdef SET_STATUS_FLAGS
18063d1a8abSmrg if (((y & SNAN_MASK64) == SNAN_MASK64))
18163d1a8abSmrg __set_status_flags (pfpsf, INVALID_EXCEPTION);
18263d1a8abSmrg #endif
18363d1a8abSmrg BID_RETURN (res);
18463d1a8abSmrg }
18563d1a8abSmrg // otherwise return +/-Inf
18663d1a8abSmrg {
18763d1a8abSmrg res = coefficient_x;
18863d1a8abSmrg BID_RETURN (res);
18963d1a8abSmrg }
19063d1a8abSmrg }
19163d1a8abSmrg // x is 0
19263d1a8abSmrg {
19363d1a8abSmrg if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
19463d1a8abSmrg if (exponent_y <= exponent_x) {
19563d1a8abSmrg res = y;
19663d1a8abSmrg BID_RETURN (res);
19763d1a8abSmrg }
19863d1a8abSmrg }
19963d1a8abSmrg }
20063d1a8abSmrg
20163d1a8abSmrg }
20263d1a8abSmrg if (!valid_y) {
20363d1a8abSmrg // y is Inf. or NaN?
20463d1a8abSmrg if (((y & INFINITY_MASK64) == INFINITY_MASK64)) {
20563d1a8abSmrg #ifdef SET_STATUS_FLAGS
20663d1a8abSmrg if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
20763d1a8abSmrg __set_status_flags (pfpsf, INVALID_EXCEPTION);
20863d1a8abSmrg #endif
20963d1a8abSmrg res = coefficient_y & QUIET_MASK64;
21063d1a8abSmrg BID_RETURN (res);
21163d1a8abSmrg }
21263d1a8abSmrg // y is 0
21363d1a8abSmrg if (!coefficient_x) { // x==0
21463d1a8abSmrg if (exponent_x <= exponent_y)
21563d1a8abSmrg res = ((UINT64) exponent_x) << 53;
21663d1a8abSmrg else
21763d1a8abSmrg res = ((UINT64) exponent_y) << 53;
21863d1a8abSmrg if (sign_x == sign_y)
21963d1a8abSmrg res |= sign_x;
22063d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
22163d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
22263d1a8abSmrg if (rnd_mode == ROUNDING_DOWN && sign_x != sign_y)
22363d1a8abSmrg res |= 0x8000000000000000ull;
22463d1a8abSmrg #endif
22563d1a8abSmrg #endif
22663d1a8abSmrg BID_RETURN (res);
22763d1a8abSmrg } else if (exponent_y >= exponent_x) {
22863d1a8abSmrg res = x;
22963d1a8abSmrg BID_RETURN (res);
23063d1a8abSmrg }
23163d1a8abSmrg }
23263d1a8abSmrg // sort arguments by exponent
23363d1a8abSmrg if (exponent_x < exponent_y) {
23463d1a8abSmrg sign_a = sign_y;
23563d1a8abSmrg exponent_a = exponent_y;
23663d1a8abSmrg coefficient_a = coefficient_y;
23763d1a8abSmrg sign_b = sign_x;
23863d1a8abSmrg exponent_b = exponent_x;
23963d1a8abSmrg coefficient_b = coefficient_x;
24063d1a8abSmrg } else {
24163d1a8abSmrg sign_a = sign_x;
24263d1a8abSmrg exponent_a = exponent_x;
24363d1a8abSmrg coefficient_a = coefficient_x;
24463d1a8abSmrg sign_b = sign_y;
24563d1a8abSmrg exponent_b = exponent_y;
24663d1a8abSmrg coefficient_b = coefficient_y;
24763d1a8abSmrg }
24863d1a8abSmrg
24963d1a8abSmrg // exponent difference
25063d1a8abSmrg diff_dec_expon = exponent_a - exponent_b;
25163d1a8abSmrg
25263d1a8abSmrg /* get binary coefficients of x and y */
25363d1a8abSmrg
25463d1a8abSmrg //--- get number of bits in the coefficients of x and y ---
25563d1a8abSmrg
25663d1a8abSmrg // version 2 (original)
25763d1a8abSmrg tempx.d = (double) coefficient_a;
25863d1a8abSmrg bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
25963d1a8abSmrg
26063d1a8abSmrg if (diff_dec_expon > MAX_FORMAT_DIGITS) {
26163d1a8abSmrg // normalize a to a 16-digit coefficient
26263d1a8abSmrg
26363d1a8abSmrg scale_ca = estimate_decimal_digits[bin_expon_ca];
26463d1a8abSmrg if (coefficient_a >= power10_table_128[scale_ca].w[0])
26563d1a8abSmrg scale_ca++;
26663d1a8abSmrg
26763d1a8abSmrg scale_k = 16 - scale_ca;
26863d1a8abSmrg
26963d1a8abSmrg coefficient_a *= power10_table_128[scale_k].w[0];
27063d1a8abSmrg
27163d1a8abSmrg diff_dec_expon -= scale_k;
27263d1a8abSmrg exponent_a -= scale_k;
27363d1a8abSmrg
27463d1a8abSmrg /* get binary coefficients of x and y */
27563d1a8abSmrg
27663d1a8abSmrg //--- get number of bits in the coefficients of x and y ---
27763d1a8abSmrg tempx.d = (double) coefficient_a;
27863d1a8abSmrg bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
27963d1a8abSmrg
28063d1a8abSmrg if (diff_dec_expon > MAX_FORMAT_DIGITS) {
28163d1a8abSmrg #ifdef SET_STATUS_FLAGS
28263d1a8abSmrg if (coefficient_b) {
28363d1a8abSmrg __set_status_flags (pfpsf, INEXACT_EXCEPTION);
28463d1a8abSmrg }
28563d1a8abSmrg #endif
28663d1a8abSmrg
28763d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
28863d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
28963d1a8abSmrg if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST
29063d1a8abSmrg {
29163d1a8abSmrg switch (rnd_mode) {
29263d1a8abSmrg case ROUNDING_DOWN:
29363d1a8abSmrg if (sign_b) {
29463d1a8abSmrg coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
29563d1a8abSmrg if (coefficient_a < 1000000000000000ull) {
29663d1a8abSmrg exponent_a--;
29763d1a8abSmrg coefficient_a = 9999999999999999ull;
29863d1a8abSmrg } else if (coefficient_a >= 10000000000000000ull) {
29963d1a8abSmrg exponent_a++;
30063d1a8abSmrg coefficient_a = 1000000000000000ull;
30163d1a8abSmrg }
30263d1a8abSmrg }
30363d1a8abSmrg break;
30463d1a8abSmrg case ROUNDING_UP:
30563d1a8abSmrg if (!sign_b) {
30663d1a8abSmrg coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
30763d1a8abSmrg if (coefficient_a < 1000000000000000ull) {
30863d1a8abSmrg exponent_a--;
30963d1a8abSmrg coefficient_a = 9999999999999999ull;
31063d1a8abSmrg } else if (coefficient_a >= 10000000000000000ull) {
31163d1a8abSmrg exponent_a++;
31263d1a8abSmrg coefficient_a = 1000000000000000ull;
31363d1a8abSmrg }
31463d1a8abSmrg }
31563d1a8abSmrg break;
31663d1a8abSmrg default: // RZ
31763d1a8abSmrg if (sign_a != sign_b) {
31863d1a8abSmrg coefficient_a--;
31963d1a8abSmrg if (coefficient_a < 1000000000000000ull) {
32063d1a8abSmrg exponent_a--;
32163d1a8abSmrg coefficient_a = 9999999999999999ull;
32263d1a8abSmrg }
32363d1a8abSmrg }
32463d1a8abSmrg break;
32563d1a8abSmrg }
32663d1a8abSmrg } else
32763d1a8abSmrg #endif
32863d1a8abSmrg #endif
32963d1a8abSmrg // check special case here
33063d1a8abSmrg if ((coefficient_a == 1000000000000000ull)
33163d1a8abSmrg && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
33263d1a8abSmrg && (sign_a ^ sign_b)
33363d1a8abSmrg && (coefficient_b > 5000000000000000ull)) {
33463d1a8abSmrg coefficient_a = 9999999999999999ull;
33563d1a8abSmrg exponent_a--;
33663d1a8abSmrg }
33763d1a8abSmrg
33863d1a8abSmrg res =
33963d1a8abSmrg fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
34063d1a8abSmrg rnd_mode, pfpsf);
34163d1a8abSmrg BID_RETURN (res);
34263d1a8abSmrg }
34363d1a8abSmrg }
34463d1a8abSmrg // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
34563d1a8abSmrg if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
34663d1a8abSmrg // coefficient_a*10^(exponent_a-exponent_b)<2^63
34763d1a8abSmrg
34863d1a8abSmrg // multiply by 10^(exponent_a-exponent_b)
34963d1a8abSmrg coefficient_a *= power10_table_128[diff_dec_expon].w[0];
35063d1a8abSmrg
35163d1a8abSmrg // sign mask
35263d1a8abSmrg sign_b = ((SINT64) sign_b) >> 63;
35363d1a8abSmrg // apply sign to coeff. of b
35463d1a8abSmrg coefficient_b = (coefficient_b + sign_b) ^ sign_b;
35563d1a8abSmrg
35663d1a8abSmrg // apply sign to coefficient a
35763d1a8abSmrg sign_a = ((SINT64) sign_a) >> 63;
35863d1a8abSmrg coefficient_a = (coefficient_a + sign_a) ^ sign_a;
35963d1a8abSmrg
36063d1a8abSmrg coefficient_a += coefficient_b;
36163d1a8abSmrg // get sign
36263d1a8abSmrg sign_s = ((SINT64) coefficient_a) >> 63;
36363d1a8abSmrg coefficient_a = (coefficient_a + sign_s) ^ sign_s;
36463d1a8abSmrg sign_s &= 0x8000000000000000ull;
36563d1a8abSmrg
36663d1a8abSmrg // coefficient_a < 10^16 ?
36763d1a8abSmrg if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
36863d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
36963d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
37063d1a8abSmrg if (rnd_mode == ROUNDING_DOWN && (!coefficient_a)
37163d1a8abSmrg && sign_a != sign_b)
37263d1a8abSmrg sign_s = 0x8000000000000000ull;
37363d1a8abSmrg #endif
37463d1a8abSmrg #endif
37563d1a8abSmrg res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
37663d1a8abSmrg BID_RETURN (res);
37763d1a8abSmrg }
37863d1a8abSmrg // otherwise rounding is necessary
37963d1a8abSmrg
38063d1a8abSmrg // already know coefficient_a<10^19
38163d1a8abSmrg // coefficient_a < 10^17 ?
38263d1a8abSmrg if (coefficient_a < power10_table_128[17].w[0])
38363d1a8abSmrg extra_digits = 1;
38463d1a8abSmrg else if (coefficient_a < power10_table_128[18].w[0])
38563d1a8abSmrg extra_digits = 2;
38663d1a8abSmrg else
38763d1a8abSmrg extra_digits = 3;
38863d1a8abSmrg
38963d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
39063d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
39163d1a8abSmrg rmode = rnd_mode;
39263d1a8abSmrg if (sign_s && (unsigned) (rmode - 1) < 2)
39363d1a8abSmrg rmode = 3 - rmode;
39463d1a8abSmrg #else
39563d1a8abSmrg rmode = 0;
39663d1a8abSmrg #endif
39763d1a8abSmrg #else
39863d1a8abSmrg rmode = 0;
39963d1a8abSmrg #endif
40063d1a8abSmrg coefficient_a += round_const_table[rmode][extra_digits];
40163d1a8abSmrg
40263d1a8abSmrg // get P*(2^M[extra_digits])/10^extra_digits
40363d1a8abSmrg __mul_64x64_to_128 (CT, coefficient_a,
40463d1a8abSmrg reciprocals10_64[extra_digits]);
40563d1a8abSmrg
40663d1a8abSmrg // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
40763d1a8abSmrg amount = short_recip_scale[extra_digits];
40863d1a8abSmrg C64 = CT.w[1] >> amount;
40963d1a8abSmrg
41063d1a8abSmrg } else {
41163d1a8abSmrg // coefficient_a*10^(exponent_a-exponent_b) is large
41263d1a8abSmrg sign_s = sign_a;
41363d1a8abSmrg
41463d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
41563d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
41663d1a8abSmrg rmode = rnd_mode;
41763d1a8abSmrg if (sign_s && (unsigned) (rmode - 1) < 2)
41863d1a8abSmrg rmode = 3 - rmode;
41963d1a8abSmrg #else
42063d1a8abSmrg rmode = 0;
42163d1a8abSmrg #endif
42263d1a8abSmrg #else
42363d1a8abSmrg rmode = 0;
42463d1a8abSmrg #endif
42563d1a8abSmrg
42663d1a8abSmrg // check whether we can take faster path
42763d1a8abSmrg scale_ca = estimate_decimal_digits[bin_expon_ca];
42863d1a8abSmrg
42963d1a8abSmrg sign_ab = sign_a ^ sign_b;
43063d1a8abSmrg sign_ab = ((SINT64) sign_ab) >> 63;
43163d1a8abSmrg
43263d1a8abSmrg // T1 = 10^(16-diff_dec_expon)
43363d1a8abSmrg T1 = power10_table_128[16 - diff_dec_expon].w[0];
43463d1a8abSmrg
43563d1a8abSmrg // get number of digits in coefficient_a
43663d1a8abSmrg if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
43763d1a8abSmrg scale_ca++;
43863d1a8abSmrg }
43963d1a8abSmrg
44063d1a8abSmrg scale_k = 16 - scale_ca;
44163d1a8abSmrg
44263d1a8abSmrg // addition
44363d1a8abSmrg saved_ca = coefficient_a - T1;
44463d1a8abSmrg coefficient_a =
44563d1a8abSmrg (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
44663d1a8abSmrg extra_digits = diff_dec_expon - scale_k;
44763d1a8abSmrg
44863d1a8abSmrg // apply sign
44963d1a8abSmrg saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
45063d1a8abSmrg // add 10^16 and rounding constant
45163d1a8abSmrg coefficient_b =
45263d1a8abSmrg saved_cb + 10000000000000000ull +
45363d1a8abSmrg round_const_table[rmode][extra_digits];
45463d1a8abSmrg
45563d1a8abSmrg // get P*(2^M[extra_digits])/10^extra_digits
45663d1a8abSmrg __mul_64x64_to_128 (CT, coefficient_b,
45763d1a8abSmrg reciprocals10_64[extra_digits]);
45863d1a8abSmrg
45963d1a8abSmrg // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
46063d1a8abSmrg amount = short_recip_scale[extra_digits];
46163d1a8abSmrg C0_64 = CT.w[1] >> amount;
46263d1a8abSmrg
46363d1a8abSmrg // result coefficient
46463d1a8abSmrg C64 = C0_64 + coefficient_a;
46563d1a8abSmrg // filter out difficult (corner) cases
46663d1a8abSmrg // this test ensures the number of digits in coefficient_a does not change
46763d1a8abSmrg // after adding (the appropriately scaled and rounded) coefficient_b
46863d1a8abSmrg if ((UINT64) (C64 - 1000000000000000ull - 1) >
46963d1a8abSmrg 9000000000000000ull - 2) {
47063d1a8abSmrg if (C64 >= 10000000000000000ull) {
47163d1a8abSmrg // result has more than 16 digits
47263d1a8abSmrg if (!scale_k) {
47363d1a8abSmrg // must divide coeff_a by 10
47463d1a8abSmrg saved_ca = saved_ca + T1;
47563d1a8abSmrg __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
47663d1a8abSmrg //reciprocals10_64[1]);
47763d1a8abSmrg coefficient_a = CA.w[1] >> 1;
47863d1a8abSmrg rem_a =
47963d1a8abSmrg saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
48063d1a8abSmrg coefficient_a = coefficient_a - T1;
48163d1a8abSmrg
48263d1a8abSmrg saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
48363d1a8abSmrg } else
48463d1a8abSmrg coefficient_a =
48563d1a8abSmrg (SINT64) (saved_ca - T1 -
48663d1a8abSmrg (T1 << 3)) * (SINT64) power10_table_128[scale_k -
48763d1a8abSmrg 1].w[0];
48863d1a8abSmrg
48963d1a8abSmrg extra_digits++;
49063d1a8abSmrg coefficient_b =
49163d1a8abSmrg saved_cb + 100000000000000000ull +
49263d1a8abSmrg round_const_table[rmode][extra_digits];
49363d1a8abSmrg
49463d1a8abSmrg // get P*(2^M[extra_digits])/10^extra_digits
49563d1a8abSmrg __mul_64x64_to_128 (CT, coefficient_b,
49663d1a8abSmrg reciprocals10_64[extra_digits]);
49763d1a8abSmrg
49863d1a8abSmrg // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
49963d1a8abSmrg amount = short_recip_scale[extra_digits];
50063d1a8abSmrg C0_64 = CT.w[1] >> amount;
50163d1a8abSmrg
50263d1a8abSmrg // result coefficient
50363d1a8abSmrg C64 = C0_64 + coefficient_a;
50463d1a8abSmrg } else if (C64 <= 1000000000000000ull) {
50563d1a8abSmrg // less than 16 digits in result
50663d1a8abSmrg coefficient_a =
50763d1a8abSmrg (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
50863d1a8abSmrg 1].w[0];
50963d1a8abSmrg //extra_digits --;
51063d1a8abSmrg exponent_b--;
51163d1a8abSmrg coefficient_b =
51263d1a8abSmrg (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
51363d1a8abSmrg round_const_table[rmode][extra_digits];
51463d1a8abSmrg
51563d1a8abSmrg // get P*(2^M[extra_digits])/10^extra_digits
51663d1a8abSmrg __mul_64x64_to_128 (CT_new, coefficient_b,
51763d1a8abSmrg reciprocals10_64[extra_digits]);
51863d1a8abSmrg
51963d1a8abSmrg // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
52063d1a8abSmrg amount = short_recip_scale[extra_digits];
52163d1a8abSmrg C0_64 = CT_new.w[1] >> amount;
52263d1a8abSmrg
52363d1a8abSmrg // result coefficient
52463d1a8abSmrg C64_new = C0_64 + coefficient_a;
52563d1a8abSmrg if (C64_new < 10000000000000000ull) {
52663d1a8abSmrg C64 = C64_new;
52763d1a8abSmrg #ifdef SET_STATUS_FLAGS
52863d1a8abSmrg CT = CT_new;
52963d1a8abSmrg #endif
53063d1a8abSmrg } else
53163d1a8abSmrg exponent_b++;
53263d1a8abSmrg }
53363d1a8abSmrg
53463d1a8abSmrg }
53563d1a8abSmrg
53663d1a8abSmrg }
53763d1a8abSmrg
53863d1a8abSmrg #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
53963d1a8abSmrg #ifndef IEEE_ROUND_NEAREST
54063d1a8abSmrg if (rmode == 0) //ROUNDING_TO_NEAREST
54163d1a8abSmrg #endif
54263d1a8abSmrg if (C64 & 1) {
54363d1a8abSmrg // check whether fractional part of initial_P/10^extra_digits is
54463d1a8abSmrg // exactly .5
54563d1a8abSmrg // this is the same as fractional part of
54663d1a8abSmrg // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
54763d1a8abSmrg
54863d1a8abSmrg // get remainder
54963d1a8abSmrg remainder_h = CT.w[1] << (64 - amount);
55063d1a8abSmrg
55163d1a8abSmrg // test whether fractional part is 0
55263d1a8abSmrg if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
55363d1a8abSmrg C64--;
55463d1a8abSmrg }
55563d1a8abSmrg }
55663d1a8abSmrg #endif
55763d1a8abSmrg
55863d1a8abSmrg #ifdef SET_STATUS_FLAGS
55963d1a8abSmrg status = INEXACT_EXCEPTION;
56063d1a8abSmrg
56163d1a8abSmrg // get remainder
56263d1a8abSmrg remainder_h = CT.w[1] << (64 - amount);
56363d1a8abSmrg
56463d1a8abSmrg switch (rmode) {
56563d1a8abSmrg case ROUNDING_TO_NEAREST:
56663d1a8abSmrg case ROUNDING_TIES_AWAY:
56763d1a8abSmrg // test whether fractional part is 0
56863d1a8abSmrg if ((remainder_h == 0x8000000000000000ull)
56963d1a8abSmrg && (CT.w[0] < reciprocals10_64[extra_digits]))
57063d1a8abSmrg status = EXACT_STATUS;
57163d1a8abSmrg break;
57263d1a8abSmrg case ROUNDING_DOWN:
57363d1a8abSmrg case ROUNDING_TO_ZERO:
57463d1a8abSmrg if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
57563d1a8abSmrg status = EXACT_STATUS;
57663d1a8abSmrg //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
57763d1a8abSmrg break;
57863d1a8abSmrg default:
57963d1a8abSmrg // round up
58063d1a8abSmrg __add_carry_out (tmp, carry, CT.w[0],
58163d1a8abSmrg reciprocals10_64[extra_digits]);
58263d1a8abSmrg if ((remainder_h >> (64 - amount)) + carry >=
58363d1a8abSmrg (((UINT64) 1) << amount))
58463d1a8abSmrg status = EXACT_STATUS;
58563d1a8abSmrg break;
58663d1a8abSmrg }
58763d1a8abSmrg __set_status_flags (pfpsf, status);
58863d1a8abSmrg
58963d1a8abSmrg #endif
59063d1a8abSmrg
59163d1a8abSmrg res =
59263d1a8abSmrg fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
59363d1a8abSmrg rnd_mode, pfpsf);
59463d1a8abSmrg BID_RETURN (res);
59563d1a8abSmrg }
596