diff options
author | Roger Sayle <sayle@gcc.gnu.org> | 2004-08-09 21:09:41 +0000 |
---|---|---|
committer | Roger Sayle <sayle@gcc.gnu.org> | 2004-08-09 21:09:41 +0000 |
commit | 067a5735c5236147708d701c247a8bee38c06a8b (patch) | |
tree | f2079c84cbe12921ff8851d7d5ad3c28a97c5a41 /libgfortran | |
parent | 9d8646d7b098cd131bd9263ec98d3c03bdf4f66b (diff) | |
download | gcc-067a5735c5236147708d701c247a8bee38c06a8b.zip gcc-067a5735c5236147708d701c247a8bee38c06a8b.tar.gz gcc-067a5735c5236147708d701c247a8bee38c06a8b.tar.bz2 |
c99_functions.c (nextafterf): New implementation that works correctly with denormalized numbers.
* intrinsics/c99_functions.c (nextafterf): New implementation that
works correctly with denormalized numbers.
From-SVN: r85724
Diffstat (limited to 'libgfortran')
-rw-r--r-- | libgfortran/ChangeLog | 8 | ||||
-rw-r--r-- | libgfortran/intrinsics/c99_functions.c | 85 |
2 files changed, 49 insertions, 44 deletions
diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index cc27e33..de54f0f 100644 --- 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 9ae84ed..eb805c3 100644 --- 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 |