Index: libgcc/config/rs6000/ibm-ldouble.c
===================================================================
--- libgcc/config/rs6000/ibm-ldouble.c	(revision 228280)
+++ libgcc/config/rs6000/ibm-ldouble.c	(working copy)
@@ -47,6 +47,10 @@ see the files COPYING3 and COPYING.RUNTIME respect
 
 #if defined (__MACH__) || defined (__powerpc__) || defined (_AIX)
 
+typedef _Bool bool;
+#define false 0
+#define true 1
+
 #define fabs(x) __builtin_fabs(x)
 #define isless(x, y) __builtin_isless (x, y)
 #define inf() __builtin_inf()
@@ -87,6 +91,58 @@ __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
 	 ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
 #endif
 
+#define FE_TONEAREST 0
+
+/* Set the rounding mode to MODE, returning the previous mode.  */
+
+static int
+swapround (int mode)
+{
+  int old_mode;
+#ifdef __NO_FPRS__
+# ifdef _SOFT_FLOAT
+  /* Not implemented.  */
+  old_mode = FE_TONEAREST;
+# else
+  /* SPE floating point.  */
+  int spefscr;
+  asm volatile ("mfspefscr %0" : "=r" (spefscr));
+  old_mode = spefscr & 3;
+  spefscr = (spefscr & ~3) | mode;
+  asm volatile ("mtspefscr %0" : : "r" (spefscr));
+# endif
+#else
+  /* Logic as in glibc.  */
+  asm volatile ("mcrfs 7,7\n\t"
+		"mfcr %0" : "=r" (old_mode) : : "cr7");
+  old_mode &= 3;
+  if ((mode ^ old_mode) & 2)
+    {
+      if (mode & 2)
+	asm volatile ("mtfsb1 30");
+      else
+	asm volatile ("mtfsb0 30");
+    }
+  if ((mode ^ old_mode) & 1)
+    {
+      if (mode & 1)
+	asm volatile ("mtfsb1 31");
+      else
+	asm volatile ("mtfsb0 31");
+    }
+#endif
+  return old_mode;
+}
+
+/* Overflowing values in each rounding mode (FE_TONEAREST,
+   FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD).  */
+
+static const long double overflow_values[2][4] =
+  {
+    { __builtin_infl (), __LDBL_MAX__, __builtin_infl (), __LDBL_MAX__ },
+    { -__builtin_infl (), -__LDBL_MAX__, -__LDBL_MAX__, -__builtin_infl () }
+  };
+
 /* Combine two 'double' values into one 'long double' and return the result.  */
 static inline long double
 pack_ldouble (double dh, double dl)
@@ -110,40 +166,74 @@ pack_ldouble (double dh, double dl)
 long double
 __gcc_qadd (double a, double aa, double c, double cc)
 {
-  double xh, xl, z, q, zz;
+  double xh, xhs, xl, z, q, zz;
+  bool scale = false;
+  int old_mode;
 
+  if (nonfinite (a) || nonfinite (c))
+    return a + c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p1022 || fabs (c) >= 0x1p1022)
+    {
+      scale = true;
+      if (fabs (a) >= 0x1p914)
+	{
+	  a *= 0.5;
+	  aa *= 0.5;
+	}
+      if (fabs (c) >= 0x1p914)
+	{
+	  c *= 0.5;
+	  cc *= 0.5;
+	}
+    }
+
   z = a + c;
+  q = a - z;
+  zz = q + c + (a - (q + z)) + aa + cc;
 
-  if (nonfinite (z))
+  /* Keep -0 result.  */
+  if (zz == 0.0)
     {
-      if (fabs (z) != inf())
-	return z;
-      z = cc + aa + c + a;
-      if (nonfinite (z))
-	return z;
-      xh = z;  /* Will always be DBL_MAX.  */
-      zz = aa + cc;
-      if (fabs(a) > fabs(c))
-	xl = a - z + c + zz;
-      else
-	xl = c - z + a + zz;
+      double ret = scale ? 2.0 * z : z;
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (ret) == inf ())
+	    return overflow_values[ret < 0][old_mode];
+	  if (ret == 0)
+	    {
+	      asm volatile ("" : "+m" (a));
+	      ret = a + c;
+	    }
+	}
+      return ret;
     }
-  else
+  xh = z + zz;
+  xhs = xh;
+  if (scale)
+    xhs *= 2.0;
+  if (nonfinite (xhs))
     {
-      q = a - z;
-      zz = q + c + (a - (q + z)) + aa + cc;
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (xhs) == inf ())
+	    return overflow_values[xhs < 0][old_mode];
+	}
+      return xhs;
+    }
 
-      /* Keep -0 result.  */
-      if (zz == 0.0)
-	return z;
+  xl = z - xh + zz;
+  if (scale)
+    xl *= 2.0;
 
-      xh = z + zz;
-      if (nonfinite (xh))
-	return xh;
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
 
-      xl = z - xh + zz;
-    }
-  return pack_ldouble (xh, xl);
+  return pack_ldouble (xhs, xl);
 }
 
 long double
@@ -159,13 +249,41 @@ static double fmsub (double, double, double);
 long double
 __gcc_qmul (double a, double b, double c, double d)
 {
-  double xh, xl, t, tau, u, v, w;
+  double xh, xl, t, tau, u, us, v, w;
+  bool scale = false;
+  int old_mode;
+
+  if (nonfinite (a) || nonfinite (c))
+    return a * c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p512)
+    {
+      scale = true;
+      a *= 0.5;
+      b *= 0.5;
+    }
+  else if (fabs (c) >= 0x1p512)
+    {
+      scale = true;
+      c *= 0.5;
+      d *= 0.5;
+    }
   
   t = a * c;			/* Highest order double term.  */
 
   if (unlikely (t == 0)		/* Preserve -0.  */
       || nonfinite (t))
-    return t;
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (t) == inf ())
+	    return overflow_values[t < 0][old_mode];
+	}
+      return t;
+    }
 
   /* Sum terms of two highest orders. */
   
@@ -179,12 +297,29 @@ __gcc_qmul (double a, double b, double c, double d
   w = b*c;
   tau += v + w;	    /* Add in other second-order terms.	 */
   u = t + tau;
+  us = u;
+  if (scale)
+    us *= 2.0;
 
   /* Construct long double result.  */
-  if (nonfinite (u))
-    return u;
-  xh = u;
+  if (nonfinite (us))
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (us) == inf ())
+	    return overflow_values[us < 0][old_mode];
+	}
+      return us;
+    }
+  xh = us;
   xl = (t - u) + tau;
+  if (scale)
+    xl *= 2.0;
+
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
+
   return pack_ldouble (xh, xl);
 }
 
@@ -191,13 +326,41 @@ __gcc_qmul (double a, double b, double c, double d
 long double
 __gcc_qdiv (double a, double b, double c, double d)
 {
-  double xh, xl, s, sigma, t, tau, u, v, w;
+  double xh, xl, s, sigma, t, tau, u, us, v, w;
+  bool scale = false;
+  int old_mode;
+
+  if (nonfinite (a) || nonfinite (c) || c == 0)
+    return a / c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p512)
+    {
+      scale = true;
+      a *= 0.5;
+      b *= 0.5;
+    }
+  else if (fabs (c) <= 0x1p-511)
+    {
+      scale = true;
+      c *= 2.0;
+      d *= 2.0;
+    }
   
   t = a / c;                    /* highest order double term */
   
   if (unlikely (t == 0)		/* Preserve -0.  */
       || nonfinite (t))
-    return t;
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (t) == inf ())
+	    return overflow_values[t < 0][old_mode];
+	}
+      return t;
+    }
 
   /* Finite nonzero result requires corrections to the highest order
      term.  These corrections require the low part of c * t to be
@@ -224,12 +387,29 @@ __gcc_qdiv (double a, double b, double c, double d
   
   tau = ((v-sigma)+w)/c;   /* Correction to t.  */
   u = t + tau;
+  us = u;
+  if (scale)
+    us *= 2.0;
 
   /* Construct long double result.  */
-  if (nonfinite (u))
-    return u;
-  xh = u;
+  if (nonfinite (us))
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (us) == inf ())
+	    return overflow_values[us < 0][old_mode];
+	}
+      return us;
+    }
+  xh = us;
   xl = (t - u) + tau;
+  if (scale)
+    xl *= 2.0;
+
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
+
   return pack_ldouble (xh, xl);
 }
 
