/* 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/>. */
|
|
|
#define BID_128RES
|
#define BID_128RES
|
#include "bid_div_macros.h"
|
#include "bid_div_macros.h"
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#include <fenv.h>
|
#include <fenv.h>
|
|
|
#define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
|
#define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
|
#endif
|
#endif
|
|
|
extern UINT32 convert_table[5][128][2];
|
extern UINT32 convert_table[5][128][2];
|
extern SINT8 factors[][2];
|
extern SINT8 factors[][2];
|
extern UINT8 packed_10000_zeros[];
|
extern UINT8 packed_10000_zeros[];
|
|
|
BID128_FUNCTION_ARG2 (bid128_div, x, y)
|
BID128_FUNCTION_ARG2 (bid128_div, x, y)
|
|
|
UINT256 CA4, CA4r, P256;
|
UINT256 CA4, CA4r, P256;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
valid_y;
|
valid_y;
|
int_float fx, fy, f64;
|
int_float fx, fy, f64;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
digits_q, amount;
|
digits_q, amount;
|
int nzeros, i, j, k, d5;
|
int nzeros, i, j, k, d5;
|
unsigned rmode;
|
unsigned rmode;
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
fexcept_t binaryflags = 0;
|
fexcept_t binaryflags = 0;
|
#endif
|
#endif
|
|
|
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
|
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
|
|
|
// unpack arguments, check for NaN or Infinity
|
// unpack arguments, check for NaN or Infinity
|
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
|
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
|
// test if x is NaN
|
// test if x is NaN
|
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
|
if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
|
(y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
|
(y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] = (CX.w[1]) & QUIET_MASK64;
|
res.w[1] = (CX.w[1]) & QUIET_MASK64;
|
res.w[0] = CX.w[0];
|
res.w[0] = CX.w[0];
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// x is Infinity?
|
// x is Infinity?
|
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// check if y is Inf.
|
// check if y is Inf.
|
if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
// 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.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is NaN?
|
// y is NaN?
|
if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
|
if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
|
// return NaN
|
// return NaN
|
{
|
{
|
// return +/-Inf
|
// return +/-Inf
|
res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
|
res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
|
0x7800000000000000ull;
|
0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// x is 0
|
// x is 0
|
if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
|
if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
|
if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
|
if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
// x=y=0, return NaN
|
// x=y=0, return NaN
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// return 0
|
// return 0
|
res.w[1] = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
|
res.w[1] = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
else if (exponent_x < 0)
|
else if (exponent_x < 0)
|
exponent_x = 0;
|
exponent_x = 0;
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
if (!valid_y) {
|
if (!valid_y) {
|
// y is Inf. or NaN
|
// y is Inf. or NaN
|
|
|
// test if y is NaN
|
// test if y is NaN
|
if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
|
if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] = CY.w[1] & QUIET_MASK64;
|
res.w[1] = CY.w[1] & QUIET_MASK64;
|
res.w[0] = CY.w[0];
|
res.w[0] = CY.w[0];
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is Infinity?
|
// y is Infinity?
|
if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// return +/-0
|
// return +/-0
|
res.w[1] = sign_x ^ sign_y;
|
res.w[1] = sign_x ^ sign_y;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is 0, return +/-Inf
|
// y is 0, return +/-Inf
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] =
|
res.w[1] =
|
((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
|
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
// CX < CY
|
// CX < CY
|
|
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
|
|
// fx ~ CX, fy ~ CY
|
// fx ~ CX, fy ~ CY
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
// expon_cy - expon_cx
|
// expon_cy - expon_cx
|
bin_index = (fy.i - fx.i) >> 23;
|
bin_index = (fy.i - fx.i) >> 23;
|
|
|
if (CX.w[1]) {
|
if (CX.w[1]) {
|
T = power10_index_binexp_128[bin_index].w[0];
|
T = power10_index_binexp_128[bin_index].w[0];
|
__mul_64x128_short (CA, T, CX);
|
__mul_64x128_short (CA, T, CX);
|
} else {
|
} else {
|
T128 = power10_index_binexp_128[bin_index];
|
T128 = power10_index_binexp_128[bin_index];
|
__mul_64x128_short (CA, CX.w[0], T128);
|
__mul_64x128_short (CA, CX.w[0], T128);
|
}
|
}
|
|
|
ed2 = 33;
|
ed2 = 33;
|
if (__unsigned_compare_gt_128 (CY, CA))
|
if (__unsigned_compare_gt_128 (CY, CA))
|
ed2++;
|
ed2++;
|
|
|
T128 = power10_table_128[ed2];
|
T128 = power10_table_128[ed2];
|
__mul_128x128_to_256 (CA4, CA, T128);
|
__mul_128x128_to_256 (CA4, CA, T128);
|
|
|
ed2 += estimate_decimal_digits[bin_index];
|
ed2 += estimate_decimal_digits[bin_index];
|
CQ.w[0] = CQ.w[1] = 0;
|
CQ.w[0] = CQ.w[1] = 0;
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
|
|
} else {
|
} else {
|
// get CQ = CX/CY
|
// get CQ = CX/CY
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
|
|
if (!CR.w[1] && !CR.w[0]) {
|
if (!CR.w[1] && !CR.w[0]) {
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
pfpsf);
|
pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// get number of decimal digits in CQ
|
// get number of decimal digits in CQ
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
// binary expon. of CQ
|
// binary expon. of CQ
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
|
|
digits_q = estimate_decimal_digits[bin_expon];
|
digits_q = estimate_decimal_digits[bin_expon];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
digits_q++;
|
digits_q++;
|
|
|
ed2 = 34 - digits_q;
|
ed2 = 34 - digits_q;
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
__mul_128x128_to_256 (CA4, CR, T128);
|
__mul_128x128_to_256 (CA4, CR, T128);
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
__mul_128x128_low (CQ, CQ, T128);
|
__mul_128x128_low (CQ, CQ, T128);
|
|
|
}
|
}
|
|
|
__div_256_by_128 (&CQ, &CA4, CY);
|
__div_256_by_128 (&CQ, &CA4, CY);
|
|
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
else
|
else
|
#endif
|
#endif
|
#else
|
#else
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
if (!CA4.w[0] && !CA4.w[1])
|
if (!CA4.w[0] && !CA4.w[1])
|
#endif
|
#endif
|
#endif
|
#endif
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
// check whether result is exact
|
// check whether result is exact
|
{
|
{
|
// check whether CX, CY are short
|
// check whether CX, CY are short
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
i = (int) CY.w[0] - 1;
|
i = (int) CY.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
// difference in powers of 2 factors for Y and X
|
// difference in powers of 2 factors for Y and X
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
// difference in powers of 5 factors
|
// difference in powers of 5 factors
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
if (d5 < nzeros)
|
if (d5 < nzeros)
|
nzeros = d5;
|
nzeros = d5;
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128_long (CQ, Qh, amount);
|
__shr_128_long (CQ, Qh, amount);
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
// decompose Q as Qh*10^17 + Ql
|
// decompose Q as Qh*10^17 + Ql
|
//T128 = reciprocals10_128[17];
|
//T128 = reciprocals10_128[17];
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
__mul_128x128_to_256 (P256, CQ, T128);
|
__mul_128x128_to_256 (P256, CQ, T128);
|
//amount = recip_scale[17];
|
//amount = recip_scale[17];
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
|
|
if (!Q_low) {
|
if (!Q_low) {
|
diff_expon += 17;
|
diff_expon += 17;
|
|
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_high >> 26;
|
QX = Q_high >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
amount = short_recip_scale[nzeros];
|
amount = short_recip_scale[nzeros];
|
CQ.w[0] = CQ.w[1] >> amount;
|
CQ.w[0] = CQ.w[1] >> amount;
|
} else
|
} else
|
CQ.w[0] = Q_high;
|
CQ.w[0] = Q_high;
|
CQ.w[1] = 0;
|
CQ.w[1] = 0;
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_low >> 26;
|
QX = Q_low >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
|
|
//now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
//now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128 (CQ, Qh, amount);
|
__shr_128 (CQ, Qh, amount);
|
}
|
}
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
|
|
}
|
}
|
}
|
}
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#endif
|
#endif
|
|
|
if (diff_expon >= 0) {
|
if (diff_expon >= 0) {
|
#ifdef IEEE_ROUND_NEAREST
|
#ifdef IEEE_ROUND_NEAREST
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
switch (rmode) {
|
switch (rmode) {
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_TIES_AWAY:
|
case ROUNDING_TIES_AWAY:
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
case ROUNDING_TO_ZERO:
|
case ROUNDING_TO_ZERO:
|
break;
|
break;
|
default: // rounding up
|
default: // rounding up
|
CQ.w[0]++;
|
CQ.w[0]++;
|
if (!CQ.w[0])
|
if (!CQ.w[0])
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
}
|
}
|
#endif
|
#endif
|
#endif
|
#endif
|
|
|
} else {
|
} else {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#endif
|
#endif
|
|
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
|
|
}
|
}
|
|
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
|
|
|
|
//#define LEAVE_TRAILING_ZEROS
|
//#define LEAVE_TRAILING_ZEROS
|
|
|
TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128, bid128dd_div, UINT64, x,
|
TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128, bid128dd_div, UINT64, x,
|
UINT64, y)
|
UINT64, y)
|
|
|
UINT256 CA4, CA4r, P256;
|
UINT256 CA4, CA4r, P256;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
valid_y;
|
valid_y;
|
int_float fx, fy, f64;
|
int_float fx, fy, f64;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
digits_q, amount;
|
digits_q, amount;
|
int nzeros, i, j, k, d5;
|
int nzeros, i, j, k, d5;
|
unsigned rmode;
|
unsigned rmode;
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
fexcept_t binaryflags = 0;
|
fexcept_t binaryflags = 0;
|
#endif
|
#endif
|
|
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
|
|
|
// unpack arguments, check for NaN or Infinity
|
// unpack arguments, check for NaN or Infinity
|
CX.w[1] = 0;
|
CX.w[1] = 0;
|
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
|
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
|
if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
|
|
// 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
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
|
res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
|
res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// x is Infinity?
|
// x is Infinity?
|
if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// check if y is Inf.
|
// check if y is Inf.
|
if ((((y) & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
if ((((y) & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
// 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.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
if ((((y) & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
|
if ((((y) & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
|
// otherwise return +/-Inf
|
// otherwise return +/-Inf
|
res.w[1] =
|
res.w[1] =
|
(((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
(((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// x is 0
|
// x is 0
|
if ((((y) & 0x7800000000000000ull) != 0x7800000000000000ull)) {
|
if ((((y) & 0x7800000000000000ull) != 0x7800000000000000ull)) {
|
if(!CY.w[0]) {
|
if(!CY.w[0]) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
// x=y=0, return NaN
|
// x=y=0, return NaN
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// return 0
|
// return 0
|
res.w[1] = ((x) ^ (y)) & 0x8000000000000000ull;
|
res.w[1] = ((x) ^ (y)) & 0x8000000000000000ull;
|
if (((y) & 0x6000000000000000ull) == 0x6000000000000000ull)
|
if (((y) & 0x6000000000000000ull) == 0x6000000000000000ull)
|
exponent_y = ((UINT32) ((y) >> 51)) & 0x3ff;
|
exponent_y = ((UINT32) ((y) >> 51)) & 0x3ff;
|
else
|
else
|
exponent_y = ((UINT32) ((y) >> 53)) & 0x3ff;
|
exponent_y = ((UINT32) ((y) >> 53)) & 0x3ff;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
else if (exponent_x < 0)
|
else if (exponent_x < 0)
|
exponent_x = 0;
|
exponent_x = 0;
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
|
|
CY.w[1] = 0;
|
CY.w[1] = 0;
|
if (!valid_y) {
|
if (!valid_y) {
|
// y is Inf. or NaN
|
// y is Inf. or NaN
|
|
|
// test if y is NaN
|
// test if y is NaN
|
if ((y & NAN_MASK64) == NAN_MASK64) {
|
if ((y & NAN_MASK64) == NAN_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.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
|
res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
|
res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is Infinity?
|
// y is Infinity?
|
if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// return +/-0
|
// return +/-0
|
res.w[1] = sign_x ^ sign_y;
|
res.w[1] = sign_x ^ sign_y;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is 0, return +/-Inf
|
// y is 0, return +/-Inf
|
res.w[1] =
|
res.w[1] =
|
(((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
(((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
|
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
// CX < CY
|
// CX < CY
|
|
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
|
|
// fx ~ CX, fy ~ CY
|
// fx ~ CX, fy ~ CY
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
// expon_cy - expon_cx
|
// expon_cy - expon_cx
|
bin_index = (fy.i - fx.i) >> 23;
|
bin_index = (fy.i - fx.i) >> 23;
|
|
|
if (CX.w[1]) {
|
if (CX.w[1]) {
|
T = power10_index_binexp_128[bin_index].w[0];
|
T = power10_index_binexp_128[bin_index].w[0];
|
__mul_64x128_short (CA, T, CX);
|
__mul_64x128_short (CA, T, CX);
|
} else {
|
} else {
|
T128 = power10_index_binexp_128[bin_index];
|
T128 = power10_index_binexp_128[bin_index];
|
__mul_64x128_short (CA, CX.w[0], T128);
|
__mul_64x128_short (CA, CX.w[0], T128);
|
}
|
}
|
|
|
ed2 = 33;
|
ed2 = 33;
|
if (__unsigned_compare_gt_128 (CY, CA))
|
if (__unsigned_compare_gt_128 (CY, CA))
|
ed2++;
|
ed2++;
|
|
|
T128 = power10_table_128[ed2];
|
T128 = power10_table_128[ed2];
|
__mul_128x128_to_256 (CA4, CA, T128);
|
__mul_128x128_to_256 (CA4, CA, T128);
|
|
|
ed2 += estimate_decimal_digits[bin_index];
|
ed2 += estimate_decimal_digits[bin_index];
|
CQ.w[0] = CQ.w[1] = 0;
|
CQ.w[0] = CQ.w[1] = 0;
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
|
|
} else {
|
} else {
|
// get CQ = CX/CY
|
// get CQ = CX/CY
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
|
|
if (!CR.w[1] && !CR.w[0]) {
|
if (!CR.w[1] && !CR.w[0]) {
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
pfpsf);
|
pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// get number of decimal digits in CQ
|
// get number of decimal digits in CQ
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
// binary expon. of CQ
|
// binary expon. of CQ
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
|
|
digits_q = estimate_decimal_digits[bin_expon];
|
digits_q = estimate_decimal_digits[bin_expon];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
digits_q++;
|
digits_q++;
|
|
|
ed2 = 34 - digits_q;
|
ed2 = 34 - digits_q;
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
__mul_128x128_to_256 (CA4, CR, T128);
|
__mul_128x128_to_256 (CA4, CR, T128);
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
__mul_128x128_low (CQ, CQ, T128);
|
__mul_128x128_low (CQ, CQ, T128);
|
|
|
}
|
}
|
|
|
__div_256_by_128 (&CQ, &CA4, CY);
|
__div_256_by_128 (&CQ, &CA4, CY);
|
|
|
|
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
else
|
else
|
#endif
|
#endif
|
#else
|
#else
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
if (!CA4.w[0] && !CA4.w[1])
|
if (!CA4.w[0] && !CA4.w[1])
|
#endif
|
#endif
|
#endif
|
#endif
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
// check whether result is exact
|
// check whether result is exact
|
{
|
{
|
// check whether CX, CY are short
|
// check whether CX, CY are short
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
i = (int) CY.w[0] - 1;
|
i = (int) CY.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
// difference in powers of 2 factors for Y and X
|
// difference in powers of 2 factors for Y and X
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
// difference in powers of 5 factors
|
// difference in powers of 5 factors
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
if (d5 < nzeros)
|
if (d5 < nzeros)
|
nzeros = d5;
|
nzeros = d5;
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128_long (CQ, Qh, amount);
|
__shr_128_long (CQ, Qh, amount);
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
// decompose Q as Qh*10^17 + Ql
|
// decompose Q as Qh*10^17 + Ql
|
//T128 = reciprocals10_128[17];
|
//T128 = reciprocals10_128[17];
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
__mul_128x128_to_256 (P256, CQ, T128);
|
__mul_128x128_to_256 (P256, CQ, T128);
|
//amount = recip_scale[17];
|
//amount = recip_scale[17];
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
|
|
if (!Q_low) {
|
if (!Q_low) {
|
diff_expon += 17;
|
diff_expon += 17;
|
|
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_high >> 26;
|
QX = Q_high >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
amount = short_recip_scale[nzeros];
|
amount = short_recip_scale[nzeros];
|
CQ.w[0] = CQ.w[1] >> amount;
|
CQ.w[0] = CQ.w[1] >> amount;
|
} else
|
} else
|
CQ.w[0] = Q_high;
|
CQ.w[0] = Q_high;
|
CQ.w[1] = 0;
|
CQ.w[1] = 0;
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_low >> 26;
|
QX = Q_low >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128 (CQ, Qh, amount);
|
__shr_128 (CQ, Qh, amount);
|
}
|
}
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
|
|
}
|
}
|
}
|
}
|
get_BID128(&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
|
get_BID128(&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#endif
|
#endif
|
|
|
if (diff_expon >= 0) {
|
if (diff_expon >= 0) {
|
#ifdef IEEE_ROUND_NEAREST
|
#ifdef IEEE_ROUND_NEAREST
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
switch (rmode) {
|
switch (rmode) {
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_TIES_AWAY:
|
case ROUNDING_TIES_AWAY:
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
case ROUNDING_TO_ZERO:
|
case ROUNDING_TO_ZERO:
|
break;
|
break;
|
default: // rounding up
|
default: // rounding up
|
CQ.w[0]++;
|
CQ.w[0]++;
|
if (!CQ.w[0])
|
if (!CQ.w[0])
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
}
|
}
|
#endif
|
#endif
|
#endif
|
#endif
|
|
|
} else {
|
} else {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#endif
|
#endif
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
|
|
}
|
}
|
|
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
|
|
|
|
BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div, UINT64, x, y)
|
BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div, UINT64, x, y)
|
UINT256 CA4, CA4r, P256;
|
UINT256 CA4, CA4r, P256;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, valid_y,
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, valid_y,
|
PD;
|
PD;
|
int_float fx, fy, f64;
|
int_float fx, fy, f64;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
digits_q, amount;
|
digits_q, amount;
|
int nzeros, i, j, k, d5;
|
int nzeros, i, j, k, d5;
|
unsigned rmode;
|
unsigned rmode;
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
fexcept_t binaryflags = 0;
|
fexcept_t binaryflags = 0;
|
#endif
|
#endif
|
|
|
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
|
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
|
|
|
// unpack arguments, check for NaN or Infinity
|
// unpack arguments, check for NaN or Infinity
|
CX.w[1] = 0;
|
CX.w[1] = 0;
|
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
|
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
|
if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
|
|
// 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
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
|
res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
|
res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// x is Infinity?
|
// x is Infinity?
|
if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// check if y is Inf.
|
// check if y is Inf.
|
if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
// 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.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
|
if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
|
// otherwise return +/-Inf
|
// otherwise return +/-Inf
|
res.w[1] =
|
res.w[1] =
|
((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// x is 0
|
// x is 0
|
if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
|
if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
|
if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
|
if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
// x=y=0, return NaN
|
// x=y=0, return NaN
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// return 0
|
// return 0
|
res.w[1] = (x ^ y.w[1]) & 0x8000000000000000ull;
|
res.w[1] = (x ^ y.w[1]) & 0x8000000000000000ull;
|
exponent_x = exponent_x - exponent_y + (DECIMAL_EXPONENT_BIAS_128<<1) - DECIMAL_EXPONENT_BIAS;
|
exponent_x = exponent_x - exponent_y + (DECIMAL_EXPONENT_BIAS_128<<1) - DECIMAL_EXPONENT_BIAS;
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
else if (exponent_x < 0)
|
else if (exponent_x < 0)
|
exponent_x = 0;
|
exponent_x = 0;
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS);
|
exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS);
|
|
|
if (!valid_y) {
|
if (!valid_y) {
|
// y is Inf. or NaN
|
// y is Inf. or NaN
|
|
|
// test if y is NaN
|
// test if y is NaN
|
if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
|
if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] = CY.w[1] & QUIET_MASK64;
|
res.w[1] = CY.w[1] & QUIET_MASK64;
|
res.w[0] = CY.w[0];
|
res.w[0] = CY.w[0];
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is Infinity?
|
// y is Infinity?
|
if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// return +/-0
|
// return +/-0
|
res.w[1] = sign_x ^ sign_y;
|
res.w[1] = sign_x ^ sign_y;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is 0, return +/-Inf
|
// y is 0, return +/-Inf
|
res.w[1] =
|
res.w[1] =
|
((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
|
|
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
// CX < CY
|
// CX < CY
|
|
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
|
|
// fx ~ CX, fy ~ CY
|
// fx ~ CX, fy ~ CY
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
// expon_cy - expon_cx
|
// expon_cy - expon_cx
|
bin_index = (fy.i - fx.i) >> 23;
|
bin_index = (fy.i - fx.i) >> 23;
|
|
|
if (CX.w[1]) {
|
if (CX.w[1]) {
|
T = power10_index_binexp_128[bin_index].w[0];
|
T = power10_index_binexp_128[bin_index].w[0];
|
__mul_64x128_short (CA, T, CX);
|
__mul_64x128_short (CA, T, CX);
|
} else {
|
} else {
|
T128 = power10_index_binexp_128[bin_index];
|
T128 = power10_index_binexp_128[bin_index];
|
__mul_64x128_short (CA, CX.w[0], T128);
|
__mul_64x128_short (CA, CX.w[0], T128);
|
}
|
}
|
|
|
ed2 = 33;
|
ed2 = 33;
|
if (__unsigned_compare_gt_128 (CY, CA))
|
if (__unsigned_compare_gt_128 (CY, CA))
|
ed2++;
|
ed2++;
|
|
|
T128 = power10_table_128[ed2];
|
T128 = power10_table_128[ed2];
|
__mul_128x128_to_256 (CA4, CA, T128);
|
__mul_128x128_to_256 (CA4, CA, T128);
|
|
|
ed2 += estimate_decimal_digits[bin_index];
|
ed2 += estimate_decimal_digits[bin_index];
|
CQ.w[0] = CQ.w[1] = 0;
|
CQ.w[0] = CQ.w[1] = 0;
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
|
|
} else {
|
} else {
|
// get CQ = CX/CY
|
// get CQ = CX/CY
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
|
|
if (!CR.w[1] && !CR.w[0]) {
|
if (!CR.w[1] && !CR.w[0]) {
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
pfpsf);
|
pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// get number of decimal digits in CQ
|
// get number of decimal digits in CQ
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
// binary expon. of CQ
|
// binary expon. of CQ
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
|
|
digits_q = estimate_decimal_digits[bin_expon];
|
digits_q = estimate_decimal_digits[bin_expon];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
digits_q++;
|
digits_q++;
|
|
|
ed2 = 34 - digits_q;
|
ed2 = 34 - digits_q;
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
__mul_128x128_to_256 (CA4, CR, T128);
|
__mul_128x128_to_256 (CA4, CR, T128);
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
__mul_128x128_low (CQ, CQ, T128);
|
__mul_128x128_low (CQ, CQ, T128);
|
|
|
}
|
}
|
|
|
__div_256_by_128 (&CQ, &CA4, CY);
|
__div_256_by_128 (&CQ, &CA4, CY);
|
|
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
else
|
else
|
#endif
|
#endif
|
#else
|
#else
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
if (!CA4.w[0] && !CA4.w[1])
|
if (!CA4.w[0] && !CA4.w[1])
|
#endif
|
#endif
|
#endif
|
#endif
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
// check whether result is exact
|
// check whether result is exact
|
{
|
{
|
//printf("ed2=%d,nz=%d,a=%d,CQ="LX16","LX16", RH="LX16", RL="LX16"\n",ed2,nzeros,amount,CQ.w[1],CQ.w[0],reciprocals10_128[nzeros].w[1],reciprocals10_128[nzeros].w[0]);fflush(stdout);
|
//printf("ed2=%d,nz=%d,a=%d,CQ="LX16","LX16", RH="LX16", RL="LX16"\n",ed2,nzeros,amount,CQ.w[1],CQ.w[0],reciprocals10_128[nzeros].w[1],reciprocals10_128[nzeros].w[0]);fflush(stdout);
|
// check whether CX, CY are short
|
// check whether CX, CY are short
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
i = (int) CY.w[0] - 1;
|
i = (int) CY.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
// difference in powers of 2 factors for Y and X
|
// difference in powers of 2 factors for Y and X
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
// difference in powers of 5 factors
|
// difference in powers of 5 factors
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
if (d5 < nzeros)
|
if (d5 < nzeros)
|
nzeros = d5;
|
nzeros = d5;
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128_long (CQ, Qh, amount);
|
__shr_128_long (CQ, Qh, amount);
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
// decompose Q as Qh*10^17 + Ql
|
// decompose Q as Qh*10^17 + Ql
|
//T128 = reciprocals10_128[17];
|
//T128 = reciprocals10_128[17];
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
__mul_128x128_to_256 (P256, CQ, T128);
|
__mul_128x128_to_256 (P256, CQ, T128);
|
//amount = recip_scale[17];
|
//amount = recip_scale[17];
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
|
|
if (!Q_low) {
|
if (!Q_low) {
|
diff_expon += 17;
|
diff_expon += 17;
|
|
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_high >> 26;
|
QX = Q_high >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
//printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
|
//printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
amount = short_recip_scale[nzeros];
|
amount = short_recip_scale[nzeros];
|
CQ.w[0] = CQ.w[1] >> amount;
|
CQ.w[0] = CQ.w[1] >> amount;
|
} else
|
} else
|
CQ.w[0] = Q_high;
|
CQ.w[0] = Q_high;
|
CQ.w[1] = 0;
|
CQ.w[1] = 0;
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_low >> 26;
|
QX = Q_low >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
//printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
|
//printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128 (CQ, Qh, amount);
|
__shr_128 (CQ, Qh, amount);
|
}
|
}
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
|
|
}
|
}
|
}
|
}
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
pfpsf);
|
pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#endif
|
#endif
|
|
|
if (diff_expon >= 0) {
|
if (diff_expon >= 0) {
|
#ifdef IEEE_ROUND_NEAREST
|
#ifdef IEEE_ROUND_NEAREST
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
switch (rmode) {
|
switch (rmode) {
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_TIES_AWAY:
|
case ROUNDING_TIES_AWAY:
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
case ROUNDING_TO_ZERO:
|
case ROUNDING_TO_ZERO:
|
break;
|
break;
|
default: // rounding up
|
default: // rounding up
|
CQ.w[0]++;
|
CQ.w[0]++;
|
if (!CQ.w[0])
|
if (!CQ.w[0])
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
}
|
}
|
#endif
|
#endif
|
#endif
|
#endif
|
|
|
} else {
|
} else {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#endif
|
#endif
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
|
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
|
|
}
|
}
|
|
|
|
|
BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div, x, UINT64, y)
|
BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div, x, UINT64, y)
|
UINT256 CA4, CA4r, P256;
|
UINT256 CA4, CA4r, P256;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
|
valid_y;
|
valid_y;
|
int_float fx, fy, f64;
|
int_float fx, fy, f64;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
|
digits_q, amount;
|
digits_q, amount;
|
int nzeros, i, j, k, d5, rmode;
|
int nzeros, i, j, k, d5, rmode;
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
fexcept_t binaryflags = 0;
|
fexcept_t binaryflags = 0;
|
#endif
|
#endif
|
|
|
|
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
|
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
|
// unpack arguments, check for NaN or Infinity
|
// unpack arguments, check for NaN or Infinity
|
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
|
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
|
// test if x is NaN
|
// test if x is NaN
|
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
|
if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
|
(y & 0x7e00000000000000ull) == 0x7e00000000000000ull)
|
(y & 0x7e00000000000000ull) == 0x7e00000000000000ull)
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] = (CX.w[1]) & QUIET_MASK64;
|
res.w[1] = (CX.w[1]) & QUIET_MASK64;
|
res.w[0] = CX.w[0];
|
res.w[0] = CX.w[0];
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// x is Infinity?
|
// x is Infinity?
|
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
|
// check if y is Inf.
|
// check if y is Inf.
|
if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull))
|
// 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.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is NaN?
|
// y is NaN?
|
if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull))
|
if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull))
|
// return NaN
|
// return NaN
|
{
|
{
|
// return +/-Inf
|
// return +/-Inf
|
res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull) |
|
res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull) |
|
0x7800000000000000ull;
|
0x7800000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
// x is 0
|
// x is 0
|
if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
|
if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
|
if (!CY.w[0]) {
|
if (!CY.w[0]) {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
__set_status_flags (pfpsf, INVALID_EXCEPTION);
|
#endif
|
#endif
|
// x=y=0, return NaN
|
// x=y=0, return NaN
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[1] = 0x7c00000000000000ull;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// return 0
|
// return 0
|
res.w[1] = (x.w[1] ^ y) & 0x8000000000000000ull;
|
res.w[1] = (x.w[1] ^ y) & 0x8000000000000000ull;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
|
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
if (exponent_x > DECIMAL_MAX_EXPON_128)
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
exponent_x = DECIMAL_MAX_EXPON_128;
|
else if (exponent_x < 0)
|
else if (exponent_x < 0)
|
exponent_x = 0;
|
exponent_x = 0;
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[1] |= (((UINT64) exponent_x) << 49);
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
}
|
}
|
CY.w[1] = 0;
|
CY.w[1] = 0;
|
if (!valid_y) {
|
if (!valid_y) {
|
// y is Inf. or NaN
|
// y is Inf. or NaN
|
|
|
// test if y is NaN
|
// test if y is NaN
|
if ((y & NAN_MASK64) == NAN_MASK64) {
|
if ((y & NAN_MASK64) == NAN_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.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
|
res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
|
res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
|
res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is Infinity?
|
// y is Infinity?
|
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
|
if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
|
// return +/-0
|
// return +/-0
|
res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull);
|
res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull);
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// y is 0
|
// y is 0
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
__set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
|
#endif
|
#endif
|
res.w[1] = (sign_x ^ sign_y) | INFINITY_MASK64;
|
res.w[1] = (sign_x ^ sign_y) | INFINITY_MASK64;
|
res.w[0] = 0;
|
res.w[0] = 0;
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
|
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
|
|
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
if (__unsigned_compare_gt_128 (CY, CX)) {
|
// CX < CY
|
// CX < CY
|
|
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
|
|
// fx ~ CX, fy ~ CY
|
// fx ~ CX, fy ~ CY
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
|
// expon_cy - expon_cx
|
// expon_cy - expon_cx
|
bin_index = (fy.i - fx.i) >> 23;
|
bin_index = (fy.i - fx.i) >> 23;
|
|
|
if (CX.w[1]) {
|
if (CX.w[1]) {
|
T = power10_index_binexp_128[bin_index].w[0];
|
T = power10_index_binexp_128[bin_index].w[0];
|
__mul_64x128_short (CA, T, CX);
|
__mul_64x128_short (CA, T, CX);
|
} else {
|
} else {
|
T128 = power10_index_binexp_128[bin_index];
|
T128 = power10_index_binexp_128[bin_index];
|
__mul_64x128_short (CA, CX.w[0], T128);
|
__mul_64x128_short (CA, CX.w[0], T128);
|
}
|
}
|
|
|
ed2 = 33;
|
ed2 = 33;
|
if (__unsigned_compare_gt_128 (CY, CA))
|
if (__unsigned_compare_gt_128 (CY, CA))
|
ed2++;
|
ed2++;
|
|
|
T128 = power10_table_128[ed2];
|
T128 = power10_table_128[ed2];
|
__mul_128x128_to_256 (CA4, CA, T128);
|
__mul_128x128_to_256 (CA4, CA, T128);
|
|
|
ed2 += estimate_decimal_digits[bin_index];
|
ed2 += estimate_decimal_digits[bin_index];
|
CQ.w[0] = CQ.w[1] = 0;
|
CQ.w[0] = CQ.w[1] = 0;
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
|
|
} else {
|
} else {
|
// get CQ = CX/CY
|
// get CQ = CX/CY
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
__div_128_by_128 (&CQ, &CR, CX, CY);
|
|
|
if (!CR.w[1] && !CR.w[0]) {
|
if (!CR.w[1] && !CR.w[0]) {
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
|
pfpsf);
|
pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
// get number of decimal digits in CQ
|
// get number of decimal digits in CQ
|
// 2^64
|
// 2^64
|
f64.i = 0x5f800000;
|
f64.i = 0x5f800000;
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
|
// binary expon. of CQ
|
// binary expon. of CQ
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
bin_expon = (fx.i - 0x3f800000) >> 23;
|
|
|
digits_q = estimate_decimal_digits[bin_expon];
|
digits_q = estimate_decimal_digits[bin_expon];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
if (__unsigned_compare_ge_128 (CQ, TP128))
|
digits_q++;
|
digits_q++;
|
|
|
ed2 = 34 - digits_q;
|
ed2 = 34 - digits_q;
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[0] = power10_table_128[ed2].w[0];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
T128.w[1] = power10_table_128[ed2].w[1];
|
__mul_128x128_to_256 (CA4, CR, T128);
|
__mul_128x128_to_256 (CA4, CR, T128);
|
diff_expon = diff_expon - ed2;
|
diff_expon = diff_expon - ed2;
|
__mul_128x128_low (CQ, CQ, T128);
|
__mul_128x128_low (CQ, CQ, T128);
|
|
|
}
|
}
|
|
|
__div_256_by_128 (&CQ, &CA4, CY);
|
__div_256_by_128 (&CQ, &CA4, CY);
|
|
|
|
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
else
|
else
|
#endif
|
#endif
|
#else
|
#else
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
if (!CA4.w[0] && !CA4.w[1])
|
if (!CA4.w[0] && !CA4.w[1])
|
#endif
|
#endif
|
#endif
|
#endif
|
#ifndef LEAVE_TRAILING_ZEROS
|
#ifndef LEAVE_TRAILING_ZEROS
|
// check whether result is exact
|
// check whether result is exact
|
{
|
{
|
// check whether CX, CY are short
|
// check whether CX, CY are short
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
|
i = (int) CY.w[0] - 1;
|
i = (int) CY.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
j = (int) CX.w[0] - 1;
|
// difference in powers of 2 factors for Y and X
|
// difference in powers of 2 factors for Y and X
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
nzeros = ed2 - factors[i][0] + factors[j][0];
|
// difference in powers of 5 factors
|
// difference in powers of 5 factors
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
d5 = ed2 - factors[i][1] + factors[j][1];
|
if (d5 < nzeros)
|
if (d5 < nzeros)
|
nzeros = d5;
|
nzeros = d5;
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128_long (CQ, Qh, amount);
|
__shr_128_long (CQ, Qh, amount);
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
// decompose Q as Qh*10^17 + Ql
|
// decompose Q as Qh*10^17 + Ql
|
//T128 = reciprocals10_128[17];
|
//T128 = reciprocals10_128[17];
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[0] = 0x44909befeb9fad49ull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
T128.w[1] = 0x000b877aa3236a4bull;
|
__mul_128x128_to_256 (P256, CQ, T128);
|
__mul_128x128_to_256 (P256, CQ, T128);
|
//amount = recip_scale[17];
|
//amount = recip_scale[17];
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
|
|
|
if (!Q_low) {
|
if (!Q_low) {
|
diff_expon += 17;
|
diff_expon += 17;
|
|
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[0] = Q_high & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_high >> 26;
|
QX = Q_high >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
__mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
|
|
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
// now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
|
amount = short_recip_scale[nzeros];
|
amount = short_recip_scale[nzeros];
|
CQ.w[0] = CQ.w[1] >> amount;
|
CQ.w[0] = CQ.w[1] >> amount;
|
} else
|
} else
|
CQ.w[0] = Q_high;
|
CQ.w[0] = Q_high;
|
CQ.w[1] = 0;
|
CQ.w[1] = 0;
|
|
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
} else {
|
} else {
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[0] = Q_low & 0x3ffffff;
|
tdigit[1] = 0;
|
tdigit[1] = 0;
|
QX = Q_low >> 26;
|
QX = Q_low >> 26;
|
QX32 = QX;
|
QX32 = QX;
|
nzeros = 0;
|
nzeros = 0;
|
|
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
for (j = 0; QX32; j++, QX32 >>= 7) {
|
k = (QX32 & 127);
|
k = (QX32 & 127);
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[0] += convert_table[j][k][0];
|
tdigit[1] += convert_table[j][k][1];
|
tdigit[1] += convert_table[j][k][1];
|
if (tdigit[0] >= 100000000) {
|
if (tdigit[0] >= 100000000) {
|
tdigit[0] -= 100000000;
|
tdigit[0] -= 100000000;
|
tdigit[1]++;
|
tdigit[1]++;
|
}
|
}
|
}
|
}
|
|
|
if (tdigit[1] >= 100000000) {
|
if (tdigit[1] >= 100000000) {
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
if (tdigit[1] >= 100000000)
|
if (tdigit[1] >= 100000000)
|
tdigit[1] -= 100000000;
|
tdigit[1] -= 100000000;
|
}
|
}
|
|
|
digit = tdigit[0];
|
digit = tdigit[0];
|
if (!digit && !tdigit[1])
|
if (!digit && !tdigit[1])
|
nzeros += 16;
|
nzeros += 16;
|
else {
|
else {
|
if (!digit) {
|
if (!digit) {
|
nzeros += 8;
|
nzeros += 8;
|
digit = tdigit[1];
|
digit = tdigit[1];
|
}
|
}
|
// decompose digit
|
// decompose digit
|
PD = (UINT64) digit *0x068DB8BBull;
|
PD = (UINT64) digit *0x068DB8BBull;
|
digit_h = (UINT32) (PD >> 40);
|
digit_h = (UINT32) (PD >> 40);
|
digit_low = digit - digit_h * 10000;
|
digit_low = digit - digit_h * 10000;
|
|
|
if (!digit_low)
|
if (!digit_low)
|
nzeros += 4;
|
nzeros += 4;
|
else
|
else
|
digit_h = digit_low;
|
digit_h = digit_low;
|
|
|
if (!(digit_h & 1))
|
if (!(digit_h & 1))
|
nzeros +=
|
nzeros +=
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
|
(digit_h & 7));
|
(digit_h & 7));
|
}
|
}
|
|
|
if (nzeros) {
|
if (nzeros) {
|
// get P*(2^M[extra_digits])/10^extra_digits
|
// get P*(2^M[extra_digits])/10^extra_digits
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
__mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
|
|
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
|
amount = recip_scale[nzeros];
|
amount = recip_scale[nzeros];
|
__shr_128 (CQ, Qh, amount);
|
__shr_128 (CQ, Qh, amount);
|
}
|
}
|
diff_expon += nzeros;
|
diff_expon += nzeros;
|
|
|
}
|
}
|
}
|
}
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
}
|
}
|
#endif
|
#endif
|
|
|
if (diff_expon >= 0) {
|
if (diff_expon >= 0) {
|
#ifdef IEEE_ROUND_NEAREST
|
#ifdef IEEE_ROUND_NEAREST
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
|
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
|
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
#else
|
#else
|
rmode = rnd_mode;
|
rmode = rnd_mode;
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
|
rmode = 3 - rmode;
|
rmode = 3 - rmode;
|
switch (rmode) {
|
switch (rmode) {
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
case ROUNDING_TO_NEAREST: // round to nearest code
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_TIES_AWAY:
|
case ROUNDING_TIES_AWAY:
|
// rounding
|
// rounding
|
// 2*CA4 - CY
|
// 2*CA4 - CY
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
CA4r.w[0] = CA4.w[0] + CA4.w[0];
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
__sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
|
CQ.w[0] += carry64;
|
CQ.w[0] += carry64;
|
if (CQ.w[0] < carry64)
|
if (CQ.w[0] < carry64)
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
case ROUNDING_DOWN:
|
case ROUNDING_DOWN:
|
case ROUNDING_TO_ZERO:
|
case ROUNDING_TO_ZERO:
|
break;
|
break;
|
default: // rounding up
|
default: // rounding up
|
CQ.w[0]++;
|
CQ.w[0]++;
|
if (!CQ.w[0])
|
if (!CQ.w[0])
|
CQ.w[1]++;
|
CQ.w[1]++;
|
break;
|
break;
|
}
|
}
|
#endif
|
#endif
|
#endif
|
#endif
|
|
|
} else {
|
} else {
|
#ifdef SET_STATUS_FLAGS
|
#ifdef SET_STATUS_FLAGS
|
if (CA4.w[0] || CA4.w[1]) {
|
if (CA4.w[0] || CA4.w[1]) {
|
// set status flags
|
// set status flags
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
|
}
|
}
|
#endif
|
#endif
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
|
|
}
|
}
|
|
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
|
#endif
|
#endif
|
BID_RETURN (res);
|
BID_RETURN (res);
|
|
|
}
|
}
|
|
|
