forked from KolibriOS/kolibrios
482 lines
15 KiB
C
482 lines
15 KiB
C
/* Copyright (C) 2007-2015 Free Software Foundation, Inc.
|
|
|
|
This file is part of GCC.
|
|
|
|
GCC is free software; you can redistribute it and/or modify it under
|
|
the terms of the GNU General Public License as published by the Free
|
|
Software Foundation; either version 3, or (at your option) any later
|
|
version.
|
|
|
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
Under Section 7 of GPL version 3, you are granted additional
|
|
permissions described in the GCC Runtime Library Exception, version
|
|
3.1, as published by the Free Software Foundation.
|
|
|
|
You should have received a copy of the GNU General Public License and
|
|
a copy of the GCC Runtime Library Exception along with this program;
|
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
#include "bid_internal.h"
|
|
|
|
/*****************************************************************************
|
|
* BID64 nextup
|
|
****************************************************************************/
|
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
void
|
|
bid64_nextup (UINT64 * pres,
|
|
UINT64 *
|
|
px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
UINT64 x = *px;
|
|
#else
|
|
UINT64
|
|
bid64_nextup (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
_EXC_INFO_PARAM) {
|
|
#endif
|
|
|
|
UINT64 res;
|
|
UINT64 x_sign;
|
|
UINT64 x_exp;
|
|
BID_UI64DOUBLE tmp1;
|
|
int x_nr_bits;
|
|
int q1, ind;
|
|
UINT64 C1; // C1 represents x_signif (UINT64)
|
|
|
|
// check for NaNs and infinities
|
|
if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
|
|
if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
|
|
x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
|
|
else
|
|
x = x & 0xfe03ffffffffffffull; // clear G6-G12
|
|
if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
|
|
// set invalid flag
|
|
*pfpsf |= INVALID_EXCEPTION;
|
|
// return quiet (SNaN)
|
|
res = x & 0xfdffffffffffffffull;
|
|
} else { // QNaN
|
|
res = x;
|
|
}
|
|
BID_RETURN (res);
|
|
} else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
|
|
if (!(x & 0x8000000000000000ull)) { // x is +inf
|
|
res = 0x7800000000000000ull;
|
|
} else { // x is -inf
|
|
res = 0xf7fb86f26fc0ffffull; // -MAXFP = -999...99 * 10^emax
|
|
}
|
|
BID_RETURN (res);
|
|
}
|
|
// unpack the argument
|
|
x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
|
// if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
|
|
if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
|
|
x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
|
|
C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
|
|
if (C1 > 9999999999999999ull) { // non-canonical
|
|
x_exp = 0;
|
|
C1 = 0;
|
|
}
|
|
} else {
|
|
x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
|
|
C1 = x & MASK_BINARY_SIG1;
|
|
}
|
|
|
|
// check for zeros (possibly from non-canonical values)
|
|
if (C1 == 0x0ull) {
|
|
// x is 0
|
|
res = 0x0000000000000001ull; // MINFP = 1 * 10^emin
|
|
} else { // x is not special and is not zero
|
|
if (x == 0x77fb86f26fc0ffffull) {
|
|
// x = +MAXFP = 999...99 * 10^emax
|
|
res = 0x7800000000000000ull; // +inf
|
|
} else if (x == 0x8000000000000001ull) {
|
|
// x = -MINFP = 1...99 * 10^emin
|
|
res = 0x8000000000000000ull; // -0
|
|
} else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
// can add/subtract 1 ulp to the significand
|
|
|
|
// Note: we could check here if x >= 10^16 to speed up the case q1 =16
|
|
// q1 = nr. of decimal digits in x (1 <= q1 <= 54)
|
|
// determine first the nr. of bits in x
|
|
if (C1 >= MASK_BINARY_OR2) { // x >= 2^53
|
|
// split the 64-bit value in two 32-bit halves to avoid rounding errors
|
|
if (C1 >= 0x0000000100000000ull) { // x >= 2^32
|
|
tmp1.d = (double) (C1 >> 32); // exact conversion
|
|
x_nr_bits =
|
|
33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
} else { // x < 2^32
|
|
tmp1.d = (double) C1; // exact conversion
|
|
x_nr_bits =
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
}
|
|
} else { // if x < 2^53
|
|
tmp1.d = (double) C1; // exact conversion
|
|
x_nr_bits =
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
}
|
|
q1 = nr_digits[x_nr_bits - 1].digits;
|
|
if (q1 == 0) {
|
|
q1 = nr_digits[x_nr_bits - 1].digits1;
|
|
if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
|
|
q1++;
|
|
}
|
|
// if q1 < P16 then pad the significand with zeros
|
|
if (q1 < P16) {
|
|
if (x_exp > (UINT64) (P16 - q1)) {
|
|
ind = P16 - q1; // 1 <= ind <= P16 - 1
|
|
// pad with P16 - q1 zeros, until exponent = emin
|
|
// C1 = C1 * 10^ind
|
|
C1 = C1 * ten2k64[ind];
|
|
x_exp = x_exp - ind;
|
|
} else { // pad with zeros until the exponent reaches emin
|
|
ind = x_exp;
|
|
C1 = C1 * ten2k64[ind];
|
|
x_exp = EXP_MIN;
|
|
}
|
|
}
|
|
if (!x_sign) { // x > 0
|
|
// add 1 ulp (add 1 to the significand)
|
|
C1++;
|
|
if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
|
|
C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
|
|
x_exp++;
|
|
}
|
|
// Ok, because MAXFP = 999...99 * 10^emax was caught already
|
|
} else { // x < 0
|
|
// subtract 1 ulp (subtract 1 from the significand)
|
|
C1--;
|
|
if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
|
|
C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
|
|
x_exp--;
|
|
}
|
|
}
|
|
// assemble the result
|
|
// if significand has 54 bits
|
|
if (C1 & MASK_BINARY_OR2) {
|
|
res =
|
|
x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
|
|
MASK_BINARY_SIG2);
|
|
} else { // significand fits in 53 bits
|
|
res = x_sign | (x_exp << 53) | C1;
|
|
}
|
|
} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
} // end x is not special and is not zero
|
|
BID_RETURN (res);
|
|
}
|
|
|
|
/*****************************************************************************
|
|
* BID64 nextdown
|
|
****************************************************************************/
|
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
void
|
|
bid64_nextdown (UINT64 * pres,
|
|
UINT64 *
|
|
px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
UINT64 x = *px;
|
|
#else
|
|
UINT64
|
|
bid64_nextdown (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
_EXC_INFO_PARAM) {
|
|
#endif
|
|
|
|
UINT64 res;
|
|
UINT64 x_sign;
|
|
UINT64 x_exp;
|
|
BID_UI64DOUBLE tmp1;
|
|
int x_nr_bits;
|
|
int q1, ind;
|
|
UINT64 C1; // C1 represents x_signif (UINT64)
|
|
|
|
// check for NaNs and infinities
|
|
if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
|
|
if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
|
|
x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
|
|
else
|
|
x = x & 0xfe03ffffffffffffull; // clear G6-G12
|
|
if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
|
|
// set invalid flag
|
|
*pfpsf |= INVALID_EXCEPTION;
|
|
// return quiet (SNaN)
|
|
res = x & 0xfdffffffffffffffull;
|
|
} else { // QNaN
|
|
res = x;
|
|
}
|
|
BID_RETURN (res);
|
|
} else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
|
|
if (x & 0x8000000000000000ull) { // x is -inf
|
|
res = 0xf800000000000000ull;
|
|
} else { // x is +inf
|
|
res = 0x77fb86f26fc0ffffull; // +MAXFP = +999...99 * 10^emax
|
|
}
|
|
BID_RETURN (res);
|
|
}
|
|
// unpack the argument
|
|
x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
|
// if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
|
|
if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
|
|
x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
|
|
C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
|
|
if (C1 > 9999999999999999ull) { // non-canonical
|
|
x_exp = 0;
|
|
C1 = 0;
|
|
}
|
|
} else {
|
|
x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
|
|
C1 = x & MASK_BINARY_SIG1;
|
|
}
|
|
|
|
// check for zeros (possibly from non-canonical values)
|
|
if (C1 == 0x0ull) {
|
|
// x is 0
|
|
res = 0x8000000000000001ull; // -MINFP = -1 * 10^emin
|
|
} else { // x is not special and is not zero
|
|
if (x == 0xf7fb86f26fc0ffffull) {
|
|
// x = -MAXFP = -999...99 * 10^emax
|
|
res = 0xf800000000000000ull; // -inf
|
|
} else if (x == 0x0000000000000001ull) {
|
|
// x = +MINFP = 1...99 * 10^emin
|
|
res = 0x0000000000000000ull; // -0
|
|
} else { // -MAXFP + 1ulp <= x <= -MINFP OR MINFP + 1 ulp <= x <= MAXFP
|
|
// can add/subtract 1 ulp to the significand
|
|
|
|
// Note: we could check here if x >= 10^16 to speed up the case q1 =16
|
|
// q1 = nr. of decimal digits in x (1 <= q1 <= 16)
|
|
// determine first the nr. of bits in x
|
|
if (C1 >= 0x0020000000000000ull) { // x >= 2^53
|
|
// split the 64-bit value in two 32-bit halves to avoid
|
|
// rounding errors
|
|
if (C1 >= 0x0000000100000000ull) { // x >= 2^32
|
|
tmp1.d = (double) (C1 >> 32); // exact conversion
|
|
x_nr_bits =
|
|
33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
} else { // x < 2^32
|
|
tmp1.d = (double) C1; // exact conversion
|
|
x_nr_bits =
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
}
|
|
} else { // if x < 2^53
|
|
tmp1.d = (double) C1; // exact conversion
|
|
x_nr_bits =
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
|
}
|
|
q1 = nr_digits[x_nr_bits - 1].digits;
|
|
if (q1 == 0) {
|
|
q1 = nr_digits[x_nr_bits - 1].digits1;
|
|
if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
|
|
q1++;
|
|
}
|
|
// if q1 < P16 then pad the significand with zeros
|
|
if (q1 < P16) {
|
|
if (x_exp > (UINT64) (P16 - q1)) {
|
|
ind = P16 - q1; // 1 <= ind <= P16 - 1
|
|
// pad with P16 - q1 zeros, until exponent = emin
|
|
// C1 = C1 * 10^ind
|
|
C1 = C1 * ten2k64[ind];
|
|
x_exp = x_exp - ind;
|
|
} else { // pad with zeros until the exponent reaches emin
|
|
ind = x_exp;
|
|
C1 = C1 * ten2k64[ind];
|
|
x_exp = EXP_MIN;
|
|
}
|
|
}
|
|
if (x_sign) { // x < 0
|
|
// add 1 ulp (add 1 to the significand)
|
|
C1++;
|
|
if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
|
|
C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
|
|
x_exp++;
|
|
// Ok, because -MAXFP = -999...99 * 10^emax was caught already
|
|
}
|
|
} else { // x > 0
|
|
// subtract 1 ulp (subtract 1 from the significand)
|
|
C1--;
|
|
if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
|
|
C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
|
|
x_exp--;
|
|
}
|
|
}
|
|
// assemble the result
|
|
// if significand has 54 bits
|
|
if (C1 & MASK_BINARY_OR2) {
|
|
res =
|
|
x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
|
|
MASK_BINARY_SIG2);
|
|
} else { // significand fits in 53 bits
|
|
res = x_sign | (x_exp << 53) | C1;
|
|
}
|
|
} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
|
} // end x is not special and is not zero
|
|
BID_RETURN (res);
|
|
}
|
|
|
|
/*****************************************************************************
|
|
* BID64 nextafter
|
|
****************************************************************************/
|
|
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
void
|
|
bid64_nextafter (UINT64 * pres, UINT64 * px,
|
|
UINT64 *
|
|
py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
|
UINT64 x = *px;
|
|
UINT64 y = *py;
|
|
#else
|
|
UINT64
|
|
bid64_nextafter (UINT64 x,
|
|
UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
|
_EXC_INFO_PARAM) {
|
|
#endif
|
|
|
|
UINT64 res;
|
|
UINT64 tmp1, tmp2;
|
|
FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
|
|
int res1, res2;
|
|
|
|
// check for NaNs or infinities
|
|
if (((x & MASK_SPECIAL) == MASK_SPECIAL) ||
|
|
((y & MASK_SPECIAL) == MASK_SPECIAL)) {
|
|
// x is NaN or infinity or y is NaN or infinity
|
|
|
|
if ((x & MASK_NAN) == MASK_NAN) { // x is NAN
|
|
if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
|
|
x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
|
|
else
|
|
x = x & 0xfe03ffffffffffffull; // clear G6-G12
|
|
if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
|
// set invalid flag
|
|
*pfpsf |= INVALID_EXCEPTION;
|
|
// return quiet (x)
|
|
res = x & 0xfdffffffffffffffull;
|
|
} else { // x is QNaN
|
|
if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
|
// set invalid flag
|
|
*pfpsf |= INVALID_EXCEPTION;
|
|
}
|
|
// return x
|
|
res = x;
|
|
}
|
|
BID_RETURN (res);
|
|
} else if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
|
|
if ((y & 0x0003ffffffffffffull) > 999999999999999ull)
|
|
y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
|
|
else
|
|
y = y & 0xfe03ffffffffffffull; // clear G6-G12
|
|
if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
|
// set invalid flag
|
|
*pfpsf |= INVALID_EXCEPTION;
|
|
// return quiet (y)
|
|
res = y & 0xfdffffffffffffffull;
|
|
} else { // y is QNaN
|
|
// return y
|
|
res = y;
|
|
}
|
|
BID_RETURN (res);
|
|
} else { // at least one is infinity
|
|
if ((x & MASK_ANY_INF) == MASK_INF) { // x = inf
|
|
x = x & (MASK_SIGN | MASK_INF);
|
|
}
|
|
if ((y & MASK_ANY_INF) == MASK_INF) { // y = inf
|
|
y = y & (MASK_SIGN | MASK_INF);
|
|
}
|
|
}
|
|
}
|
|
// neither x nor y is NaN
|
|
|
|
// if not infinity, check for non-canonical values x (treated as zero)
|
|
if ((x & MASK_ANY_INF) != MASK_INF) { // x != inf
|
|
// unpack x
|
|
if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
|
|
// if the steering bits are 11 (condition will be 0), then
|
|
// the exponent is G[0:w+1]
|
|
if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
|
|
9999999999999999ull) {
|
|
// non-canonical
|
|
x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
|
|
}
|
|
} else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) x is unch.
|
|
; // canonical
|
|
}
|
|
}
|
|
// no need to check for non-canonical y
|
|
|
|
// neither x nor y is NaN
|
|
tmp_fpsf = *pfpsf; // save fpsf
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
bid64_quiet_equal (&res1, px,
|
|
py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
bid64_quiet_greater (&res2, px,
|
|
py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#else
|
|
res1 =
|
|
bid64_quiet_equal (x,
|
|
y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
res2 =
|
|
bid64_quiet_greater (x,
|
|
y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#endif
|
|
*pfpsf = tmp_fpsf; // restore fpsf
|
|
if (res1) { // x = y
|
|
// return x with the sign of y
|
|
res = (y & 0x8000000000000000ull) | (x & 0x7fffffffffffffffull);
|
|
} else if (res2) { // x > y
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
bid64_nextdown (&res,
|
|
px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#else
|
|
res =
|
|
bid64_nextdown (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#endif
|
|
} else { // x < y
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
bid64_nextup (&res, px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#else
|
|
res = bid64_nextup (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
|
#endif
|
|
}
|
|
// if the operand x is finite but the result is infinite, signal
|
|
// overflow and inexact
|
|
if (((x & MASK_INF) != MASK_INF) && ((res & MASK_INF) == MASK_INF)) {
|
|
// set the inexact flag
|
|
*pfpsf |= INEXACT_EXCEPTION;
|
|
// set the overflow flag
|
|
*pfpsf |= OVERFLOW_EXCEPTION;
|
|
}
|
|
// if the result is in (-10^emin, 10^emin), and is different from the
|
|
// operand x, signal underflow and inexact
|
|
tmp1 = 0x00038d7ea4c68000ull; // +100...0[16] * 10^emin
|
|
tmp2 = res & 0x7fffffffffffffffull;
|
|
tmp_fpsf = *pfpsf; // save fpsf
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
|
bid64_quiet_greater (&res1, &tmp1,
|
|
&tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
_EXC_INFO_ARG);
|
|
bid64_quiet_not_equal (&res2, &x,
|
|
&res _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
_EXC_INFO_ARG);
|
|
#else
|
|
res1 =
|
|
bid64_quiet_greater (tmp1,
|
|
tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
_EXC_INFO_ARG);
|
|
res2 =
|
|
bid64_quiet_not_equal (x,
|
|
res _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
|
_EXC_INFO_ARG);
|
|
#endif
|
|
*pfpsf = tmp_fpsf; // restore fpsf
|
|
if (res1 && res2) {
|
|
// if (bid64_quiet_greater (tmp1, tmp2, &tmp_fpsf) &&
|
|
// bid64_quiet_not_equal (x, res, &tmp_fpsf)) {
|
|
// set the inexact flag
|
|
*pfpsf |= INEXACT_EXCEPTION;
|
|
// set the underflow flag
|
|
*pfpsf |= UNDERFLOW_EXCEPTION;
|
|
}
|
|
BID_RETURN (res);
|
|
}
|