]> oss.titaniummirror.com Git - msp430-gcc.git/blobdiff - libgcc/config/libbid/bid128_next.c
Imported gcc-4.4.3
[msp430-gcc.git] / libgcc / config / libbid / bid128_next.c
diff --git a/libgcc/config/libbid/bid128_next.c b/libgcc/config/libbid/bid128_next.c
new file mode 100644 (file)
index 0000000..11688c1
--- /dev/null
@@ -0,0 +1,643 @@
+/* Copyright (C) 2007, 2009  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);
+}