--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libgcc/config/libbid/bid128_next.c	Fri Jul 17 14:47:48 2009 +0900
@@ -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);
+}