OpenE2K/gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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

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

From-SVN: r85724
This commit is contained in:
Roger Sayle 2004-08-09 21:09:41 +00:00
parent 9d8646d7b0
commit 067a5735c5
2 changed files with 49 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
libgfortran/ChangeLog
View File

@ -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:

85
libgfortran/intrinsics/c99_functions.c
View File

@ -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