[PATCH/gfortran] Fix rounding in constant folding

[Steve Kargl]

On Sat, Apr 09, 2005 at 07:45:01PM -0700, Steve Kargl wrote:
> On Sun, Apr 10, 2005 at 01:49:31AM +0200, Tobias Schl?ter wrote:
> > !   direction = mpfr_cmp_ui (s->value.real, 0);
> 
>       direction = mpfr_sgn (s->value.real);
> 
> > !   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);
> 
> How many times do you go through the loop?  I would think that
> you would want to add another ulp each time through.  Right now
> you have ulp, 2*ulp, 4*ulp, 8*ulp, etc.
> 
> BTW, the algorithm that you're using is what I would do.
> Also, note GFC_RND_MODE gives round-to-nearest behavior.
> You can set this to one of
> 
>    * `GMP_RNDN': round to nearest
>    * `GMP_RNDZ': round towards zero
>    * `GMP_RNDU': round towards plus infinity
>    * `GMP_RNDD': round towards minus infinity

I used mpfr_add_one_ulp and mpfr_sub_one_ulp, I believe 
MPFR is doing the right.  I instrumented gfc_simplify_nearest
with 

  mpfr_print_binary (x->value.real);
  printf("\n");
  mpfr_print_binary (result->value.real);
  printf("\n");

which yields 

0.100000000000000000000000[00000000]E0
0.111111111111111111111110[00000000]E-1
0.100000000000000000000000[00000000]E0
0.111111111111111111111110[00000000]E-1
0.100000000000000000000000[00000000]E0
0.111111111111111111111110[00000000]E-1
 
Your program is giving

kargl[238] ./z
  0.5000000     T           0
  0.4999999     T           0
  0.5000000     T           0

so there may be a rounding problem with IO in the runtime library.
I've attached my diff for your amusement.

-- 
Steve
-------------- next part --------------
Index: simplify.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/fortran/simplify.c,v
retrieving revision 1.21
diff -c -p -r1.21 simplify.c
*** simplify.c	7 Apr 2005 18:26:37 -0000	1.21
--- simplify.c	10 Apr 2005 05:52:56 -0000
*************** gfc_expr *
*** 2311,2369 ****
  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);
--- 2263,2283 ----
  gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
  {
    gfc_expr *result;
!   int cmp;
  
!   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);
  
!   cmp = mpfr_sgn (s->value.real);
  
!   if (cmp > 0)
!     mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
!   else if (cmp < 0)
!     mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
    else
      {
        gfc_error ("Invalid second argument of NEAREST at %L", &s->where);