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

diff --git a/lib/ldexpl.c b/lib/ldexpl.c
index a4cdae19e..5eb9ff3ca 100644
--- a/lib/ldexpl.c
+++ b/lib/ldexpl.c
@@ -1,7 +1,7 @@
 /* Emulation for ldexpl.
    Contributed by Paolo Bonzini
 
-   Copyright 2002, 2003 Free Software Foundation, Inc.
+   Copyright 2002, 2003, 2007 Free Software Foundation, Inc.
 
    This file is part of gnulib.
 
@@ -19,32 +19,50 @@
    with this program; if not, write to the Free Software Foundation,
    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
 
-#include <float.h>
+#include <config.h>
+
+/* Specification.  */
 #include <math.h>
 
-#include "mathl.h"
+#include <float.h>
+#include "fpucw.h"
+#include "isnanl.h"
 
 long double
 ldexpl(long double x, int exp)
 {
   long double factor;
   int bit;
+  DECL_LONG_DOUBLE_ROUNDING
 
-  /* Check for zero, nan and infinity. */
-  if (x != 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;
 }