From a5ac7fc326df28a68d0147a5203839762d5de30c Mon Sep 17 00:00:00 2001 From: tobi Date: Mon, 11 Apr 2005 21:48:27 +0000 Subject: [PATCH] fortran/ * simplify.c (gfc_simplify_nearest): Overhaul. testsuite/ * gfortran.dg/fold_nearest.f90: New test. git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@97987 138bc75d-0d04-0410-961f-82ee72b054a4 --- gcc/fortran/ChangeLog | 4 ++ gcc/fortran/simplify.c | 100 +++++++++++++++++------------ gcc/testsuite/ChangeLog | 4 ++ gcc/testsuite/gfortran.dg/fold_nearest.f90 | 27 ++++++++ 4 files changed, 94 insertions(+), 41 deletions(-) create mode 100644 gcc/testsuite/gfortran.dg/fold_nearest.f90 diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index 04cab329d18..5df549eadd1 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,7 @@ +2005-04-11 Tobias Schl"uter + + * simplify.c (gfc_simplify_nearest): Overhaul. + 2005-04-10 Kazu Hirata * interface.c: Fix a comment typo. diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c index e4b891697bb..1ca5b52bdfb 100644 --- a/gcc/fortran/simplify.c +++ b/gcc/fortran/simplify.c @@ -2263,64 +2263,82 @@ gfc_expr * gfc_simplify_nearest (gfc_expr * x, gfc_expr * s) { gfc_expr *result; - float rval; - double val, eps; - int p, i, k, match_float; - - /* FIXME: This implementation is dopey and probably not quite right, - but it's a start. */ + mpfr_t tmp; + int direction, sgn; - if (x->expr_type != EXPR_CONSTANT) + if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT) return NULL; - k = gfc_validate_kind (x->ts.type, x->ts.kind, false); - - result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where); + gfc_set_model_kind (x->ts.kind); + result = gfc_copy_expr (x); - val = mpfr_get_d (x->value.real, GFC_RND_MODE); - p = gfc_real_kinds[k].digits; + direction = mpfr_sgn (s->value.real); - eps = 1.; - for (i = 1; i < p; ++i) + if (direction == 0) { - eps = eps / 2.; + gfc_error ("Second argument of NEAREST at %L may not be zero", + &s->where); + gfc_free (result); + return &gfc_bad_expr; } - /* TODO we should make sure that 'float' matches kind 4 */ - match_float = gfc_real_kinds[k].kind == 4; - if (mpfr_cmp_ui (s->value.real, 0) > 0) + /* TODO: Use mpfr_nextabove and mpfr_nextbelow once we move to a + newer version of mpfr. */ + + sgn = mpfr_sgn (x->value.real); + + if (sgn == 0) { - if (match_float) - { - rval = (float) val; - rval = rval + eps; - mpfr_set_d (result->value.real, rval, GFC_RND_MODE); - } + int k = gfc_validate_kind (BT_REAL, x->ts.kind, 0); + + if (direction > 0) + mpfr_add (result->value.real, + x->value.real, gfc_real_kinds[k].tiny, GFC_RND_MODE); else - { - val = val + eps; - mpfr_set_d (result->value.real, val, GFC_RND_MODE); - } + mpfr_sub (result->value.real, + x->value.real, gfc_real_kinds[k].tiny, GFC_RND_MODE); + +#if 0 + /* FIXME: This gives an arithmetic error because we compare + against tiny when range-checking. Also, it doesn't give the + right value. */ + /* TINY is the smallest model number, we want the smallest + machine representable number. Therefore we have to shift the + value to the right by the number of digits - 1. */ + mpfr_div_2ui (result->value.real, result->value.real, + gfc_real_kinds[k].precision - 1, GFC_RND_MODE); +#endif } - else if (mpfr_cmp_ui (s->value.real, 0) < 0) + else { - if (match_float) + if (sgn < 0) { - rval = (float) val; - rval = rval - eps; - mpfr_set_d (result->value.real, rval, GFC_RND_MODE); + direction = -direction; + mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE); } + + if (direction > 0) + mpfr_add_one_ulp (result->value.real, GFC_RND_MODE); else { - val = val - eps; - mpfr_set_d (result->value.real, val, GFC_RND_MODE); + /* In this case the exponent can shrink, which makes us skip + over one number because we substract one ulp with the + larger exponent. Thus we need to compensate for this. */ + mpfr_init_set (tmp, result->value.real, GFC_RND_MODE); + + mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE); + mpfr_add_one_ulp (result->value.real, GFC_RND_MODE); + + /* If we're back to where we started, the spacing is one + ulp, and we get the correct result by subtracting. */ + if (mpfr_cmp (tmp, result->value.real) == 0) + mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE); + + mpfr_clear (tmp); } - } - else - { - gfc_error ("Invalid second argument of NEAREST at %L", &s->where); - gfc_free (result); - return &gfc_bad_expr; + + if (sgn < 0) + mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE); } return range_check (result, "NEAREST"); diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index ff03dbee86b..2b10695cfe6 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,7 @@ +2005-04-11 Tobias Schl"uter + + * gfortran.dg/fold_nearest.f90: New test. + 2005-04-11 Andrew Pinski * gcc.dg/tree-ssa/alias-1.c: New test. diff --git a/gcc/testsuite/gfortran.dg/fold_nearest.f90 b/gcc/testsuite/gfortran.dg/fold_nearest.f90 new file mode 100644 index 00000000000..743e2023ab6 --- /dev/null +++ b/gcc/testsuite/gfortran.dg/fold_nearest.f90 @@ -0,0 +1,27 @@ +! { dg-do run } +! Tests for the constant folding of the NEAREST intrinsic +! We compare against the results of the runtime implementation, +! thereby making sure that they remain consistent +REAL, PARAMETER :: x(10) = (/ 1., 0.49999997, 0.5, 8388609.0, -1., & + -0.49999997, -0.5, -8388609.0, & + 0., 0. /), & + dir(10) = (/ -1., +1., -1., -1., +1., & + -1., +1., +1., & + +1.,-1./) +REAL :: a(10) + +a = x +if (nearest (x(1), dir(1)) /= nearest (a(1), dir(1))) call abort () +if (nearest (x(2), dir(2)) /= nearest (a(2), dir(2))) call abort () +if (nearest (x(3), dir(3)) /= nearest (a(3), dir(3))) call abort () +if (nearest (x(4), dir(4)) /= nearest (a(4), dir(4))) call abort () +if (nearest (x(5), dir(5)) /= nearest (a(5), dir(5))) call abort () +if (nearest (x(6), dir(6)) /= nearest (a(6), dir(6))) call abort () +if (nearest (x(7), dir(7)) /= nearest (a(7), dir(7))) call abort () +if (nearest (x(8), dir(8)) /= nearest (a(8), dir(8))) call abort () +! These last two tests are commented out because mpfr provides no support +! for denormals, and therefore we get TINY instead of the correct result. +!if (nearest (x(9), dir(9)) /= nearest (a(9), dir(9))) call abort () +!if (nearest (x(10), dir(10)) /= nearest (a(10), dir(10))) call abort () + +end -- 2.11.0