| /* Copyright (C) 2007-2021 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); |
| } |