c99_functions.c (nextafterf): New implementation that works correctly with denormaliz...

	* intrinsics/c99_functions.c (nextafterf): New implementation that
	works correctly with denormalized numbers.

From-SVN: r85724

diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog
index cc27e33325c05cf97103f5c7bcd49200f4d345b6..de54f0f21155ac885a6f3af6fad2b1724e675230 100644 (file)
--- a/libgfortran/ChangeLog
+++ b/libgfortran/ChangeLog
@@ -1,4 +1,10 @@
-2004-09-09  Victor Leikehman  <lei@il.ibm.com>
+2004-08-09  Richard Henderson  <rth@redhat.com>
+           Roger Sayle  <roger@eyesopen.com>
+
+       * intrinsics/c99_functions.c (nextafterf): New implementation that
+       works correctly with denormalized numbers.
+
+2004-08-09  Victor Leikehman  <lei@il.ibm.com>
 
        * m4/matmul.m4, m4/matmull.m4, intrinsics/eoshift0.c,
        intrinsics/eoshift2.c, intrinsics/transpose_generic.c:

diff --git a/libgfortran/intrinsics/c99_functions.c b/libgfortran/intrinsics/c99_functions.c
index 9ae84ed93550eef177d19c4b46b99e4fc7f79eea..eb805c36f088a76b31d678b2f185846f4bc9fbf1 100644 (file)
--- a/libgfortran/intrinsics/c99_functions.c
+++ b/libgfortran/intrinsics/c99_functions.c
@@ -188,61 +188,60 @@ tanhf(float x)
 #ifndef HAVE_NEXTAFTERF
 /* This is a portable implementation of nextafterf that is intended to be
    independent of the floating point format or its in memory representation.
-   This implementation skips denormalized values, for example returning
-   FLT_MIN as the next value after zero, as many target's frexpf, scalbnf
-   and ldexpf functions don't work as expected with denormalized values.  */
+   This implementation works correctly with denormalized values.  */
 float
 nextafterf(float x, float y)
 {
-  int origexp, newexp;
+  /* This variable is marked volatile to avoid excess precision problems
+     on some platforms, including IA-32.  */
+  volatile float delta;
+  float absx, denorm_min;
 
   if (isnan(x) || isnan(y))
-    return x+y;
+    return x + y;
   if (x == y)
     return x;
 
-  if (x == 0.0f)
-    return y > 0.0f ? FLT_MIN : -FLT_MIN;
+  /* absx = fabsf (x);  */
+  absx = (x < 0.0) ? -x : x;
 
-  frexpf(x, &origexp);
-  if (x >= 0.0)
-    {
-      if (y > x)
-       {
-         if (x < FLT_MIN)
-           return FLT_MIN;
-         return x + scalbnf(FLT_EPSILON, origexp-1);
-       }
-      else if (x > FLT_MIN)
-       {
-         float temp = x - scalbnf(FLT_EPSILON, origexp-1);
-         frexpf(temp, &newexp);
-         if (newexp == origexp)
-           return temp;
-         return x - scalbnf(FLT_EPSILON, origexp-2);
-       }
-      else
-       return 0.0f;
-    }
+  /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
+  if (__FLT_DENORM_MIN__ == 0.0f)
+    denorm_min = __FLT_MIN__;
+  else
+    denorm_min = __FLT_DENORM_MIN__;
+
+  if (absx < __FLT_MIN__)
+    delta = denorm_min;
   else
     {
-      if (y < x)
-       {
-         if (x > -FLT_MIN)
-           return -FLT_MIN;
-         return x - scalbnf(FLT_EPSILON, origexp-1);
-       }
-      else if (x < -FLT_MIN)
-       {
-         float temp = x + scalbnf(FLT_EPSILON, origexp-1);
-         frexpf(temp, &newexp);
-         if (newexp == origexp)
-           return temp;
-         return x + scalbnf(FLT_EPSILON, origexp-2);
-       }
-      else
-       return 0.0f;
+      float frac;
+      int exp;
+
+      /* Discard the fraction from x.  */
+      frac = frexpf (absx, &exp);
+      delta = scalbnf (0.5f, exp);
+
+      /* Scale x by the epsilon of the representation.  By rights we should
+        have been able to combine this with scalbnf, but some targets don't
+        get that correct with denormals.  */
+      delta *= __FLT_EPSILON__;
+
+      /* If we're going to be reducing the absolute value of X, and doing so
+        would reduce the exponent of X, then the delta to be applied is
+        one exponent smaller.  */
+      if (frac == 0.5f && (y < x) == (x > 0))
+       delta *= 0.5f;
+
+      /* If that underflows to zero, then we're back to the minimum.  */
+      if (delta == 0.0f)
+       delta = denorm_min;
     }
+
+  if (y < x)
+    delta = -delta;
+
+  return x + delta;
 }
 #endif