Problem with nearest(tiny(o),-1.0)/2.0

[Steve Kargl]

On Tue, Apr 04, 2006 at 11:57:12AM +0200, Dominique Dhumieres wrote:
> The following program
> 
> program chop
>   real o, t, td, tu, x, y
>   o = 1.
>   t = tiny(o)
>   td = nearest(t,-1.0)
>   x = td/2.0
>   y = nearest(tiny(o),-1.0)/2.0
>   print *, x, y, x - y
> end program chop
> 
> gives
> 
>   5.8774718E-39  5.8774704E-39  1.4012985E-45
> 
> on OSX 10.3.9 and GNU Fortran 95 (GCC) 4.2.0 20060401 (experimental).
> I.e, nearest(tiny(o),-1.0)/2.0 does not give the same result as
> 
>   td = nearest(t,-1.0)
>   x = td/2.0
> 

With the attached patch, I get

troutmask:sgk[273] ./z
  5.8774718E-39  5.8774718E-39   0.000000    

This patch requires MPFR 2.2.0 or higher.  I'm still testing 
this patch, so it may change.  I also need to update configure
to look for a correct version of MPFR.

-- 
Steve
-------------- next part --------------
Index: simplify.c
===================================================================
--- simplify.c	(revision 112712)
+++ simplify.c	(working copy)
@@ -2385,71 +2385,42 @@ gfc_simplify_nearest (gfc_expr * x, gfc_
 {
   gfc_expr *result;
   mpfr_t tmp;
-  int direction, sgn;
+  int sgn;
+  mp_exp_t emin, emax;
 
   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)
+  if (mpfr_sgn (s->value.real) == 0)
     {
-      gfc_error ("Second argument of NEAREST at %L may not be zero",
-		 &s->where);
-      gfc_free (result);
+      gfc_error ("Second argument of NEAREST at %L shall not be zero", &s->where);
       return &gfc_bad_expr;
     }
 
-  /* TODO: Use mpfr_nextabove and mpfr_nextbelow once we move to a
-     newer version of mpfr.  */
+  gfc_set_model_kind (x->ts.kind);
 
-  sgn = mpfr_sgn (x->value.real);
+  sgn = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
 
-  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].subnormal, GFC_RND_MODE);
-      else
-	mpfr_sub (result->value.real,
-		  x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
-    }
-  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 subtract 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);
+  result = gfc_copy_expr (x);
+       
+  /* Save current values of emin and emax.  */
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
+  /* Set emin and emax for the current model number.  */
+  mpfr_set_emin ((mp_exp_t) gfc_real_kinds[sgn].min_exponent - 1);
+  mpfr_set_emax ((mp_exp_t) gfc_real_kinds[sgn].max_exponent - 1);
+
+  sgn = mpfr_sgn (s->value.real); 
+  mpfr_init (tmp);
+  mpfr_set_inf (tmp, sgn);
+  mpfr_nexttoward (result->value.real, tmp);
+  mpfr_subnormalize (result->value.real, 0, GFC_RND_MODE);
 
-	  mpfr_clear (tmp);
-	}
+  mpfr_set_emin (emin);
+  mpfr_set_emax (emax);
 
-      if (sgn < 0)
-	mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
-    }
+  mpfr_clear(tmp);
 
   return range_check (result, "NEAREST");
 }