author roland Tue, 31 Mar 1992 01:42:50 +0000 (01:42 +0000) committer roland Tue, 31 Mar 1992 01:42:50 +0000 (01:42 +0000)

index c63c050..033cbf9 100644 (file)
@@ -37,6 +37,7 @@ Cambridge, MA 02139, USA.  */

#include <ansidecl.h>
#include <math.h>

#include <ansidecl.h>
#include <math.h>
+#include <float.h>
#include "ieee754.h"

double
#include "ieee754.h"

double
@@ -44,38 +45,97 @@ DEFUN(ldexp, (x, exp),
double x AND int exp)
{
union ieee754_double u;
double x AND int exp)
{
union ieee754_double u;
+  unsigned int exponent;
+
u.d = x;
#define        x u.d

u.d = x;
#define        x u.d

-  if (__finite (x))
+  exponent = u.ieee.exponent;
+
+  /* The order of the tests is carefully chosen to handle
+     the usual case first, with no branches taken.  */
+
+  if (exponent != 0)
{
{
-      unsigned int k;
+      /* X is nonzero and not denormalized.  */

-      if (exp < -2100)
-       return 3e-39 * 3e-39;   /* Cause underflow.  */
-      else if (exp > 2100)
-       return 1.7e308 + 1.7e308; /* Cause overflow.  */
+      if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
+       {
+         /* X is finite.  When EXP < 0, overflow is actually underflow.  */

-      if (u.ieee.exponent == 0)
-       return ldexp (x * ldexp (1.0, 54), exp - 54);
+         exponent += exp;
+
+         if (exponent != 0)
+           {
+             if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
+               {
+                 /* In range.  */
+                 u.ieee.exponent = exponent;
+                 return x;
+               }
+
+             if (exp >= 0)
+             overflow:
+               {
+                 CONST int negative = u.ieee.negative;
+                 u.d = HUGE_VAL;
+                 u.ieee.negative = negative;
+                 errno = ERANGE;
+                 return u.d;
+               }
+
+             if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
+               {
+                 /* Underflow.  */
+                 CONST int negative = u.ieee.negative;
+                 u.d = 0.0;
+                 u.ieee.negative = negative;
+                 errno = ERANGE;
+                 return u.d;
+               }
+           }

-      k = u.ieee.exponent + exp;
-      if (k > 0)
-       {
-         if (k < 2047)
-           u.ieee.exponent = k;
-         else
-           return 1.7e308 + 1.7e308; /* Cause overflow.  */
-       }
-      else if (k > -54)
-       {
u.ieee.exponent = 1;
u.ieee.exponent = 1;
-         x *= ldexp (1.0, k - 1);
+         u.d *= ldexp (1.0, (int) exponent - 1);
+         if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+           /* Underflow.  */
+           errno = ERANGE;
+         return u.d;
+       }
+
+      /* X is +-infinity or NaN.  */
+      if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+       {
+         /* X is +-infinity.  */
+         if (exp >= 0)
+           goto overflow;
+         else
+           {
+             /* (infinity * number < 1).  With infinite precision,
+                (infinity / finite) would be infinity, but otherwise it's
+                safest to regard (infinity / 2) as indeterminate.  The
+                infinity might be (2 * finite).  */
+             CONST int negative = u.ieee.negative;
+             u.d = NAN;
+             u.ieee.negative = negative;
+             errno = EDOM;
+             return u.d;
+           }
}
}
-      else
-       return 3e-39 * 3e-39;   /* Cause underflow.  */
+
+      /* X is NaN.  */
+      errno = EDOM;
+      return u.d;
}

}

-  return x;
+  /* X is zero or denormalized.  */
+  if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+    /* X is +-0.0. */
+    return x;
+
+  /* X is denormalized.
+     Multiplying by 2 ** DBL_MANT_DIG normalizes it;
+     we then subtract the DBL_MANT_DIG we added to the exponent.  */
+  return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
}
}