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. */
+ mpfr_t tmp;
+ int direction, sgn;
- if (x->expr_type != EXPR_CONSTANT)
+ if (x->expr_type != EXPR_CONSTANT || s->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);
+ gfc_set_model_kind (x->ts.kind);
+ result = gfc_copy_expr (x);
- val = mpfr_get_d (x->value.real, GFC_RND_MODE);
- p = gfc_real_kinds[k].digits;
+ direction = mpfr_sgn (s->value.real);
- eps = 1.;
- for (i = 1; i < p; ++i)
+ if (direction == 0)
{
- eps = eps / 2.;
+ gfc_error ("Second argument of NEAREST at %L may not be zero",
+ &s->where);
+ gfc_free (result);
+ return &gfc_bad_expr;
}
- /* 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)
+ /* 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)
{
- if (match_float)
- {
- rval = (float) val;
- rval = rval + eps;
- mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
- }
+ 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
- {
- val = val + eps;
- mpfr_set_d (result->value.real, val, GFC_RND_MODE);
- }
+ 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 (mpfr_cmp_ui (s->value.real, 0) < 0)
+ else
{
- if (match_float)
+ if (sgn < 0)
{
- rval = (float) val;
- rval = rval - eps;
- mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
+ 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
{
- val = val - eps;
- mpfr_set_d (result->value.real, val, GFC_RND_MODE);
+ /* 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);
+
+ /* 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);
+
+ mpfr_clear (tmp);
}
- }
- else
- {
- gfc_error ("Invalid second argument of NEAREST at %L", &s->where);
- gfc_free (result);
- return &gfc_bad_expr;
+
+ if (sgn < 0)
+ mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
}
return range_check (result, "NEAREST");
--- /dev/null
+! { 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 ()
+! These last two tests are commented out because mpfr provides no support
+! for denormals, and therefore we get TINY instead of the correct result.
+!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