X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Fldexpl.c;h=5eb9ff3ca18cff5a6290ef293881b79fab098fd9;hb=dd6ceef65e424d30db8706273f88c20b015012ce;hp=f85847b24e9d2e7c119a4344806b0696dd426971;hpb=55eed6bbc87934f2db04981e6f0f1773ff77974f;p=gnulib.git

diff --git a/lib/ldexpl.c b/lib/ldexpl.c
index f85847b24..5eb9ff3ca 100644
--- a/lib/ldexpl.c
+++ b/lib/ldexpl.c
@@ -25,6 +25,7 @@
 #include <math.h>
 
 #include <float.h>
+#include "fpucw.h"
 #include "isnanl.h"
 
 long double
@@ -32,22 +33,36 @@ ldexpl(long double x, int exp)
 {
   long double factor;
   int bit;
+  DECL_LONG_DOUBLE_ROUNDING
 
-  /* Check for zero, nan and infinity. */
-  if (isnanl (x) || x + x == x)
-    return x;
+  BEGIN_LONG_DOUBLE_ROUNDING ();
 
-  if (exp < 0)
+  /* Check for zero, nan and infinity. */
+  if (!(isnanl (x) || x + x == x))
     {
-      exp = -exp;
-      factor = 0.5L;
+      if (exp < 0)
+	{
+	  exp = -exp;
+	  factor = 0.5L;
+	}
+      else
+	factor = 2.0L;
+
+      if (exp > 0)
+	for (bit = 1;;)
+	  {
+	    /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i,
+	       and bit <= exp.  */
+	    if (exp & bit)
+	      x *= factor;
+	    bit <<= 1;
+	    if (bit > exp)
+	      break;
+	    factor = factor * factor;
+	  }
     }
-  else
-    factor = 2.0L;
 
-  for (bit = 1; bit <= exp; bit <<= 1, factor *= factor)
-    if (exp & bit)
-      x *= factor;
+  END_LONG_DOUBLE_ROUNDING ();
 
   return x;
 }