Viewing file:  bid64_next.c (15.33 KB) -rw-r--r--
Select action/file-type:
 (+) |  (+) |  (+) | Code (+) | Session (+) |  (+) | SDB (+) |  (+) |  (+) |  (+) |  (+) |  (+) |
/* Copyright (C) 2007-2022 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);
}
|
