Fix _Decimal128 arithmetic error under FE_UPWARD.

libgcc/config/libbid/ChangeLog:

	PR target/120691
	* bid128_div.c: Fix _Decimal128 arithmetic error under
	FE_UPWARD.
	* bid128_rem.c: Ditto.
	* bid128_sqrt.c: Ditto.
	* bid64_div.c (bid64_div): Ditto.
	* bid64_sqrt.c (bid64_sqrt): Ditto.

gcc/testsuite/ChangeLog:

	* gcc.target/i386/pr120691.c: New test.

6 files changed, 586 insertions, 215 deletions

diff --git a/gcc/testsuite/gcc.target/i386/pr120691.c b/gcc/testsuite/gcc.target/i386/pr120691.c
new file mode 100644
index 0000000..241a34f
--- /dev/null
+++ b/gcc/testsuite/gcc.target/i386/pr120691.c
@@ -0,0 +1,18 @@
+/* { dg-do run { target { ! ia32 } } } */
+/* { dg-options "-O0 -mfpmath=sse" } */
+/* { dg-require-effective-target fenv } */
+/* { dg-require-effective-target dfp } */
+
+#include <fenv.h>
+
+int   main()   {
+  fesetround( FE_UPWARD );
+  _Decimal128   x1 =  9825,  x2 =  10000 ;
+
+  double c = (double) (x1 / x2);
+
+  if (c != 0.9825)
+    __builtin_abort ();
+
+  return   0 ;
+}
diff --git a/libgcc/config/libbid/bid128_div.c b/libgcc/config/libbid/bid128_div.c
index 925bf14..04955f2 100644..100755
--- a/libgcc/config/libbid/bid128_div.c
+++ b/libgcc/config/libbid/bid128_div.c
@@ -24,15 +24,16 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 #define BID_128RES
 #include "bid_div_macros.h"
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
-#include <fenv.h>
-
 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
 #endif
+#include <stdio.h>
+#include <fenv.h>
 
 extern UINT32 convert_table[5][128][2];
 extern SINT8 factors[][2];
 extern UINT8 packed_10000_zeros[];
 
+
 BID128_FUNCTION_ARG2 (bid128_div, x, y)
 
      UINT256 CA4, CA4r, P256;
@@ -45,10 +46,18 @@ BID128_FUNCTION_ARG2 (bid128_div, x, y)
        digits_q, amount;
      int nzeros, i, j, k, d5;
      unsigned rmode;
+     int old_rm, rm_changed=0;
+
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+  
 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
 
   // unpack arguments, check for NaN or Infinity
@@ -62,6 +71,8 @@ if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 #endif
   res.w[1] = (CX.w[1]) & QUIET_MASK64;
   res.w[0] = CX.w[0];
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
     // x is Infinity?
@@ -75,6 +86,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
 #endif
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is NaN?
@@ -85,6 +98,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
     res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
       0x7800000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -97,6 +112,8 @@ if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
     // x=y=0, return NaN
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // return 0
@@ -108,6 +125,8 @@ if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
     exponent_x = 0;
   res.w[1] |= (((UINT64) exponent_x) << 49);
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 }
@@ -122,6 +141,8 @@ if (!valid_y) {
 #endif
     res.w[1] = CY.w[1] & QUIET_MASK64;
     res.w[0] = CY.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is Infinity?
@@ -129,6 +150,8 @@ if (!valid_y) {
     // return +/-0
     res.w[1] = sign_x ^ sign_y;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -138,6 +161,8 @@ if (!valid_y) {
   res.w[1] =
     ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -186,6 +211,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // get number of decimal digits in CQ
@@ -379,6 +406,8 @@ if (!CA4.w[0] && !CA4.w[1])
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #endif
@@ -469,6 +498,8 @@ if (diff_expon >= 0) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 
 }
@@ -477,6 +508,8 @@ get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 
@@ -496,10 +529,17 @@ TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128, bid128dd_div, UINT64, x,
        digits_q, amount;
      int nzeros, i, j, k, d5;
      unsigned rmode;
+     int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
 
 	// unpack arguments, check for NaN or Infinity
@@ -519,6 +559,8 @@ if ((x & NAN_MASK64) == NAN_MASK64) {
   res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
   __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
   res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 	   // x is Infinity?
@@ -532,13 +574,17 @@ if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
 #endif
   res.w[1] = 0x7c00000000000000ull;
   res.w[0] = 0;
-    BID_RETURN (res);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
+  BID_RETURN (res);
   }
   if ((((y) & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
   // otherwise return +/-Inf
   res.w[1] =
     (((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
   }
 }
@@ -551,6 +597,8 @@ if ((((y) & 0x7800000000000000ull) != 0x7800000000000000ull)) {
   // x=y=0, return NaN
   res.w[1] = 0x7c00000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 	   // return 0
@@ -566,6 +614,8 @@ else if (exponent_x < 0)
   exponent_x = 0;
 res.w[1] |= (((UINT64) exponent_x) << 49);
 res.w[0] = 0;
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 }
@@ -583,13 +633,17 @@ if (!valid_y) {
   res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
   __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
   res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
-    BID_RETURN (res);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
+  BID_RETURN (res);
   }
   // y is Infinity?
   if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
     // return +/-0
     res.w[1] = sign_x ^ sign_y;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -599,6 +653,8 @@ if (!valid_y) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -647,6 +703,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // get number of decimal digits in CQ
@@ -843,6 +901,8 @@ __div_256_by_128 (&CQ, &CA4, CY);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 #endif
@@ -932,6 +992,8 @@ if (diff_expon >= 0) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 
 }
@@ -940,6 +1002,8 @@ get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 
@@ -955,10 +1019,17 @@ BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div, UINT64, x, y)
        digits_q, amount;
      int nzeros, i, j, k, d5;
      unsigned rmode;
+     int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
 
 	// unpack arguments, check for NaN or Infinity
@@ -978,6 +1049,8 @@ if ((x & NAN_MASK64) == NAN_MASK64) {
   res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
   __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
   res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 	   // x is Infinity?
@@ -991,6 +1064,8 @@ if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
 #endif
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
@@ -998,6 +1073,8 @@ if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
   res.w[1] =
     ((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
   }
 }
@@ -1010,6 +1087,8 @@ if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
     // x=y=0, return NaN
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // return 0
@@ -1021,6 +1100,8 @@ if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
     exponent_x = 0;
   res.w[1] |= (((UINT64) exponent_x) << 49);
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 }
@@ -1037,6 +1118,8 @@ if (!valid_y) {
 #endif
     res.w[1] = CY.w[1] & QUIET_MASK64;
     res.w[0] = CY.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is Infinity?
@@ -1044,6 +1127,8 @@ if (!valid_y) {
     // return +/-0
     res.w[1] = sign_x ^ sign_y;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -1053,6 +1138,8 @@ if (!valid_y) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -1101,6 +1188,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // get number of decimal digits in CQ
@@ -1300,6 +1389,8 @@ __div_256_by_128 (&CQ, &CA4, CY);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 #endif
@@ -1389,6 +1480,8 @@ if (diff_expon >= 0) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 
@@ -1396,6 +1489,8 @@ get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 
 }
@@ -1411,10 +1506,16 @@ BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div, x, UINT64, y)
      int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
        digits_q, amount;
      int nzeros, i, j, k, d5, rmode;
+     int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
 
 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
 	// unpack arguments, check for NaN or Infinity
@@ -1428,6 +1529,8 @@ if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 #endif
   res.w[1] = (CX.w[1]) & QUIET_MASK64;
   res.w[0] = CX.w[0];
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
     // x is Infinity?
@@ -1441,6 +1544,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
 #endif
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is NaN?
@@ -1451,6 +1556,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
     res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull) |
       0x7800000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -1463,6 +1570,8 @@ if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
     // x=y=0, return NaN
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // return 0
@@ -1474,6 +1583,8 @@ if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
     exponent_x = 0;
   res.w[1] |= (((UINT64) exponent_x) << 49);
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 }
@@ -1490,13 +1601,17 @@ if (!valid_y) {
   res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
   __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
   res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
-    BID_RETURN (res);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
+  BID_RETURN (res);
   }
   // y is Infinity?
   if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
     // return +/-0
     res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull);
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0
@@ -1505,6 +1620,8 @@ if (!valid_y) {
 #endif
   res.w[1] = (sign_x ^ sign_y) | INFINITY_MASK64;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -1553,6 +1670,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // get number of decimal digits in CQ
@@ -1749,6 +1868,8 @@ __div_256_by_128 (&CQ, &CA4, CY);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 #endif
@@ -1838,6 +1959,8 @@ if (diff_expon >= 0) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 
 }
@@ -1846,6 +1969,8 @@ get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 
 }
diff --git a/libgcc/config/libbid/bid128_rem.c b/libgcc/config/libbid/bid128_rem.c
index f229b28..68a88d2 100644..100755
--- a/libgcc/config/libbid/bid128_rem.c
+++ b/libgcc/config/libbid/bid128_rem.c
@@ -23,6 +23,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 
 #define BID_128RES
 #include "bid_div_macros.h"
+#include <fenv.h>
 
 
 BID128_FUNCTION_ARG2_NORND_CUSTOMRESTYPE (UINT128, bid128_rem, x, y)
@@ -34,9 +35,16 @@ BID128_FUNCTION_ARG2_NORND_CUSTOMRESTYPE (UINT128, bid128_rem, x, y)
      int_float f64, fx;
      int exponent_x, exponent_y, diff_expon, bin_expon_cx, scale,
        scale0;
+     int old_rm, rm_changed=0;
 
   // unpack arguments, check for NaN or Infinity
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
 
 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
@@ -52,6 +60,8 @@ if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 #endif
   res.w[1] = CX.w[1] & QUIET_MASK64;
   res.w[0] = CX.w[0];
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
     // x is Infinity?
@@ -66,6 +76,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
 #endif
     res.w[1] = 0x7c00000000000000ull;
     res.w[0] = 0;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 
@@ -79,6 +91,8 @@ if ((!CY.w[1]) && (!CY.w[0])) {
   // x=y=0, return NaN
   res.w[1] = 0x7c00000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 if (valid_y || ((y.w[1] & NAN_MASK64) == INFINITY_MASK64)) {
@@ -89,6 +103,8 @@ if (valid_y || ((y.w[1] & NAN_MASK64) == INFINITY_MASK64)) {
 
   res.w[1] = sign_x | (((UINT64) exponent_x) << 49);
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 }
@@ -103,6 +119,8 @@ if (!valid_y) {
 #endif
     res.w[1] = CY.w[1] & QUIET_MASK64;
     res.w[0] = CY.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is Infinity?
@@ -110,6 +128,8 @@ if (!valid_y) {
     // return x
     res.w[1] = x.w[1];
     res.w[0] = x.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0
@@ -119,6 +139,8 @@ if (!valid_y) {
 #endif
   res.w[1] = 0x7c00000000000000ull;
   res.w[0] = 0;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 
@@ -130,6 +152,8 @@ if (diff_expon <= 0) {
   if (diff_expon > 34) {
     // |x|<|y| in this case
     res = x;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // set exponent of y to exponent_x, scale coefficient_y
@@ -139,6 +163,8 @@ if (diff_expon <= 0) {
   if (P256.w[2] || P256.w[3]) {
     // |x|<|y| in this case
     res = x;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 
@@ -147,6 +173,8 @@ if (diff_expon <= 0) {
   if (__unsigned_compare_ge_128 (P256, CX2)) {
     // |x|<|y| in this case
     res = x;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 
@@ -164,6 +192,8 @@ if (diff_expon <= 0) {
   }
 
   get_BID128_very_fast (&res, sign_x, exponent_x, CR);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
   // 2^64
@@ -200,6 +230,8 @@ while (diff_expon > 0) {
   // check for remainder == 0
   if (!CX.w[1] && !CX.w[0]) {
     get_BID128_very_fast (&res, sign_x, exponent_y, CX);
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -213,5 +245,7 @@ if ((__unsigned_compare_gt_128 (CX2, CY))
 }
 
 get_BID128_very_fast (&res, sign_x, exponent_y, CX);
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
diff --git a/libgcc/config/libbid/bid128_sqrt.c b/libgcc/config/libbid/bid128_sqrt.c
index b280383..577ece9 100644..100755
--- a/libgcc/config/libbid/bid128_sqrt.c
+++ b/libgcc/config/libbid/bid128_sqrt.c
@@ -24,9 +24,8 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 #define BID_128RES
 #include "bid_internal.h"
 #include "bid_sqrt_macros.h"
-#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #include <fenv.h>
-
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
 #endif
 
@@ -39,10 +38,17 @@ BID128_FUNCTION_ARG1 (bid128_sqrt, x)
      int_float fx, f64;
      int exponent_x, bin_expon_cx;
      int digits, scale, exponent_q;
+     int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
   // unpack arguments, check for NaN or Infinity
 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
 res.w[1] = CX.w[1];
@@ -54,6 +60,8 @@ if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
     __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
   res.w[1] = CX.w[1] & QUIET_MASK64;
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
     // x is Infinity?
@@ -66,6 +74,8 @@ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
     __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
   }
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
     // x is 0 otherwise
@@ -74,6 +84,8 @@ res.w[1] =
   sign_x |
   ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) << 49);
 res.w[0] = 0;
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 if (sign_x) {
@@ -82,6 +94,8 @@ if (sign_x) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -117,6 +131,8 @@ if (CS.w[0] * CS.w[0] == A10.w[0]) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -288,6 +304,8 @@ get_BID128_fast (&res, 0,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 
@@ -302,10 +320,18 @@ BID128_FUNCTION_ARGTYPE1 (bid128d_sqrt, UINT64, x)
      int_float fx, f64;
      int exponent_x, bin_expon_cx;
      int digits, scale, exponent_q;
+     int old_rm, rm_changed=0;
+
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
      fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 	// unpack arguments, check for NaN or Infinity
    // unpack arguments, check for NaN or Infinity
 CX.w[1] = 0;
@@ -321,6 +347,8 @@ if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
   res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
   __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
   res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 	   // x is Infinity?
@@ -332,6 +360,8 @@ if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
     __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
   }
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 	   // x is 0 otherwise
@@ -342,6 +372,8 @@ res.w[1] =
   sign_x | ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1)
 	    << 49);
 res.w[0] = 0;
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 }
 if (sign_x) {
@@ -350,6 +382,8 @@ if (sign_x) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -387,6 +421,8 @@ if (CS.w[0] * CS.w[0] == A10.w[0]) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -558,7 +594,8 @@ get_BID128_fast (&res, 0, (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 
-
 }
diff --git a/libgcc/config/libbid/bid64_div.c b/libgcc/config/libbid/bid64_div.c
index 6975848..7f34556 100644..100755
--- a/libgcc/config/libbid/bid64_div.c
+++ b/libgcc/config/libbid/bid64_div.c
@@ -55,9 +55,8 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 
 #include "bid_internal.h"
 #include "bid_div_macros.h"
-#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #include <fenv.h>
-
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
 #endif
 
@@ -65,7 +64,6 @@ extern UINT32 convert_table[5][128][2];
 extern SINT8 factors[][2];
 extern UINT8 packed_10000_zeros[];
 
-
 #if DECIMAL_CALL_BY_REFERENCE
 
 void
@@ -94,6 +92,7 @@ bid64_div (UINT64 x,
   int rmode, amount;
   int nzeros, i, j, k, d5;
   UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
+  int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   fexcept_t binaryflags = 0;
 #endif
@@ -106,6 +105,12 @@ bid64_div (UINT64 x,
   y = *py;
 #endif
 
+  // Set it to round-to-nearest (if different)
+  if ((old_rm=fegetround()) != FE_TONEAREST) {
+    rm_changed=1;
+    fesetround(FE_TONEAREST);
+  }
+
   valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
   valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
 
@@ -123,6 +128,8 @@ bid64_div (UINT64 x,
       if ((x & SNAN_MASK64) == SNAN_MASK64)	// sNaN
 	__set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (coefficient_x & QUIET_MASK64);
     }
     // x is Infinity?
@@ -134,10 +141,14 @@ bid64_div (UINT64 x,
 #ifdef SET_STATUS_FLAGS
 	  __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+          // restore the rounding mode back if it has been changed
+          if (rm_changed) fesetround(old_rm);
 	  BID_RETURN (NAN_MASK64);
 	}
       } else {
 	// otherwise return +/-Inf
+        // restore the rounding mode back if it has been changed
+        if (rm_changed) fesetround(old_rm);
 	BID_RETURN (((x ^ y) & 0x8000000000000000ull) |
 		    INFINITY_MASK64);
       }
@@ -149,6 +160,8 @@ bid64_div (UINT64 x,
 #ifdef SET_STATUS_FLAGS
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (NAN_MASK64);
     }
     if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
@@ -163,6 +176,8 @@ bid64_div (UINT64 x,
 	exponent_x = DECIMAL_MAX_EXPON_64;
       else if (exponent_x < 0)
 	exponent_x = 0;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
     }
 
@@ -176,17 +191,23 @@ bid64_div (UINT64 x,
       if ((y & SNAN_MASK64) == SNAN_MASK64)	// sNaN
 	__set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (coefficient_y & QUIET_MASK64);
     }
     // y is Infinity?
     if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
       // return +/-0
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (((x ^ y) & 0x8000000000000000ull));
     }
     // y is 0
 #ifdef SET_STATUS_FLAGS
     __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN ((sign_x ^ sign_y) | INFINITY_MASK64);
   }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -255,6 +276,8 @@ bid64_div (UINT64 x,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
     // get decimal digits of Q
@@ -424,6 +447,8 @@ bid64_div (UINT64 x,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
   }
@@ -494,6 +519,8 @@ bid64_div (UINT64 x,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   } else {
     // UF occurs
@@ -510,6 +537,8 @@ bid64_div (UINT64 x,
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
 
   }
@@ -528,10 +557,17 @@ int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
   digits_q, amount;
 int nzeros, i, j, k, d5, done = 0;
 unsigned rmode;
+int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
 
 	// unpack arguments, check for NaN or Infinity
@@ -545,6 +581,8 @@ if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
   // test if x is NaN
   if (((x) & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
     res = CX.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res & QUIET_MASK64);
   }
   // x is Infinity?
@@ -557,14 +595,18 @@ if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
       res = 0x7c00000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
+    }
+    if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
+      // otherwise return +/-Inf
+      res =
+	(((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
-	if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
-    // otherwise return +/-Inf
-    res =
-      (((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
-    BID_RETURN (res);
-	}
   }
   // x is 0
   if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
@@ -574,6 +616,8 @@ if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
 #endif
       // x=y=0, return NaN
       res = 0x7c00000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
     // return 0
@@ -584,6 +628,8 @@ if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
     else if (exponent_x < 0)
       exponent_x = 0;
     res |= (((UINT64) exponent_x) << 53);
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -604,12 +650,16 @@ if (!valid_y) {
     amount = recip_scale[18];
     __shr_128 (Tmp, Qh, amount);
     res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is Infinity?
   if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
     // return +/-0
     res = sign_x ^ sign_y;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -618,6 +668,8 @@ if (!valid_y) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -680,6 +732,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
 
@@ -816,15 +870,17 @@ if (!done) {
 
       }
     }
-	if(diff_expon>=0){
-    res =
-      fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
-			       rnd_mode, pfpsf);
+    if(diff_expon>=0){
+      res =
+	fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
+				 rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
-    BID_RETURN (res);
-	}
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
+    }
   }
 #endif
 
@@ -905,6 +961,8 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   } else {
     // UF occurs
@@ -921,8 +979,9 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
-
   }
 
 }
@@ -942,10 +1001,18 @@ int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
   digits_q, amount;
 int nzeros, i, j, k, d5, done = 0;
 unsigned rmode;
+int old_rm, rm_changed=0;
+
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], (y));
 
 	// unpack arguments, check for NaN or Infinity
@@ -964,7 +1031,9 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
       amount = recip_scale[18];
       __shr_128 (Tmp, Qh, amount);
       res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
-    BID_RETURN (res);
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
   }
   // x is Infinity?
   if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
@@ -976,14 +1045,18 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
       res = 0x7c00000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
+    }
+    if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
+      // otherwise return +/-Inf
+      res =
+	((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
-	if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
-    // otherwise return +/-Inf
-    res =
-      ((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
-    BID_RETURN (res);
-	}
   }
   // x is 0
   if (((y & INFINITY_MASK64) != INFINITY_MASK64) &&
@@ -993,17 +1066,21 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
 #endif
     // x=y=0, return NaN
     res = 0x7c00000000000000ull;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // return 0
   if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)) {
-	  if (!CY.w[0]) {
+    if (!CY.w[0]) {
 #ifdef SET_STATUS_FLAGS
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
       res = 0x7c00000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
-	  }
+    }
     exponent_x =
       exponent_x - exponent_y - DECIMAL_EXPONENT_BIAS_128 +
       (DECIMAL_EXPONENT_BIAS << 1);
@@ -1012,6 +1089,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
     else if (exponent_x < 0)
       exponent_x = 0;
     res = (sign_x ^ sign_y) | (((UINT64) exponent_x) << 53);
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 }
@@ -1025,12 +1104,16 @@ if (!valid_y) {
     if ((y & SNAN_MASK64) == SNAN_MASK64)	// sNaN
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (CY.w[0] & QUIET_MASK64);
   }
   // y is Infinity?
   if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
     // return +/-0
     res = sign_x ^ sign_y;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -1039,6 +1122,8 @@ if (!valid_y) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -1103,6 +1188,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
 
@@ -1159,102 +1246,104 @@ if (!done) {
 #endif
 #ifndef LEAVE_TRAILING_ZEROS
     // check whether result is exact
-  {
-	  if(!done) {
-    // check whether CX, CY are short
-    if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
-      i = (int) CY.w[0] - 1;
-      j = (int) CX.w[0] - 1;
-      // difference in powers of 2 factors for Y and X
-      nzeros = ed2 - factors[i][0] + factors[j][0];
-      // difference in powers of 5 factors
-      d5 = ed2 - factors[i][1] + factors[j][1];
-      if (d5 < nzeros)
-		nzeros = d5;
-      // get P*(2^M[extra_digits])/10^extra_digits
-      __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
-      //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
-
-      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
-      amount = recip_scale[nzeros];
-      __shr_128_long (CQ, Qh, amount);
-
-      diff_expon += nzeros;
-    } else {
-      // decompose Q as Qh*10^17 + Ql
-      //T128 = reciprocals10_128[17];
-      Q_low = CQ.w[0];
-
-      {
-	tdigit[0] = Q_low & 0x3ffffff;
-	tdigit[1] = 0;
-	QX = Q_low >> 26;
-	QX32 = QX;
-	nzeros = 0;
-
-	for (j = 0; QX32; j++, QX32 >>= 7) {
-	  k = (QX32 & 127);
-	  tdigit[0] += convert_table[j][k][0];
-	  tdigit[1] += convert_table[j][k][1];
-	  if (tdigit[0] >= 100000000) {
-	    tdigit[0] -= 100000000;
-	    tdigit[1]++;
-	  }
-	}
-
-	if (tdigit[1] >= 100000000) {
-	  tdigit[1] -= 100000000;
-	  if (tdigit[1] >= 100000000)
-	    tdigit[1] -= 100000000;
-	}
-
-	digit = tdigit[0];
-	if (!digit && !tdigit[1])
-	  nzeros += 16;
-	else {
-	  if (!digit) {
-	    nzeros += 8;
-	    digit = tdigit[1];
-	  }
-	  // decompose digit
-	  PD = (UINT64) digit *0x068DB8BBull;
-	  digit_h = (UINT32) (PD >> 40);
-	  digit_low = digit - digit_h * 10000;
-
-	  if (!digit_low)
-	    nzeros += 4;
-	  else
-	    digit_h = digit_low;
-
-	  if (!(digit_h & 1))
-	    nzeros +=
-	      3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
-			    (digit_h & 7));
-	}
-
-	if (nzeros) {
+    {
+      if(!done) {
+	// check whether CX, CY are short
+	if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
+	  i = (int) CY.w[0] - 1;
+	  j = (int) CX.w[0] - 1;
+	  // difference in powers of 2 factors for Y and X
+	  nzeros = ed2 - factors[i][0] + factors[j][0];
+	  // difference in powers of 5 factors
+	  d5 = ed2 - factors[i][1] + factors[j][1];
+	  if (d5 < nzeros)
+	    nzeros = d5;
 	  // get P*(2^M[extra_digits])/10^extra_digits
 	  __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
+	  //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
 
 	  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 	  amount = recip_scale[nzeros];
-	  __shr_128 (CQ, Qh, amount);
-	}
-	diff_expon += nzeros;
+	  __shr_128_long (CQ, Qh, amount);
+
+	  diff_expon += nzeros;
+	} else {
+	  // decompose Q as Qh*10^17 + Ql
+	  //T128 = reciprocals10_128[17];
+	  Q_low = CQ.w[0];
+
+	  {
+	    tdigit[0] = Q_low & 0x3ffffff;
+	    tdigit[1] = 0;
+	    QX = Q_low >> 26;
+	    QX32 = QX;
+	    nzeros = 0;
+
+	    for (j = 0; QX32; j++, QX32 >>= 7) {
+	      k = (QX32 & 127);
+	      tdigit[0] += convert_table[j][k][0];
+	      tdigit[1] += convert_table[j][k][1];
+	      if (tdigit[0] >= 100000000) {
+		tdigit[0] -= 100000000;
+		tdigit[1]++;
+	      }
+	    }
+
+	    if (tdigit[1] >= 100000000) {
+	      tdigit[1] -= 100000000;
+	      if (tdigit[1] >= 100000000)
+		tdigit[1] -= 100000000;
+	    }
+
+	    digit = tdigit[0];
+	    if (!digit && !tdigit[1])
+	      nzeros += 16;
+	    else {
+	      if (!digit) {
+		nzeros += 8;
+		digit = tdigit[1];
+	      }
+	      // decompose digit
+	      PD = (UINT64) digit *0x068DB8BBull;
+	      digit_h = (UINT32) (PD >> 40);
+	      digit_low = digit - digit_h * 10000;
+
+	      if (!digit_low)
+		nzeros += 4;
+	      else
+		digit_h = digit_low;
+
+	      if (!(digit_h & 1))
+		nzeros +=
+		  3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
+				(digit_h & 7));
+	    }
+	
+	    if (nzeros) {
+	      // get P*(2^M[extra_digits])/10^extra_digits
+	      __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
+
+	      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
+	      amount = recip_scale[nzeros];
+	      __shr_128 (CQ, Qh, amount);
+	    }
+	    diff_expon += nzeros;
 
-      }
-    }
 	  }
-	if(diff_expon>=0){
-    res =
-      fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
-			       rnd_mode, pfpsf);
+	}
+      }
+      if(diff_expon>=0){
+	res =
+	  fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
+				   rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
-    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+	(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
-    BID_RETURN (res);
-	}
-  }
+        // restore the rounding mode back if it has been changed
+	if (rm_changed) fesetround(old_rm);
+	BID_RETURN (res);
+      }
+    }
 #endif
 
   if (diff_expon >= 0) {
@@ -1337,6 +1426,8 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   } else {
     // UF occurs
@@ -1353,8 +1444,9 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
-
   }
 
 }
@@ -1379,10 +1471,17 @@ int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
   digits_q, amount;
 int nzeros, i, j, k, d5, done = 0;
 unsigned rmode;
+int old_rm, rm_changed=0;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 fexcept_t binaryflags = 0;
 #endif
 
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
+
 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
 
 	// unpack arguments, check for NaN or Infinity
@@ -1401,7 +1500,9 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
       amount = recip_scale[18];
       __shr_128 (Tmp, Qh, amount);
       res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
-    BID_RETURN (res);
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
   }
   // x is Infinity?
   if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
@@ -1413,15 +1514,19 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
       res = 0x7c00000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
+    }
+    if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
+      // otherwise return +/-Inf
+      res =
+	((x.w[1] ^ y.
+	  w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
-	if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
-    // otherwise return +/-Inf
-    res =
-      ((x.w[1] ^ y.
-	w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
-    BID_RETURN (res);
-	}
   }
   // x is 0
   if (((y.w[1] & 0x7800000000000000ull) != 0x7800000000000000ull)) {
@@ -1431,6 +1536,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
 #endif
     // x=y=0, return NaN
     res = 0x7c00000000000000ull;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // return 0
@@ -1441,6 +1548,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
   else if (exponent_x < 0)
     exponent_x = 0;
   res |= (((UINT64) exponent_x) << 53);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
   }
 }
@@ -1460,12 +1569,16 @@ if (!valid_y) {
       amount = recip_scale[18];
       __shr_128 (Tmp, Qh, amount);
       res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
-    BID_RETURN (res);
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
   }
   // y is Infinity?
   if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
     // return +/-0
     res = sign_x ^ sign_y;
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // y is 0, return +/-Inf
@@ -1474,6 +1587,8 @@ if (!valid_y) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -1536,6 +1651,8 @@ if (__unsigned_compare_gt_128 (CY, CX)) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
 
@@ -1594,100 +1711,102 @@ if (!done) {
 #ifndef LEAVE_TRAILING_ZEROS
     // check whether result is exact
   {
-	  if(!done) {
-    // check whether CX, CY are short
-    if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
-      i = (int) CY.w[0] - 1;
-      j = (int) CX.w[0] - 1;
-      // difference in powers of 2 factors for Y and X
-      nzeros = ed2 - factors[i][0] + factors[j][0];
-      // difference in powers of 5 factors
-      d5 = ed2 - factors[i][1] + factors[j][1];
-      if (d5 < nzeros)
-	nzeros = d5;
-      // get P*(2^M[extra_digits])/10^extra_digits
-      __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
-      //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
-
-      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
-      amount = recip_scale[nzeros];
-      __shr_128_long (CQ, Qh, amount);
-
-      diff_expon += nzeros;
-    } else {
-      // decompose Q as Qh*10^17 + Ql
-      //T128 = reciprocals10_128[17];
-      Q_low = CQ.w[0];
-
-      {
-	tdigit[0] = Q_low & 0x3ffffff;
-	tdigit[1] = 0;
-	QX = Q_low >> 26;
-	QX32 = QX;
-	nzeros = 0;
-
-	for (j = 0; QX32; j++, QX32 >>= 7) {
-	  k = (QX32 & 127);
-	  tdigit[0] += convert_table[j][k][0];
-	  tdigit[1] += convert_table[j][k][1];
-	  if (tdigit[0] >= 100000000) {
-	    tdigit[0] -= 100000000;
-	    tdigit[1]++;
+    if(!done) {
+      // check whether CX, CY are short
+      if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
+	i = (int) CY.w[0] - 1;
+	j = (int) CX.w[0] - 1;
+	// difference in powers of 2 factors for Y and X
+	nzeros = ed2 - factors[i][0] + factors[j][0];
+	// difference in powers of 5 factors
+	d5 = ed2 - factors[i][1] + factors[j][1];
+	if (d5 < nzeros)
+	  nzeros = d5;
+	// get P*(2^M[extra_digits])/10^extra_digits
+	__mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
+	//__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
+	
+	// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
+	amount = recip_scale[nzeros];
+	__shr_128_long (CQ, Qh, amount);
+	
+	diff_expon += nzeros;
+      } else {
+	// decompose Q as Qh*10^17 + Ql
+	//T128 = reciprocals10_128[17];
+	Q_low = CQ.w[0];
+      
+	{
+	  tdigit[0] = Q_low & 0x3ffffff;
+	  tdigit[1] = 0;
+	  QX = Q_low >> 26;
+	  QX32 = QX;
+	  nzeros = 0;
+
+	  for (j = 0; QX32; j++, QX32 >>= 7) {
+	    k = (QX32 & 127);
+	    tdigit[0] += convert_table[j][k][0];
+	    tdigit[1] += convert_table[j][k][1];
+	    if (tdigit[0] >= 100000000) {
+	      tdigit[0] -= 100000000;
+	      tdigit[1]++;
+	    }
 	  }
-	}
 
-	if (tdigit[1] >= 100000000) {
-	  tdigit[1] -= 100000000;
-	  if (tdigit[1] >= 100000000)
+	  if (tdigit[1] >= 100000000) {
 	    tdigit[1] -= 100000000;
-	}
-
-	digit = tdigit[0];
-	if (!digit && !tdigit[1])
-	  nzeros += 16;
-	else {
-	  if (!digit) {
-	    nzeros += 8;
-	    digit = tdigit[1];
+	    if (tdigit[1] >= 100000000)
+	      tdigit[1] -= 100000000;
 	  }
-	  // decompose digit
-	  PD = (UINT64) digit *0x068DB8BBull;
-	  digit_h = (UINT32) (PD >> 40);
-	  digit_low = digit - digit_h * 10000;
 
-	  if (!digit_low)
-	    nzeros += 4;
-	  else
-	    digit_h = digit_low;
+	  digit = tdigit[0];
+	  if (!digit && !tdigit[1])
+	    nzeros += 16;
+	  else {
+	    if (!digit) {
+	      nzeros += 8;
+	      digit = tdigit[1];
+	    }
+	    // decompose digit
+	    PD = (UINT64) digit *0x068DB8BBull;
+	    digit_h = (UINT32) (PD >> 40);
+	    digit_low = digit - digit_h * 10000;
+
+	    if (!digit_low)
+	      nzeros += 4;
+	    else
+	      digit_h = digit_low;
+
+	    if (!(digit_h & 1))
+	      nzeros +=
+		3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
+			      (digit_h & 7));
+	  }
 
-	  if (!(digit_h & 1))
-	    nzeros +=
-	      3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
-			    (digit_h & 7));
-	}
+	  if (nzeros) {
+	    // get P*(2^M[extra_digits])/10^extra_digits
+	    __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
 
-	if (nzeros) {
-	  // get P*(2^M[extra_digits])/10^extra_digits
-	  __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
+	    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
+	    amount = recip_scale[nzeros];
+	    __shr_128 (CQ, Qh, amount);
+	  }
+	  diff_expon += nzeros;
 
-	  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
-	  amount = recip_scale[nzeros];
-	  __shr_128 (CQ, Qh, amount);
 	}
-	diff_expon += nzeros;
-
       }
     }
-	  }
-	if(diff_expon>=0){
-    res =
-      fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
-			       rnd_mode, pfpsf);
+    if(diff_expon>=0){
+      res =
+	fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
+				 rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
-    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+      (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
-    BID_RETURN (res);
-	}
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
+      BID_RETURN (res);
+    }
   }
 #endif
 
@@ -1772,6 +1891,8 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   } else {
     // UF occurs
@@ -1788,8 +1909,9 @@ if (!done) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
-
   }
 
 }
diff --git a/libgcc/config/libbid/bid64_sqrt.c b/libgcc/config/libbid/bid64_sqrt.c
index 29f4cf1..6309c42 100644..100755
--- a/libgcc/config/libbid/bid64_sqrt.c
+++ b/libgcc/config/libbid/bid64_sqrt.c
@@ -41,9 +41,8 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 
 #include "bid_internal.h"
 #include "bid_sqrt_macros.h"
-#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #include <fenv.h>
-
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
 #endif
 
@@ -73,6 +72,8 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
   int exponent_x, exponent_q, bin_expon_cx;
   int digits_x;
   int scale;
+  int old_rm, rm_changed=0;
+
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   fexcept_t binaryflags = 0;
 #endif
@@ -84,6 +85,12 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
   x = *px;
 #endif
 
+  // Set it to round-to-nearest (if different)
+  if ((old_rm=fegetround()) != FE_TONEAREST) {
+    rm_changed=1;
+    fesetround(FE_TONEAREST);
+  }
+  
   // unpack arguments, check for NaN or Infinity
   if (!unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x)) {
     // x is Inf. or NaN or 0
@@ -100,11 +107,15 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
       if ((x & SNAN_MASK64) == SNAN_MASK64)	// sNaN
 	__set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res & QUIET_MASK64);
     }
     // x is 0
     exponent_x = (exponent_x + DECIMAL_EXPONENT_BIAS) >> 1;
     res = sign_x | (((UINT64) exponent_x) << 53);
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // x<0?
@@ -113,6 +124,8 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
 #ifdef SET_STATUS_FLAGS
     __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -141,6 +154,8 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
     (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // if exponent is odd, scale coefficient by 10
@@ -206,6 +221,8 @@ bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
   (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 
@@ -223,6 +240,13 @@ int digits, scale, exponent_q = 0, exact = 1, amount, extra_digits;
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 fexcept_t binaryflags = 0;
 #endif
+int old_rm, rm_changed=0;
+
+// Set it to round-to-nearest (if different)
+if ((old_rm=fegetround()) != FE_TONEAREST) {
+  rm_changed=1;
+  fesetround(FE_TONEAREST);
+}
 
 	// unpack arguments, check for NaN or Infinity
 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
@@ -240,6 +264,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
     amount = recip_scale[18];
     __shr_128 (Tmp, Qh, amount);
     res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // x is Infinity?
@@ -251,6 +277,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
     }
+    // restore the rounding mode back if it has been changed
+    if (rm_changed) fesetround(old_rm);
     BID_RETURN (res);
   }
   // x is 0 otherwise
@@ -264,6 +292,8 @@ if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
     exponent_x = DECIMAL_MAX_EXPON_64;
   //res= sign_x | (((UINT64)exponent_x)<<53);
   res = get_BID64 (sign_x, exponent_x, 0, rnd_mode, pfpsf);
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 if (sign_x) {
@@ -271,6 +301,8 @@ if (sign_x) {
 #ifdef SET_STATUS_FLAGS
   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 #endif
+  // restore the rounding mode back if it has been changed
+  if (rm_changed) fesetround(old_rm);
   BID_RETURN (res);
 }
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
@@ -312,6 +344,8 @@ if (CS.w[0] < 10000000000000000ull) {
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
       (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+      // restore the rounding mode back if it has been changed
+      if (rm_changed) fesetround(old_rm);
       BID_RETURN (res);
     }
   }
@@ -546,7 +580,8 @@ res = get_BID64 (0, exponent_q, CS.w[0], rnd_mode, pfpsf);
 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
 #endif
+// restore the rounding mode back if it has been changed
+if (rm_changed) fesetround(old_rm);
 BID_RETURN (res);
 
-
 }