/* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
/* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
|
|
This file is part of GCC.
|
This file is part of GCC.
|
|
|
GCC is free software; you can redistribute it and/or modify it under
|
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
|
the terms of the GNU General Public License as published by the Free
|
Software Foundation; either version 3, or (at your option) any later
|
Software Foundation; either version 3, or (at your option) any later
|
version.
|
version.
|
|
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
for more details.
|
for more details.
|
|
|
Under Section 7 of GPL version 3, you are granted additional
|
Under Section 7 of GPL version 3, you are granted additional
|
permissions described in the GCC Runtime Library Exception, version
|
permissions described in the GCC Runtime Library Exception, version
|
3.1, as published by the Free Software Foundation.
|
3.1, as published by the Free Software Foundation.
|
|
|
You should have received a copy of the GNU General Public License and
|
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;
|
a copy of the GCC Runtime Library Exception along with this program;
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
<http://www.gnu.org/licenses/>. */
|
<http://www.gnu.org/licenses/>. */
|
|
|
/*****************************************************************************
|
/*****************************************************************************
|
* BID64 add
|
* BID64 add
|
*****************************************************************************
|
*****************************************************************************
|
*
|
*
|
* Algorithm description:
|
* Algorithm description:
|
*
|
*
|
* if(exponent_a < exponent_b)
|
* if(exponent_a < exponent_b)
|
* switch a, b
|
* switch a, b
|
* diff_expon = exponent_a - exponent_b
|
* diff_expon = exponent_a - exponent_b
|
* if(diff_expon > 16)
|
* if(diff_expon > 16)
|
* return normalize(a)
|
* return normalize(a)
|
* if(coefficient_a*10^diff_expon guaranteed below 2^62)
|
* if(coefficient_a*10^diff_expon guaranteed below 2^62)
|
* S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
|
* S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
|
* if(|S|<10^16)
|
* if(|S|<10^16)
|
* return get_BID64(sign(S),exponent_b,|S|)
|
* return get_BID64(sign(S),exponent_b,|S|)
|
* else
|
* else
|
* determine number of extra digits in S (1, 2, or 3)
|
* determine number of extra digits in S (1, 2, or 3)
|
* return rounded result
|
* return rounded result
|
* else // large exponent difference
|
* else // large exponent difference
|
* if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
|
* if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
|
* guaranteed the same as
|
* guaranteed the same as
|
* number_digits(coefficient_a*10^diff_expon) )
|
* number_digits(coefficient_a*10^diff_expon) )
|
* S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-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^16 + (sign_a^sign_b)*coefficient_b
|
* corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
|
* corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
|
* return get_BID64(sign_a,exponent(S),S+rounded(corr))
|
* return get_BID64(sign_a,exponent(S),S+rounded(corr))
|
* else
|
* else
|
* add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
|
* add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
|
* in 128-bit integer arithmetic, then round to 16 decimal digits
|
* in 128-bit integer arithmetic, then round to 16 decimal digits
|
*
|
*
|
*
|
*
|
****************************************************************************/
|
****************************************************************************/
|
|
|
#include "bid_internal.h"
|
#include "bid_internal.h"
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
#if DECIMAL_CALL_BY_REFERENCE
|
void bid64_add (UINT64 * pres, UINT64 * px,
|
void bid64_add (UINT64 * pres, UINT64 * px,
|
UINT64 *
|
UINT64 *
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
_EXC_INFO_PARAM);
|
_EXC_INFO_PARAM);
|
#else
|
#else
|
UINT64 bid64_add (UINT64 x,
|
UINT64 bid64_add (UINT64 x,
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM);
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM);
|
#endif
|
#endif
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
|
void
|
void
|
bid64_sub (UINT64 * pres, UINT64 * px,
|
bid64_sub (UINT64 * pres, UINT64 * px,
|
UINT64 *
|
UINT64 *
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
_EXC_INFO_PARAM) {
|
_EXC_INFO_PARAM) {
|
UINT64 y = *py;
|
UINT64 y = *py;
|
#if !DECIMAL_GLOBAL_ROUNDING
|
#if !DECIMAL_GLOBAL_ROUNDING
|
_IDEC_round rnd_mode = *prnd_mode;
|
_IDEC_round rnd_mode = *prnd_mode;
|
#endif
|
#endif
|
// check if y is not NaN
|
// check if y is not NaN
|
if (((y & NAN_MASK64) != NAN_MASK64))
|
if (((y & NAN_MASK64) != NAN_MASK64))
|
y ^= 0x8000000000000000ull;
|
y ^= 0x8000000000000000ull;
|
bid64_add (pres, px,
|
bid64_add (pres, px,
|
&y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
&y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
_EXC_INFO_ARG);
|
_EXC_INFO_ARG);
|
}
|
}
|
#else
|
#else
|
|
|
UINT64
|
UINT64
|
bid64_sub (UINT64 x,
|
bid64_sub (UINT64 x,
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
// check if y is not NaN
|
// check if y is not NaN
|
if (((y & NAN_MASK64) != NAN_MASK64))
|
if (((y & NAN_MASK64) != NAN_MASK64))
|
y ^= 0x8000000000000000ull;
|
y ^= 0x8000000000000000ull;
|
|
|
return bid64_add (x,
|
return bid64_add (x,
|
y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
_EXC_INFO_ARG);
|
_EXC_INFO_ARG);
|
}
|
}
|
#endif
|
#endif
|
|
|
|
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
|
void
|
void
|
bid64_add (UINT64 * pres, UINT64 * px,
|
bid64_add (UINT64 * pres, UINT64 * px,
|
UINT64 *
|
UINT64 *
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
_EXC_INFO_PARAM) {
|
_EXC_INFO_PARAM) {
|
UINT64 x, y;
|
UINT64 x, y;
|
#else
|
#else
|
|
|
UINT64
|
UINT64
|
bid64_add (UINT64 x,
|
bid64_add (UINT64 x,
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
#endif
|
#endif
|
|
|
UINT128 CA, CT, CT_new;
|
UINT128 CA, CT, CT_new;
|
UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
|
UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
|
UINT64 valid_x, valid_y;
|
UINT64 valid_x, valid_y;
|
UINT64 res;
|
UINT64 res;
|
UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
|
UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
|
rem_a;
|
rem_a;
|
UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
|
UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
|
int_double tempx;
|
int_double tempx;
|
int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
|
int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
|
int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
|
int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
|
unsigned rmode, status;
|
unsigned rmode, status;
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
#if DECIMAL_CALL_BY_REFERENCE
|
#if !DECIMAL_GLOBAL_ROUNDING
|
#if !DECIMAL_GLOBAL_ROUNDING
|
_IDEC_round rnd_mode = *prnd_mode;
|
_IDEC_round rnd_mode = *prnd_mode;
|
#endif
|
#endif
|
x = *px;
|
x = *px;
|
y = *py;
|
y = *py;
|
#endif
|
#endif
|
|
|
valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
|
valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
|
|
|
// unpack arguments, check for NaN or Infinity
|
// unpack arguments, check for NaN or Infinity
|
if (!valid_x) {
|
if (!valid_x) {
|
// x is Inf. or NaN
|
// x is Inf. or NaN
|
|
|
// test if x is NaN
|
// test if x is NaN
|
if ((x & NAN_MASK64) == NAN_MASK64) {
|
if ((x & NAN_MASK64) == NAN_MASK64) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
if (((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
|| ((y & SNAN_MASK64) == SNAN_MASK64))
|
|| ((y & SNAN_MASK64) == SNAN_MASK64))
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res = coefficient_x & QUIET_MASK64;
|
res = coefficient_x & QUIET_MASK64;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// x is Infinity?
|
// x is Infinity?
|
if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
|
if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
|
// check if y is Inf
|
// check if y is Inf
|
if (((y & NAN_MASK64) == INFINITY_MASK64)) {
|
if (((y & NAN_MASK64) == INFINITY_MASK64)) {
|
if (sign_x == (y & 0x8000000000000000ull)) {
|
if (sign_x == (y & 0x8000000000000000ull)) {
|
res = coefficient_x;
|
res = coefficient_x;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// return NaN
|
// return NaN
|
{
|
{
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res = NAN_MASK64;
|
res = NAN_MASK64;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// check if y is NaN
|
// check if y is NaN
|
if (((y & NAN_MASK64) == NAN_MASK64)) {
|
if (((y & NAN_MASK64) == NAN_MASK64)) {
|
res = coefficient_y & QUIET_MASK64;
|
res = coefficient_y & QUIET_MASK64;
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (((y & SNAN_MASK64) == SNAN_MASK64))
|
if (((y & SNAN_MASK64) == SNAN_MASK64))
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// otherwise return +/-Inf
|
// otherwise return +/-Inf
|
{
|
{
|
res = coefficient_x;
|
res = coefficient_x;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// x is 0
|
// x is 0
|
{
|
{
|
if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
|
if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
|
if (exponent_y <= exponent_x) {
|
if (exponent_y <= exponent_x) {
|
res = y;
|
res = y;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
}
|
}
|
|
|
}
|
}
|
if (!valid_y) {
|
if (!valid_y) {
|
// y is Inf. or NaN?
|
// y is Inf. or NaN?
|
if (((y & INFINITY_MASK64) == INFINITY_MASK64)) {
|
if (((y & INFINITY_MASK64) == INFINITY_MASK64)) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res = coefficient_y & QUIET_MASK64;
|
res = coefficient_y & QUIET_MASK64;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is 0
|
// y is 0
|
if (!coefficient_x) { // x==0
|
if (!coefficient_x) { // x==0
|
if (exponent_x <= exponent_y)
|
if (exponent_x <= exponent_y)
|
res = ((UINT64) exponent_x) << 53;
|
res = ((UINT64) exponent_x) << 53;
|
else
|
else
|
res = ((UINT64) exponent_y) << 53;
|
res = ((UINT64) exponent_y) << 53;
|
if (sign_x == sign_y)
|
if (sign_x == sign_y)
|
res |= sign_x;
|
res |= sign_x;
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
if (rnd_mode == ROUNDING_DOWN && sign_x != sign_y)
|
if (rnd_mode == ROUNDING_DOWN && sign_x != sign_y)
|
res |= 0x8000000000000000ull;
|
res |= 0x8000000000000000ull;
|
#endif
|
#endif
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
} else if (exponent_y >= exponent_x) {
|
} else if (exponent_y >= exponent_x) {
|
res = x;
|
res = x;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// sort arguments by exponent
|
// sort arguments by exponent
|
if (exponent_x < exponent_y) {
|
if (exponent_x < exponent_y) {
|
sign_a = sign_y;
|
sign_a = sign_y;
|
exponent_a = exponent_y;
|
exponent_a = exponent_y;
|
coefficient_a = coefficient_y;
|
coefficient_a = coefficient_y;
|
sign_b = sign_x;
|
sign_b = sign_x;
|
exponent_b = exponent_x;
|
exponent_b = exponent_x;
|
coefficient_b = coefficient_x;
|
coefficient_b = coefficient_x;
|
} else {
|
} else {
|
sign_a = sign_x;
|
sign_a = sign_x;
|
exponent_a = exponent_x;
|
exponent_a = exponent_x;
|
coefficient_a = coefficient_x;
|
coefficient_a = coefficient_x;
|
sign_b = sign_y;
|
sign_b = sign_y;
|
exponent_b = exponent_y;
|
exponent_b = exponent_y;
|
coefficient_b = coefficient_y;
|
coefficient_b = coefficient_y;
|
}
|
}
|
|
|
// exponent difference
|
// exponent difference
|
diff_dec_expon = exponent_a - exponent_b;
|
diff_dec_expon = exponent_a - exponent_b;
|
|
|
/* get binary coefficients of x and y */
|
/* get binary coefficients of x and y */
|
|
|
//--- get number of bits in the coefficients of x and y ---
|
//--- get number of bits in the coefficients of x and y ---
|
|
|
// version 2 (original)
|
// version 2 (original)
|
tempx.d = (double) coefficient_a;
|
tempx.d = (double) coefficient_a;
|
bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
|
|
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
// normalize a to a 16-digit coefficient
|
// normalize a to a 16-digit coefficient
|
|
|
scale_ca = estimate_decimal_digits[bin_expon_ca];
|
scale_ca = estimate_decimal_digits[bin_expon_ca];
|
if (coefficient_a >= power10_table_128[scale_ca].w[0])
|
if (coefficient_a >= power10_table_128[scale_ca].w[0])
|
scale_ca++;
|
scale_ca++;
|
|
|
scale_k = 16 - scale_ca;
|
scale_k = 16 - scale_ca;
|
|
|
coefficient_a *= power10_table_128[scale_k].w[0];
|
coefficient_a *= power10_table_128[scale_k].w[0];
|
|
|
diff_dec_expon -= scale_k;
|
diff_dec_expon -= scale_k;
|
exponent_a -= scale_k;
|
exponent_a -= scale_k;
|
|
|
/* get binary coefficients of x and y */
|
/* get binary coefficients of x and y */
|
|
|
//--- get number of bits in the coefficients of x and y ---
|
//--- get number of bits in the coefficients of x and y ---
|
tempx.d = (double) coefficient_a;
|
tempx.d = (double) coefficient_a;
|
bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
|
|
|
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
if (diff_dec_expon > MAX_FORMAT_DIGITS) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (coefficient_b) {
|
if (coefficient_b) {
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#endif
|
#endif
|
|
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST
|
if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST
|
{
|
{
|
switch (rnd_mode) {
|
switch (rnd_mode) {
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
if (sign_b) {
|
if (sign_b) {
|
coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
|
coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
|
if (coefficient_a < 1000000000000000ull) {
|
if (coefficient_a < 1000000000000000ull) {
|
exponent_a--;
|
exponent_a--;
|
coefficient_a = 9999999999999999ull;
|
coefficient_a = 9999999999999999ull;
|
} else if (coefficient_a >= 10000000000000000ull) {
|
} else if (coefficient_a >= 10000000000000000ull) {
|
exponent_a++;
|
exponent_a++;
|
coefficient_a = 1000000000000000ull;
|
coefficient_a = 1000000000000000ull;
|
}
|
}
|
}
|
}
|
break;
|
break;
|
case ROUNDING_UP:
|
case ROUNDING_UP:
|
if (!sign_b) {
|
if (!sign_b) {
|
coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
|
coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
|
if (coefficient_a < 1000000000000000ull) {
|
if (coefficient_a < 1000000000000000ull) {
|
exponent_a--;
|
exponent_a--;
|
coefficient_a = 9999999999999999ull;
|
coefficient_a = 9999999999999999ull;
|
} else if (coefficient_a >= 10000000000000000ull) {
|
} else if (coefficient_a >= 10000000000000000ull) {
|
exponent_a++;
|
exponent_a++;
|
coefficient_a = 1000000000000000ull;
|
coefficient_a = 1000000000000000ull;
|
}
|
}
|
}
|
}
|
break;
|
break;
|
default: // RZ
|
default: // RZ
|
if (sign_a != sign_b) {
|
if (sign_a != sign_b) {
|
coefficient_a--;
|
coefficient_a--;
|
if (coefficient_a < 1000000000000000ull) {
|
if (coefficient_a < 1000000000000000ull) {
|
exponent_a--;
|
exponent_a--;
|
coefficient_a = 9999999999999999ull;
|
coefficient_a = 9999999999999999ull;
|
}
|
}
|
}
|
}
|
break;
|
break;
|
}
|
}
|
} else
|
} else
|
#endif
|
#endif
|
#endif
|
#endif
|
// check special case here
|
// check special case here
|
if ((coefficient_a == 1000000000000000ull)
|
if ((coefficient_a == 1000000000000000ull)
|
&& (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
|
&& (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
|
&& (sign_a ^ sign_b)
|
&& (sign_a ^ sign_b)
|
&& (coefficient_b > 5000000000000000ull)) {
|
&& (coefficient_b > 5000000000000000ull)) {
|
coefficient_a = 9999999999999999ull;
|
coefficient_a = 9999999999999999ull;
|
exponent_a--;
|
exponent_a--;
|
}
|
}
|
|
|
res =
|
res =
|
fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
|
fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
|
rnd_mode, pfpsf);
|
rnd_mode, pfpsf);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
|
// test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
|
if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
|
if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
|
// coefficient_a*10^(exponent_a-exponent_b)<2^63
|
// coefficient_a*10^(exponent_a-exponent_b)<2^63
|
|
|
// multiply by 10^(exponent_a-exponent_b)
|
// multiply by 10^(exponent_a-exponent_b)
|
coefficient_a *= power10_table_128[diff_dec_expon].w[0];
|
coefficient_a *= power10_table_128[diff_dec_expon].w[0];
|
|
|
// sign mask
|
// sign mask
|
sign_b = ((SINT64) sign_b) >> 63;
|
sign_b = ((SINT64) sign_b) >> 63;
|
// apply sign to coeff. of b
|
// apply sign to coeff. of b
|
coefficient_b = (coefficient_b + sign_b) ^ sign_b;
|
coefficient_b = (coefficient_b + sign_b) ^ sign_b;
|
|
|
// apply sign to coefficient a
|
// apply sign to coefficient a
|
sign_a = ((SINT64) sign_a) >> 63;
|
sign_a = ((SINT64) sign_a) >> 63;
|
coefficient_a = (coefficient_a + sign_a) ^ sign_a;
|
coefficient_a = (coefficient_a + sign_a) ^ sign_a;
|
|
|
coefficient_a += coefficient_b;
|
coefficient_a += coefficient_b;
|
// get sign
|
// get sign
|
sign_s = ((SINT64) coefficient_a) >> 63;
|
sign_s = ((SINT64) coefficient_a) >> 63;
|
coefficient_a = (coefficient_a + sign_s) ^ sign_s;
|
coefficient_a = (coefficient_a + sign_s) ^ sign_s;
|
sign_s &= 0x8000000000000000ull;
|
sign_s &= 0x8000000000000000ull;
|
|
|
// coefficient_a < 10^16 ?
|
// coefficient_a < 10^16 ?
|
if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
|
if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
if (rnd_mode == ROUNDING_DOWN && (!coefficient_a)
|
if (rnd_mode == ROUNDING_DOWN && (!coefficient_a)
|
&& sign_a != sign_b)
|
&& sign_a != sign_b)
|
sign_s = 0x8000000000000000ull;
|
sign_s = 0x8000000000000000ull;
|
#endif
|
#endif
|
#endif
|
#endif
|
res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
|
res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// otherwise rounding is necessary
|
// otherwise rounding is necessary
|
|
|
// already know coefficient_a<10^19
|
// already know coefficient_a<10^19
|
// coefficient_a < 10^17 ?
|
// coefficient_a < 10^17 ?
|
if (coefficient_a < power10_table_128[17].w[0])
|
if (coefficient_a < power10_table_128[17].w[0])
|
extra_digits = 1;
|
extra_digits = 1;
|
else if (coefficient_a < power10_table_128[18].w[0])
|
else if (coefficient_a < power10_table_128[18].w[0])
|
extra_digits = 2;
|
extra_digits = 2;
|
else
|
else
|
extra_digits = 3;
|
extra_digits = 3;
|
|
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_s && (unsigned) (rmode - 1) < 2)
|
if (sign_s && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
#else
|
#else
|
rmode = 0;
|
rmode = 0;
|
#endif
|
#endif
|
#else
|
#else
|
rmode = 0;
|
rmode = 0;
|
#endif
|
#endif
|
coefficient_a += round_const_table[rmode][extra_digits];
|
coefficient_a += round_const_table[rmode][extra_digits];
|
|
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_64x64_to_128 (CT, coefficient_a,
|
__mul_64x64_to_128 (CT, coefficient_a,
|
reciprocals10_64[extra_digits]);
|
reciprocals10_64[extra_digits]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
amount = short_recip_scale[extra_digits];
|
amount = short_recip_scale[extra_digits];
|
C64 = CT.w[1] >> amount;
|
C64 = CT.w[1] >> amount;
|
|
|
} else {
|
} else {
|
// coefficient_a*10^(exponent_a-exponent_b) is large
|
// coefficient_a*10^(exponent_a-exponent_b) is large
|
sign_s = sign_a;
|
sign_s = sign_a;
|
|
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_s && (unsigned) (rmode - 1) < 2)
|
if (sign_s && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
#else
|
#else
|
rmode = 0;
|
rmode = 0;
|
#endif
|
#endif
|
#else
|
#else
|
rmode = 0;
|
rmode = 0;
|
#endif
|
#endif
|
|
|
// check whether we can take faster path
|
// check whether we can take faster path
|
scale_ca = estimate_decimal_digits[bin_expon_ca];
|
scale_ca = estimate_decimal_digits[bin_expon_ca];
|
|
|
sign_ab = sign_a ^ sign_b;
|
sign_ab = sign_a ^ sign_b;
|
sign_ab = ((SINT64) sign_ab) >> 63;
|
sign_ab = ((SINT64) sign_ab) >> 63;
|
|
|
// T1 = 10^(16-diff_dec_expon)
|
// T1 = 10^(16-diff_dec_expon)
|
T1 = power10_table_128[16 - diff_dec_expon].w[0];
|
T1 = power10_table_128[16 - diff_dec_expon].w[0];
|
|
|
// get number of digits in coefficient_a
|
// get number of digits in coefficient_a
|
if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
|
if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
|
scale_ca++;
|
scale_ca++;
|
}
|
}
|
|
|
scale_k = 16 - scale_ca;
|
scale_k = 16 - scale_ca;
|
|
|
// addition
|
// addition
|
saved_ca = coefficient_a - T1;
|
saved_ca = coefficient_a - T1;
|
coefficient_a =
|
coefficient_a =
|
(SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
|
(SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
|
extra_digits = diff_dec_expon - scale_k;
|
extra_digits = diff_dec_expon - scale_k;
|
|
|
// apply sign
|
// apply sign
|
saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
|
saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
|
// add 10^16 and rounding constant
|
// add 10^16 and rounding constant
|
coefficient_b =
|
coefficient_b =
|
saved_cb + 10000000000000000ull +
|
saved_cb + 10000000000000000ull +
|
round_const_table[rmode][extra_digits];
|
round_const_table[rmode][extra_digits];
|
|
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_64x64_to_128 (CT, coefficient_b,
|
__mul_64x64_to_128 (CT, coefficient_b,
|
reciprocals10_64[extra_digits]);
|
reciprocals10_64[extra_digits]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
amount = short_recip_scale[extra_digits];
|
amount = short_recip_scale[extra_digits];
|
C0_64 = CT.w[1] >> amount;
|
C0_64 = CT.w[1] >> amount;
|
|
|
// result coefficient
|
// result coefficient
|
C64 = C0_64 + coefficient_a;
|
C64 = C0_64 + coefficient_a;
|
// filter out difficult (corner) cases
|
// filter out difficult (corner) cases
|
// this test ensures the number of digits in coefficient_a does not change
|
// this test ensures the number of digits in coefficient_a does not change
|
// after adding (the appropriately scaled and rounded) coefficient_b
|
// after adding (the appropriately scaled and rounded) coefficient_b
|
if ((UINT64) (C64 - 1000000000000000ull - 1) >
|
if ((UINT64) (C64 - 1000000000000000ull - 1) >
|
9000000000000000ull - 2) {
|
9000000000000000ull - 2) {
|
if (C64 >= 10000000000000000ull) {
|
if (C64 >= 10000000000000000ull) {
|
// result has more than 16 digits
|
// result has more than 16 digits
|
if (!scale_k) {
|
if (!scale_k) {
|
// must divide coeff_a by 10
|
// must divide coeff_a by 10
|
saved_ca = saved_ca + T1;
|
saved_ca = saved_ca + T1;
|
__mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
|
__mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
|
//reciprocals10_64[1]);
|
//reciprocals10_64[1]);
|
coefficient_a = CA.w[1] >> 1;
|
coefficient_a = CA.w[1] >> 1;
|
rem_a =
|
rem_a =
|
saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
|
saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
|
coefficient_a = coefficient_a - T1;
|
coefficient_a = coefficient_a - T1;
|
|
|
saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
|
saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
|
} else
|
} else
|
coefficient_a =
|
coefficient_a =
|
(SINT64) (saved_ca - T1 -
|
(SINT64) (saved_ca - T1 -
|
(T1 << 3)) * (SINT64) power10_table_128[scale_k -
|
(T1 << 3)) * (SINT64) power10_table_128[scale_k -
|
1].w[0];
|
1].w[0];
|
|
|
extra_digits++;
|
extra_digits++;
|
coefficient_b =
|
coefficient_b =
|
saved_cb + 100000000000000000ull +
|
saved_cb + 100000000000000000ull +
|
round_const_table[rmode][extra_digits];
|
round_const_table[rmode][extra_digits];
|
|
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_64x64_to_128 (CT, coefficient_b,
|
__mul_64x64_to_128 (CT, coefficient_b,
|
reciprocals10_64[extra_digits]);
|
reciprocals10_64[extra_digits]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
amount = short_recip_scale[extra_digits];
|
amount = short_recip_scale[extra_digits];
|
C0_64 = CT.w[1] >> amount;
|
C0_64 = CT.w[1] >> amount;
|
|
|
// result coefficient
|
// result coefficient
|
C64 = C0_64 + coefficient_a;
|
C64 = C0_64 + coefficient_a;
|
} else if (C64 <= 1000000000000000ull) {
|
} else if (C64 <= 1000000000000000ull) {
|
// less than 16 digits in result
|
// less than 16 digits in result
|
coefficient_a =
|
coefficient_a =
|
(SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
|
(SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
|
1].w[0];
|
1].w[0];
|
//extra_digits --;
|
//extra_digits --;
|
exponent_b--;
|
exponent_b--;
|
coefficient_b =
|
coefficient_b =
|
(saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
|
(saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
|
round_const_table[rmode][extra_digits];
|
round_const_table[rmode][extra_digits];
|
|
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_64x64_to_128 (CT_new, coefficient_b,
|
__mul_64x64_to_128 (CT_new, coefficient_b,
|
reciprocals10_64[extra_digits]);
|
reciprocals10_64[extra_digits]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
|
amount = short_recip_scale[extra_digits];
|
amount = short_recip_scale[extra_digits];
|
C0_64 = CT_new.w[1] >> amount;
|
C0_64 = CT_new.w[1] >> amount;
|
|
|
// result coefficient
|
// result coefficient
|
C64_new = C0_64 + coefficient_a;
|
C64_new = C0_64 + coefficient_a;
|
if (C64_new < 10000000000000000ull) {
|
if (C64_new < 10000000000000000ull) {
|
C64 = C64_new;
|
C64 = C64_new;
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
CT = CT_new;
|
CT = CT_new;
|
#endif
|
#endif
|
} else
|
} else
|
exponent_b++;
|
exponent_b++;
|
}
|
}
|
|
|
}
|
}
|
|
|
}
|
}
|
|
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifndef IEEE_ROUND_NEAREST
|
#ifndef IEEE_ROUND_NEAREST
|
if (rmode == 0) //ROUNDING_TO_NEAREST
|
if (rmode == 0) //ROUNDING_TO_NEAREST
|
#endif
|
#endif
|
if (C64 & 1) {
|
if (C64 & 1) {
|
// check whether fractional part of initial_P/10^extra_digits is
|
// check whether fractional part of initial_P/10^extra_digits is
|
// exactly .5
|
// exactly .5
|
// this is the same as fractional part of
|
// this is the same as fractional part of
|
// (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
|
// (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
|
|
|
// get remainder
|
// get remainder
|
remainder_h = CT.w[1] << (64 - amount);
|
remainder_h = CT.w[1] << (64 - amount);
|
|
|
// test whether fractional part is 0
|
// test whether fractional part is 0
|
if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
|
if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
|
C64--;
|
C64--;
|
}
|
}
|
}
|
}
|
#endif
|
#endif
|
|
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
status = INEXACT_EXCEPTION;
|
status = INEXACT_EXCEPTION;
|
|
|
// get remainder
|
// get remainder
|
remainder_h = CT.w[1] << (64 - amount);
|
remainder_h = CT.w[1] << (64 - amount);
|
|
|
switch (rmode) {
|
switch (rmode) {
|
case ROUNDING_TO_NEAREST:
|
case ROUNDING_TO_NEAREST:
|
case ROUNDING_TIES_AWAY:
|
case ROUNDING_TIES_AWAY:
|
// test whether fractional part is 0
|
// test whether fractional part is 0
|
if ((remainder_h == 0x8000000000000000ull)
|
if ((remainder_h == 0x8000000000000000ull)
|
&& (CT.w[0] < reciprocals10_64[extra_digits]))
|
&& (CT.w[0] < reciprocals10_64[extra_digits]))
|
status = EXACT_STATUS;
|
status = EXACT_STATUS;
|
break;
|
break;
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
case ROUNDING_TO_ZERO:
|
case ROUNDING_TO_ZERO:
|
if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
|
if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
|
status = EXACT_STATUS;
|
status = EXACT_STATUS;
|
//if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
|
//if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
|
break;
|
break;
|
default:
|
default:
|
// round up
|
// round up
|
__add_carry_out (tmp, carry, CT.w[0],
|
__add_carry_out (tmp, carry, CT.w[0],
|
reciprocals10_64[extra_digits]);
|
reciprocals10_64[extra_digits]);
|
if ((remainder_h >> (64 - amount)) + carry >=
|
if ((remainder_h >> (64 - amount)) + carry >=
|
(((UINT64) 1) << amount))
|
(((UINT64) 1) << amount))
|
status = EXACT_STATUS;
|
status = EXACT_STATUS;
|
break;
|
break;
|
}
|
}
|
__set_status_flags (pfpsf, status);
|
__set_status_flags (pfpsf, status);
|
|
|
#endif
|
#endif
|
|
|
res =
|
res =
|
fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
|
fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
|
rnd_mode, pfpsf);
|
rnd_mode, pfpsf);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
|
|
