[gfortran] Fix NEAREST constant folder

[Tobias Schlüter <tobias dot schlueter at physik dot uni-muenchen dot de>]

If the exponent shrinks during the cosntant folding of NEAREST, we get
incorrect results. E.g. NEAREST (1.,-1.) would yield too small a value.

This is fixed by the attached patch which also gives a better error for a zero
second argument to NEAREST and which fixes an ICE for a nonconstant second
argument to NEAREST.

There's one shortcoming which I enclosed into an #if 0 which I'd like some
advice on: NEAREST returns the closest machine-representable number which in
the case of zero are denormals.  range_check rejects these.  Should I simply
not range_check (and try to construct the correct value, for which my feeble
attempt didn't suffice), even though this might lead to underflows in other
places, or is there an alternative route?

Bubblestrapped and regtested on i686-linux.  Ok together with the testcase?

- Tobi

2005-04-10  Tobias Schlueter  <tobias.schlueter@physik.uni-muenchen.de>

	* simplify.c (gfc_simplify_nearest): Overhaul.

Index: simplify.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/fortran/simplify.c,v
retrieving revision 1.21
diff -c -3 -p -r1.21 simplify.c
*** simplify.c	7 Apr 2005 18:26:37 -0000	1.21
--- simplify.c	10 Apr 2005 14:17:31 -0000
*************** gfc_expr *
*** 2311,2374 ****
  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.  */
! 
!   if (x->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);
  
!   val = mpfr_get_d (x->value.real, GFC_RND_MODE);
!   p = gfc_real_kinds[k].digits;
  
!   eps = 1.;
!   for (i = 1; i < p; ++i)
      {
!       eps = eps / 2.;
      }
  
!   /* 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)
      {
!       if (match_float)
! 	{
! 	  rval = (float) val;
! 	  rval = rval + eps;
! 	  mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
! 	}
        else
! 	{
! 	  val = val + eps;
! 	  mpfr_set_d (result->value.real, val, GFC_RND_MODE);
! 	}
      }
!   else if (mpfr_cmp_ui (s->value.real, 0) < 0)
      {
!       if (match_float)
  	{
! 	  rval = (float) val;
! 	  rval = rval - eps;
! 	  mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
  	}
        else
  	{
! 	  val = val - eps;
! 	  mpfr_set_d (result->value.real, val, GFC_RND_MODE);
  	}
!     }
!   else
!     {
!       gfc_error ("Invalid second argument of NEAREST at %L", &s->where);
!       gfc_free (result);
!       return &gfc_bad_expr;
      }
  
    return range_check (result, "NEAREST");
--- 2311,2392 ----
  gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
  {
    gfc_expr *result;
!   mpfr_t tmp;
!   int direction, sgn;
  
!   if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
      return NULL;
  
!   gfc_set_model_kind (x->ts.kind);
!   result = gfc_copy_expr (x);
  
!   direction = mpfr_sgn (s->value.real);
  
!   if (direction == 0)
      {
!       gfc_error ("Second argument of NEAREST at %L may not be zero",
! 		 &s->where);
!       gfc_free (result);
!       return &gfc_bad_expr;
      }
  
!   /* 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)
      {
!       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
! 	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 (sgn < 0)
  	{
! 	  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
  	{
! 	  /* 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);
! 
! 	  /* Are we back to where we started?  */
! 	  if (mpfr_cmp (tmp, result->value.real) == 0)
! 	    mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
! 
! 	  mpfr_clear (tmp);
  	}
! 
!       if (sgn < 0)
! 	mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
      }
  
    return range_check (result, "NEAREST");

! { 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 ()
!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