simplify.c (gfc_simplify_nearest): Overhaul.

fortran/
* simplify.c (gfc_simplify_nearest): Overhaul.

testsuite/
* gfortran.dg/fold_nearest.f90: New test.

From-SVN: r97987

gcc/fortran/ChangeLog

+2005-04-11  Tobias Schl"uter  <tobias.schlueter@physik.uni-muenchen.de>
+
+	* simplify.c (gfc_simplify_nearest): Overhaul.
+
 2005-04-10  Kazu Hirata  <kazu@cs.umass.edu>
 
 	* interface.c: Fix a comment typo.

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");

gcc/testsuite/ChangeLog

+2005-04-11  Tobias Schl"uter  <tobias.schlueter@physik.uni-muenchen.de>
+
+	* gfortran.dg/fold_nearest.f90: New test.
+
 2005-04-11  Andrew Pinski  <pinskia@physics.uc.edu>
 
 	* gcc.dg/tree-ssa/alias-1.c: New test.

gcc/testsuite/gfortran.dg/fold_nearest.f90

+! { 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