123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643 |
- /* 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/>. */
- #define BID_128RES
- #include "bid_internal.h"
- /*****************************************************************************
- * BID128 nextup
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid128_nextup (UINT128 * pres,
- UINT128 *
- px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
- UINT128 x = *px;
- #else
- UINT128
- bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp;
- BID_UI64DOUBLE tmp1;
- int x_nr_bits;
- int q1, ind;
- UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
- BID_SWAP128 (x);
- // unpack the argument
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
- && (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- } else { // x is not NaN, so it must be infinity
- if (!x_sign) { // x is +inf
- res.w[1] = 0x7800000000000000ull; // +inf
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- res.w[1] = 0xdfffed09bead87c0ull; // -MAXFP = -999...99 * 10^emax
- res.w[0] = 0x378d8e63ffffffffull;
- }
- }
- BID_RETURN (res);
- }
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is +/-0
- res.w[1] = 0x0000000000000000ull; // +1 * 10^emin
- res.w[0] = 0x0000000000000001ull;
- } else { // x is not special and is not zero
- if (x.w[1] == 0x5fffed09bead87c0ull
- && x.w[0] == 0x378d8e63ffffffffull) {
- // x = +MAXFP = 999...99 * 10^emax
- res.w[1] = 0x7800000000000000ull; // +inf
- res.w[0] = 0x0000000000000000ull;
- } else if (x.w[1] == 0x8000000000000000ull
- && x.w[0] == 0x0000000000000001ull) {
- // x = -MINFP = 1...99 * 10^emin
- res.w[1] = 0x8000000000000000ull; // -0
- res.w[0] = 0x0000000000000000ull;
- } 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^34 to speed up the case q1 = 34
- // q1 = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rnd errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
- 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
- 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
- || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
- && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q1++;
- }
- // if q1 < P34 then pad the significand with zeros
- if (q1 < P34) {
- exp = (x_exp >> 49) - 6176;
- if (exp + 6176 > P34 - q1) {
- ind = P34 - q1; // 1 <= ind <= P34 - 1
- // pad with P34 - q1 zeros, until exponent = emin
- // C1 = C1 * 10^ind
- if (q1 <= 19) { // 64-bit C1
- if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 128-bit 10^ind and 64-bit C1
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- } else { // C1 is (most likely) 128-bit
- if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
- __mul_128x64_to_128 (C1, ten2k64[ind], C1);
- } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- }
- x_exp = x_exp - ((UINT64) ind << 49);
- } else { // pad with zeros until the exponent reaches emin
- ind = exp + 6176;
- // C1 = C1 * 10^ind
- if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
- if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
- __mul_128x64_to_128 (C1, ten2k64[ind], C1);
- }
- } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
- // 64-bit C1, 128-bit 10^ind
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- x_exp = EXP_MIN;
- }
- }
- if (!x_sign) { // x > 0
- // add 1 ulp (add 1 to the significand)
- C1.w[0]++;
- if (C1.w[0] == 0)
- C1.w[1]++;
- if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
- C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
- C1.w[0] = 0x38c15b0a00000000ull;
- x_exp = x_exp + EXP_P1;
- }
- } else { // x < 0
- // subtract 1 ulp (subtract 1 from the significand)
- C1.w[0]--;
- if (C1.w[0] == 0xffffffffffffffffull)
- C1.w[1]--;
- if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
- C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
- C1.w[0] = 0x378d8e63ffffffffull;
- x_exp = x_exp - EXP_P1;
- }
- }
- // assemble the result
- res.w[1] = x_sign | x_exp | C1.w[1];
- res.w[0] = C1.w[0];
- } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
- } // end x is not special and is not zero
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID128 nextdown
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid128_nextdown (UINT128 * pres,
- UINT128 *
- px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
- UINT128 x = *px;
- #else
- UINT128
- bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp;
- BID_UI64DOUBLE tmp1;
- int x_nr_bits;
- int q1, ind;
- UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
- BID_SWAP128 (x);
- // unpack the argument
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
- && (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- } else { // x is not NaN, so it must be infinity
- if (!x_sign) { // x is +inf
- res.w[1] = 0x5fffed09bead87c0ull; // +MAXFP = +999...99 * 10^emax
- res.w[0] = 0x378d8e63ffffffffull;
- } else { // x is -inf
- res.w[1] = 0xf800000000000000ull; // -inf
- res.w[0] = 0x0000000000000000ull;
- }
- }
- BID_RETURN (res);
- }
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is +/-0
- res.w[1] = 0x8000000000000000ull; // -1 * 10^emin
- res.w[0] = 0x0000000000000001ull;
- } else { // x is not special and is not zero
- if (x.w[1] == 0xdfffed09bead87c0ull
- && x.w[0] == 0x378d8e63ffffffffull) {
- // x = -MAXFP = -999...99 * 10^emax
- res.w[1] = 0xf800000000000000ull; // -inf
- res.w[0] = 0x0000000000000000ull;
- } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) { // +MINFP
- res.w[1] = 0x0000000000000000ull; // +0
- res.w[0] = 0x0000000000000000ull;
- } 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^34 to speed up the case q1 = 34
- // q1 = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rnd errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
- 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
- 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
- || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
- && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q1++;
- }
- // if q1 < P then pad the significand with zeros
- if (q1 < P34) {
- exp = (x_exp >> 49) - 6176;
- if (exp + 6176 > P34 - q1) {
- ind = P34 - q1; // 1 <= ind <= P34 - 1
- // pad with P34 - q1 zeros, until exponent = emin
- // C1 = C1 * 10^ind
- if (q1 <= 19) { // 64-bit C1
- if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 128-bit 10^ind and 64-bit C1
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- } else { // C1 is (most likely) 128-bit
- if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
- __mul_128x64_to_128 (C1, ten2k64[ind], C1);
- } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- }
- x_exp = x_exp - ((UINT64) ind << 49);
- } else { // pad with zeros until the exponent reaches emin
- ind = exp + 6176;
- // C1 = C1 * 10^ind
- if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
- if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
- __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
- } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
- __mul_128x64_to_128 (C1, ten2k64[ind], C1);
- }
- } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
- // 64-bit C1, 128-bit 10^ind
- __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
- }
- x_exp = EXP_MIN;
- }
- }
- if (x_sign) { // x < 0
- // add 1 ulp (add 1 to the significand)
- C1.w[0]++;
- if (C1.w[0] == 0)
- C1.w[1]++;
- if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
- C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
- C1.w[0] = 0x38c15b0a00000000ull;
- x_exp = x_exp + EXP_P1;
- }
- } else { // x > 0
- // subtract 1 ulp (subtract 1 from the significand)
- C1.w[0]--;
- if (C1.w[0] == 0xffffffffffffffffull)
- C1.w[1]--;
- if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
- C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
- C1.w[0] = 0x378d8e63ffffffffull;
- x_exp = x_exp - EXP_P1;
- }
- }
- // assemble the result
- res.w[1] = x_sign | x_exp | C1.w[1];
- res.w[0] = C1.w[0];
- } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
- } // end x is not special and is not zero
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID128 nextafter
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid128_nextafter (UINT128 * pres, UINT128 * px,
- UINT128 *
- py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
- {
- UINT128 x = *px;
- UINT128 y = *py;
- UINT128 xnswp = *px;
- UINT128 ynswp = *py;
- #else
- UINT128
- bid128_nextafter (UINT128 x,
- UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT128 xnswp = x;
- UINT128 ynswp = y;
- #endif
- UINT128 res;
- UINT128 tmp1, tmp2, tmp3;
- FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
- int res1, res2;
- UINT64 x_exp;
- BID_SWAP128 (x);
- BID_SWAP128 (y);
- // check for NaNs
- if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
- || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
- // x is special or y is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
- && (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- }
- }
- BID_RETURN (res)
- } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
- && (y.w[0] > 0x38c15b09ffffffffull))) {
- y.w[1] = y.w[1] & 0xffffc00000000000ull;
- y.w[0] = 0x0ull;
- }
- if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = y.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = y.w[0];
- }
- BID_RETURN (res)
- } else { // at least one is infinity
- if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
- x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
- x.w[0] = 0x0ull;
- }
- if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
- y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
- y.w[0] = 0x0ull;
- }
- }
- }
- // neither x nor y is NaN
- // if not infinity, check for non-canonical values x (treated as zero)
- if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
- x.w[0] = 0x0ull;
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
- ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
- && x.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
- x.w[0] = 0x0ull;
- } else { // 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
- bid128_quiet_equal (&res1, &xnswp,
- &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- bid128_quiet_greater (&res2, &xnswp,
- &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- #else
- res1 =
- bid128_quiet_equal (xnswp,
- ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- res2 =
- bid128_quiet_greater (xnswp,
- ynswp _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.w[1] =
- (x.w[1] & 0x7fffffffffffffffull) | (y.
- w[1] & 0x8000000000000000ull);
- res.w[0] = x.w[0];
- } else if (res2) { // x > y
- #if DECIMAL_CALL_BY_REFERENCE
- bid128_nextdown (&res,
- &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- #else
- res =
- bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- #endif
- BID_SWAP128 (res);
- } else { // x < y
- #if DECIMAL_CALL_BY_REFERENCE
- bid128_nextup (&res,
- &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
- #else
- res =
- bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
- #endif
- BID_SWAP128 (res);
- }
- // if the operand x is finite but the result is infinite, signal
- // overflow and inexact
- if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
- && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
- // 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.w[HIGH_128W] = 0x0000314dc6448d93ull;
- tmp1.w[LOW_128W] = 0x38c15b0a00000000ull; // +100...0[34] * 10^emin
- tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
- tmp2.w[LOW_128W] = res.w[0];
- tmp3.w[HIGH_128W] = res.w[1];
- tmp3.w[LOW_128W] = res.w[0];
- tmp_fpsf = *pfpsf; // save fpsf
- #if DECIMAL_CALL_BY_REFERENCE
- bid128_quiet_greater (&res1, &tmp1,
- &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- bid128_quiet_not_equal (&res2, &xnswp,
- &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- #else
- res1 =
- bid128_quiet_greater (tmp1,
- tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- res2 =
- bid128_quiet_not_equal (xnswp,
- tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
- _EXC_INFO_ARG);
- #endif
- *pfpsf = tmp_fpsf; // restore fpsf
- if (res1 && res2) {
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // set the underflow flag
- *pfpsf |= UNDERFLOW_EXCEPTION;
- }
- BID_RETURN (res);
- }
|