wrong results for SPACING and RRSPACING at compile time
Steve Kargl
sgk@troutmask.apl.washington.edu
Tue Oct 3 21:03:00 GMT 2006
On Sun, Oct 01, 2006 at 04:15:35PM +0200, Dominique Dhumieres wrote:
>
> real x, xd, xu
> x = -3.0*epsilon(x)
> xd = abs(fraction(x))*set_exponent(real(radix(x), kind=kind(x)), 1+digits(x))
> print *, -3.0*epsilon(1.0)-x
> print *, x, spacing(x), spacing(-3.0*epsilon(1.0))
> print *, x, xd, rrspacing(x), rrspacing(-3.0*epsilon(1.0))
> end
>
> gives
>
> 0.000000
> -3.5762787E-07 2.8421709E-14 5.6843419E-14
> -3.5762787E-07 1.2582912E+07 1.2582912E+07 6291456.
>
> where the last column of the two last lines computed from
> spacing(-3.0*epsilon(1.0)) and rrspacing(-3.0*epsilon(1.0))
> (hence the 'compile time' in the subject) do not match the
> values computed through a variable.
>
Dominique,
I have some fixes for your spacing/rrspacing issues. There is
one failure mode that presists
troutmask:sgk[210] gfc4x -o z z.f90
troutmask:sgk[211] ./z
0.000000
-4.2038954E-45 1.1754944E-38 1.1754944E-38
-4.2038954E-45 NaN 8.5070592E+37
troutmask:sgk[212] cat z.f90
x = -3.0*(nearest(tiny(x),1.0)-tiny(x))
print *, -3.0*(nearest(tiny(1.0),1.0)-tiny(1.0))-x
print *, x, spacing(x), spacing(-3.0*(nearest(tiny(1.0),1.0)-tiny(1.0)))
print *, x, rrspacing(x), rrspacing(-3.0*(nearest(tiny(1.0),1.0)-tiny(1.0)))
end
The NaN is coming from trans-intrinsics.c(gfc_conv_intrinsic_rrspacing),
which claims to handle subnormal numbers. This function is inlining the
rrspacing function and at the moment it is magick to me. From what
I've been able to find on the web:
http://j3-fortran.org/doc/year/05/05-153.txt
The result of rrspacing is:
if (x = 0)
rrspacing(x) = 0;
else if (0 < x < tiny(x))
rrspacing(x) = 1 / tiny(x);
else
rrspacing(x) = |x * b**(-e)| * b**p = |x|*b**(p-e)
where b, e, and p are the radix, exponent of x, and precision.
I've attached the diff, if you want to give it a spin.
--
Steve
-------------- next part --------------
Index: gcc/fortran/arith.c
===================================================================
--- gcc/fortran/arith.c (revision 117407)
+++ gcc/fortran/arith.c (working copy)
@@ -438,14 +438,16 @@ gfc_check_real_range (mpfr_t p, int kind
gfc_free (bin);
#else
mp_exp_t emin, emax;
+ int en;
/* 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[i].min_exponent - 1);
- mpfr_set_emax ((mp_exp_t) gfc_real_kinds[i].max_exponent - 1);
+ en = gfc_real_kinds[i].min_exponent - gfc_real_kinds[i].digits + 1;
+ mpfr_set_emin ((mp_exp_t) en);
+ mpfr_set_emax ((mp_exp_t) gfc_real_kinds[i].max_exponent);
mpfr_subnormalize (q, 0, GFC_RND_MODE);
/* Reset emin and emax. */
Index: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c (revision 117407)
+++ gcc/fortran/simplify.c (working copy)
@@ -2497,8 +2497,6 @@ gfc_simplify_nearest (gfc_expr * x, gfc_
int sgn;
#if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
int direction;
-#else
- mp_exp_t emin, emax;
#endif
if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
@@ -2561,25 +2559,10 @@ gfc_simplify_nearest (gfc_expr * x, gfc_
mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
}
#else
-
- /* Save current values of emin and emax. */
- emin = mpfr_get_emin ();
- emax = mpfr_get_emax ();
-
- /* Set emin and emax for the current model number. */
- sgn = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
- 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_set_emin (emin);
- mpfr_set_emax (emax);
-
mpfr_clear(tmp);
#endif
@@ -3111,6 +3094,7 @@ bad_reshape:
}
+#if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
gfc_expr *
gfc_simplify_rrspacing (gfc_expr * x)
{
@@ -3131,7 +3115,7 @@ gfc_simplify_rrspacing (gfc_expr * x)
if (mpfr_sgn (x->value.real) == 0)
{
- mpfr_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
return result;
}
@@ -3160,7 +3144,47 @@ gfc_simplify_rrspacing (gfc_expr * x)
return range_check (result, "RRSPACING");
}
+#else
+gfc_expr *
+gfc_simplify_rrspacing (gfc_expr * x)
+{
+ gfc_expr *result;
+ int i;
+ long int e, p;
+
+ if (x->expr_type != EXPR_CONSTANT)
+ return NULL;
+
+ i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
+ result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
+
+ mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
+
+ /* Special case x = 0 and 0. */
+ if (mpfr_sgn (result->value.real) == 0)
+ {
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
+ return result;
+ }
+
+ /* Special case subnormal numbers. */
+ if (mpfr_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
+ {
+ mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
+ mpfr_div (result->value.real, result->value.real,
+ gfc_real_kinds[i].tiny, GFC_RND_MODE);
+ return result;
+ }
+
+ /* | x * 2**(-e) | * 2**p = | x | * 2**(p-e). */
+ e = (long int) mpfr_get_exp (x->value.real);
+ p = (long int) gfc_real_kinds[i].digits;
+ mpfr_mul_2si (result->value.real, result->value.real, p - e, GFC_RND_MODE);
+
+ return range_check (result, "RRSPACING");
+}
+#endif
gfc_expr *
gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
@@ -3604,7 +3628,7 @@ gfc_simplify_sngl (gfc_expr * a)
return range_check (result, "SNGL");
}
-
+#if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
gfc_expr *
gfc_simplify_spacing (gfc_expr * x)
{
@@ -3624,16 +3648,16 @@ gfc_simplify_spacing (gfc_expr * x)
gfc_set_model_kind (x->ts.kind);
- if (mpfr_sgn (x->value.real) == 0)
+ /* Special case x = 0 and -0. */
+ mpfr_init (absv);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ if (mpfr_sgn (absv) == 0)
{
mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
return result;
}
mpfr_init (log2);
- mpfr_init (absv);
-
- mpfr_abs (absv, x->value.real, GFC_RND_MODE);
mpfr_log2 (log2, absv, GFC_RND_MODE);
mpfr_trunc (log2, log2);
@@ -3655,7 +3679,44 @@ gfc_simplify_spacing (gfc_expr * x)
return range_check (result, "SPACING");
}
+#else
+gfc_expr *
+gfc_simplify_spacing (gfc_expr * x)
+{
+ gfc_expr *result;
+ int i;
+ long int en, ep;
+
+ if (x->expr_type != EXPR_CONSTANT)
+ return NULL;
+
+ i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
+
+ result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
+
+ /* Special case x = 0 and -0. */
+ mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
+ if (mpfr_sgn (result->value.real) == 0)
+ {
+ mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+ return result;
+ }
+
+ /* In the Fortran 95 standard, the result is b**(e - p) where b, e, and p
+ are the radix, exponent of x, and precision. This excludes the
+ possibility of subnormal numbers. Fortran 2003 states the result is
+ b**max(e - p, emin - 1). */
+
+ ep = (long int) mpfr_get_exp (x->value.real) - gfc_real_kinds[i].digits;
+ en = (long int) gfc_real_kinds[i].min_exponent - 1;
+ en = en > ep ? en : ep;
+
+ mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
+ mpfr_mul_2si (result->value.real, result->value.real, en, GFC_RND_MODE);
+ return range_check (result, "SPACING");
+}
+#endif
gfc_expr *
gfc_simplify_sqrt (gfc_expr * e)
More information about the Fortran
mailing list