Re: [PATCH/gfortran] Fix rounding in constant folding

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

Tobias Schlüter wrote:
> gives the following output:
>   0.5000000               1
>   0.4999999               0
>   0.5000000               0
> So the result for NEAREST as evaluated by the constant folders is too small.
> It looks like mpfr provides what is needed, so I'll cook up a patch.

Unfortunately, I couldn't get it to work, neither with mpfr_{add|sub}_one_ulp,
nor with the hand-written equivalent (?) below.  The following testcase
   print *, 0.49999997, 0.49999997 < 0.5, nint (0.49999997)
   print *, nearest (0.5, -1.), nearest (0.5, -1.) < 0.5, &
        nint(nearest (0.5, -1.))
   x = 0.5
   print *, nearest (x, -1.), nearest (x, -1.) < 0.5, &
        nint(nearest(x, -1.))
   end
produces this result with both approaches:
   [tobi@marktplatz tests]$ ./a.out
     0.5000000     T           1
     0.4999999     T           0
     0.5000000     T           0
   [tobi@marktplatz tests]$
the second to last line should match the last line, so something's wrong here.

I'm going to bed now, maybe I'll see clearer tomorrow.  Is there something
obviously wrong?

- Tobi

ps I know that I'm introducing a memory leak, this is by no means meant for
inclusion.

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  9 Apr 2005 23:38:18 -0000
*************** gfc_expr *
*** 2311,2372 ****
  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;
      }
--- 2311,2357 ----
  gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
  {
    gfc_expr *result;
!   mpfr_t ulp;
!   int direction;

!   if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
      return NULL;

!   /* TODO: Use mpfr_nextabove and mpfr_nextbelow once we move to a
!      newer version of mpfr.  */

    result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+   mpfr_init (ulp);
+   mpfr_set (ulp,
+           gfc_real_kinds[gfc_validate_kind (BT_REAL, x->ts.kind, 0)].epsilon,
+           GFC_RND_MODE);

!   direction = mpfr_cmp_ui (s->value.real, 0);
!   if (direction > 0)
      {
!       /* mpfr_add_one_ulp (result->value.real);  (after initializing
!        it to x->value.real, of course) didn't do the trick.  */
!       do
        {
!         mpfr_add (result->value.real, x->value.real, ulp, GFC_RND_MODE);
!         mpfr_mul_ui (ulp, ulp, 2, GFC_RND_MODE);
        }
+       while (mpfr_cmp (result->value.real, x->value.real) == 0);
      }
!   else if (direction < 0)
      {
!       /* As above but with mpfr_sub_one_ulp  */
!       do
        {
!         mpfr_sub (result->value.real, x->value.real, ulp, GFC_RND_MODE);
!         mpfr_mul_ui (ulp, ulp, 2, GFC_RND_MODE);
        }
+       while (mpfr_cmp (result->value.real, x->value.real) == 0);
      }
    else
      {
!       gfc_error ("Second argument of NEAREST at %L may not be zero",
!                &s->where);
        gfc_free (result);
        return &gfc_bad_expr;
      }