re PR fortran/34192 (NEAREST can return wrong numbers)

2007-11-23  Tobias Burnus  <burnus@net-b.de>
            Steven G. Kargl  <kargl@gcc.gnu.org>

        PR fortran/34192
        * simplify.c (gfc_simplify_nearest): Fix NEAREST for
        subnormal numbers.

2007-11-23  Tobias Burnus  <burnus@net-b.de>

        PR fortran/34192
        * gfortran.dg/nearest_2.f90: New.


Co-Authored-By: Steven G. Kargl <kargl@gcc.gnu.org>

From-SVN: r130383

4 files changed, 194 insertions, 8 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 85f2a5a..93f775e 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,10 @@
+2007-11-23  Tobias Burnus  <burnus@net-b.de>
+	    Steven G. Kargl  <kargl@gcc.gnu.org>
+
+	PR fortran/34192
+	* simplify.c (gfc_simplify_nearest): Fix NEAREST for
+	subnormal numbers.
+
 2007-11-23  Aldy Hernandez  <aldyh@redhat.com>
 
 	* trans-expr.c (gfc_trans_string_copy): Use "void *" when building a
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index cdf1118..687e87f 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -2691,8 +2691,8 @@ gfc_expr *
 gfc_simplify_nearest (gfc_expr *x, gfc_expr *s)
 {
   gfc_expr *result;
-  mpfr_t tmp;
-  int sgn;
+  mp_exp_t emin, emax;
+  int kind;
 
   if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -2707,13 +2707,39 @@ gfc_simplify_nearest (gfc_expr *x, gfc_expr *s)
   gfc_set_model_kind (x->ts.kind);
   result = gfc_copy_expr (x);
 
-  sgn = mpfr_sgn (s->value.real); 
-  mpfr_init (tmp);
-  mpfr_set_inf (tmp, sgn);
-  mpfr_nexttoward (result->value.real, tmp);
-  mpfr_clear (tmp);
+  /* Save current values of emin and emax.  */
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
+  /* Set emin and emax for the current model number.  */
+  kind = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
+  mpfr_set_emin ((mp_exp_t) gfc_real_kinds[kind].min_exponent -
+		mpfr_get_prec(result->value.real) + 1);
+  mpfr_set_emax ((mp_exp_t) gfc_real_kinds[kind].max_exponent - 1);
+
+  if (mpfr_sgn (s->value.real) > 0)
+    {
+      mpfr_nextabove (result->value.real);
+      mpfr_subnormalize (result->value.real, 0, GMP_RNDU);
+    }
+  else
+    {
+      mpfr_nextbelow (result->value.real);
+      mpfr_subnormalize (result->value.real, 0, GMP_RNDD);
+    }
+
+  mpfr_set_emin (emin);
+  mpfr_set_emax (emax);
 
-  return range_check (result, "NEAREST");
+  /* Only NaN can occur. Do not use range check as it gives an
+     error for denormal numbers.  */
+  if (mpfr_nan_p (result->value.real) && gfc_option.flag_range_check)
+    {
+      gfc_error ("Result of NEAREST is NaN at %L", &result->where);
+      return &gfc_bad_expr;
+    }
+
+  return result;
 }
 
 
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 6bd778a..d87601f 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,9 @@
+2007-11-23  Tobias Burnus  <burnus@net-b.de>
+
+	PR fortran/34192
+	* gfortran.dg/nearest_2.f90: New.
+
+
 2007-11-23  Jakub Jelinek  <jakub@redhat.com>
 
 	PR c++/30293
diff --git a/gcc/testsuite/gfortran.dg/nearest_2.f90 b/gcc/testsuite/gfortran.dg/nearest_2.f90
new file mode 100644
index 0000000..4bdad31
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/nearest_2.f90
@@ -0,0 +1,147 @@
+! { dg-do run }
+!
+! PR fortran/34192
+!
+! Test compile-time implementation of NEAREST
+!
+program test
+  implicit none
+
+! Single precision
+
+  ! 0+ > 0
+  if (nearest(0.0, 1.0) &
+      <= 0.0) &
+    call abort()
+  ! 0++ > 0+
+  if (nearest(nearest(0.0, 1.0), 1.0) &
+      <= nearest(0.0, 1.0)) &
+    call abort()
+  ! 0+++ > 0++
+  if (nearest(nearest(nearest(0.0, 1.0), 1.0), 1.0) &
+      <= nearest(nearest(0.0, 1.0), 1.0)) &
+    call abort()
+  ! 0+- = 0
+  if (nearest(nearest(0.0, 1.0), -1.0) &
+      /= 0.0) &
+    call abort()
+  ! 0++- = 0+
+  if (nearest(nearest(nearest(0.0, 1.0), 1.0), -1.0) &
+      /= nearest(0.0, 1.0)) &
+    call abort()
+  ! 0++-- = 0
+  if (nearest(nearest(nearest(nearest(0.0, 1.0), 1.0), -1.0), -1.0) &
+      /= 0.0) &
+    call abort()
+
+  ! 0- < 0
+  if (nearest(0.0, -1.0) &
+      >= 0.0) &
+    call abort()
+  ! 0-- < 0+
+  if (nearest(nearest(0.0, -1.0), -1.0) &
+      >= nearest(0.0, -1.0)) &
+    call abort()
+  ! 0--- < 0--
+  if (nearest(nearest(nearest(0.0, -1.0), -1.0), -1.0) &
+      >= nearest(nearest(0.0, -1.0), -1.0)) &
+    call abort()
+  ! 0-+ = 0
+  if (nearest(nearest(0.0, -1.0), 1.0) &
+      /= 0.0) &
+    call abort()
+  ! 0--+ = 0-
+  if (nearest(nearest(nearest(0.0, -1.0), -1.0), 1.0) &
+      /= nearest(0.0, -1.0)) &
+    call abort()
+  ! 0--++ = 0
+  if (nearest(nearest(nearest(nearest(0.0, -1.0), -1.0), 1.0), 1.0) &
+      /= 0.0) &
+    call abort()
+
+  ! 42++ > 42+
+  if (nearest(nearest(42.0, 1.0), 1.0) &
+      <= nearest(42.0, 1.0)) &
+    call abort()
+  ! 42-- < 42-
+  if (nearest(nearest(42.0, -1.0), -1.0) &
+      >= nearest(42.0, -1.0)) &
+    call abort()
+  ! 42-+ = 42
+  if (nearest(nearest(42.0, -1.0), 1.0) &
+      /= 42.0) &
+    call abort()
+  ! 42+- = 42
+  if (nearest(nearest(42.0, 1.0), -1.0) &
+      /= 42.0) &
+    call abort()
+
+! Double precision
+
+  ! 0+ > 0
+  if (nearest(0.0d0, 1.0) &
+      <= 0.0d0) &
+    call abort()
+  ! 0++ > 0+
+  if (nearest(nearest(0.0d0, 1.0), 1.0) &
+      <= nearest(0.0d0, 1.0)) &
+    call abort()
+  ! 0+++ > 0++
+  if (nearest(nearest(nearest(0.0d0, 1.0), 1.0), 1.0) &
+      <= nearest(nearest(0.0d0, 1.0), 1.0)) &
+    call abort()
+  ! 0+- = 0
+  if (nearest(nearest(0.0d0, 1.0), -1.0) &
+      /= 0.0d0) &
+    call abort()
+  ! 0++- = 0+
+  if (nearest(nearest(nearest(0.0d0, 1.0), 1.0), -1.0) &
+      /= nearest(0.0d0, 1.0)) &
+    call abort()
+  ! 0++-- = 0
+  if (nearest(nearest(nearest(nearest(0.0d0, 1.0), 1.0), -1.0), -1.0) &
+      /= 0.0d0) &
+    call abort()
+
+  ! 0- < 0
+  if (nearest(0.0d0, -1.0) &
+      >= 0.0d0) &
+    call abort()
+  ! 0-- < 0+
+  if (nearest(nearest(0.0d0, -1.0), -1.0) &
+      >= nearest(0.0d0, -1.0)) &
+    call abort()
+  ! 0--- < 0--
+  if (nearest(nearest(nearest(0.0d0, -1.0), -1.0), -1.0) &
+      >= nearest(nearest(0.0d0, -1.0), -1.0)) &
+    call abort()
+  ! 0-+ = 0
+  if (nearest(nearest(0.0d0, -1.0), 1.0) &
+      /= 0.0d0) &
+    call abort()
+  ! 0--+ = 0-
+  if (nearest(nearest(nearest(0.0d0, -1.0), -1.0), 1.0) &
+      /= nearest(0.0d0, -1.0)) &
+    call abort()
+  ! 0--++ = 0
+  if (nearest(nearest(nearest(nearest(0.0d0, -1.0), -1.0), 1.0), 1.0) &
+      /= 0.0d0) &
+    call abort()
+
+  ! 42++ > 42+
+  if (nearest(nearest(42.0d0, 1.0), 1.0) &
+      <= nearest(42.0d0, 1.0)) &
+    call abort()
+  ! 42-- < 42-
+  if (nearest(nearest(42.0d0, -1.0), -1.0) &
+      >= nearest(42.0d0, -1.0)) &
+    call abort()
+  ! 42-+ = 42
+  if (nearest(nearest(42.0d0, -1.0), 1.0) &
+      /= 42.0d0) &
+    call abort()
+  ! 42+- = 42
+  if (nearest(nearest(42.0d0, 1.0), -1.0) &
+      /= 42.0d0) &
+    call abort()
+end program test