Hi,
the attached patch implements a few fixes and cleanups for the MOD and
MODULO intrinsics.
- When the arguments are constant, use mpfr_fmod instead of the naive
algorithms which are numerically unstable for large arguments. This
extends the PR 24518 fix to constant arguments as well, and makes the
compile-time evaluation match the runtime implementation which also
uses fmod in the same manner.
- Remove the old fallback path for the case builtin_fmod is not
available, as the builtin is AFAICS always available.
The patch does not per se fix the corner-case bug as reported in PR
49010, in fact it makes it worse in a way as with the patch the result
if the arguments are parameters is the same as the runtime result
(previously, the compile-time result was correct). But, I think we
should leave it as it is. Due to the reasons above, we're not using
the naive algorithms anyway, and IMHO -0.0 is quite a good
approximation for +0.0 anyway. One might even argue that due to the
numerical instability, specifying the naive algorithms is a bug in the
standard.
The patch adds notes to the documentation about the usage of fmod, so
users interested in corner-case behavior can look up how that function
is supposed to behave on their target. FWIW, AFAICS MPFR and glibc
fmod conform to the behavior specified in C99 Annex F.
Regtested on x86_64-unknown-linux-gnu, Ok for trunk?
2012-03-14 Janne Blomqvist <j...@gcc.gnu.org>
PR fortran/49010
PR fortran/24518
* intrinsic.texi (MOD,MODULO): Mention usage of fmod instead of
naive algorithm.
* simplify.c (gfc_simplify_mod): Use mpfr_fmod.
(gfc_simplify_modulo): Likewise.
* trans-intrinsic.c (gfc_conv_intrinsic_mod): Remove fallback as
builtin_fmod is always available.
--
Janne Blomqvist
diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi
index 294818e..171c5d2 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -8991,8 +8991,7 @@ cases, the result is of the same type and kind as @var{ARRAY}.
@table @asis
@item @emph{Description}:
-@code{MOD(A,P)} computes the remainder of the division of A by P@. It is
-calculated as @code{A - (INT(A/P) * P)}.
+@code{MOD(A,P)} computes the remainder of the division of A by P@.
@item @emph{Standard}:
Fortran 77 and later
@@ -9011,8 +9010,14 @@ equal to zero
@end multitable
@item @emph{Return value}:
-The kind of the return value is the result of cross-promoting
-the kinds of the arguments.
+The return value is the result of @code{A - (INT(A/P) * P)}. The kind
+of the return value is the result of cross-promoting the kinds of the
+arguments.
+
+@item @emph{Note}:
+The obvious algorithm as specified above is unstable for large real
+inputs. Hence, for real inputs the result is calculated by using the
+@code{fmod} function in the C math library.
@item @emph{Example}:
@smallexample
@@ -9082,6 +9087,11 @@ The type and kind of the result are those of the arguments.
@end table
In all cases, if @var{P} is zero the result is processor-dependent.
+@item @emph{Note}:
+The obvious algorithm as specified above is unstable for large real
+inputs. Hence, for real inputs the result is based on using the
+@code{fmod} function in the C math library.
+
@item @emph{Example}:
@smallexample
program test_modulo
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index 706dab4..7e3bc8c 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -4222,7 +4222,6 @@ gfc_expr *
gfc_simplify_mod (gfc_expr *a, gfc_expr *p)
{
gfc_expr *result;
- mpfr_t tmp;
int kind;
if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
@@ -4254,12 +4253,8 @@ gfc_simplify_mod (gfc_expr *a, gfc_expr *p)
}
gfc_set_model_kind (kind);
- mpfr_init (tmp);
- mpfr_div (tmp, a->value.real, p->value.real, GFC_RND_MODE);
- mpfr_trunc (tmp, tmp);
- mpfr_mul (tmp, tmp, p->value.real, GFC_RND_MODE);
- mpfr_sub (result->value.real, a->value.real, tmp, GFC_RND_MODE);
- mpfr_clear (tmp);
+ mpfr_fmod (result->value.real, a->value.real, p->value.real,
+ GFC_RND_MODE);
break;
default:
@@ -4274,7 +4269,6 @@ gfc_expr *
gfc_simplify_modulo (gfc_expr *a, gfc_expr *p)
{
gfc_expr *result;
- mpfr_t tmp;
int kind;
if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
@@ -4308,12 +4302,12 @@ gfc_simplify_modulo (gfc_expr *a, gfc_expr *p)
}
gfc_set_model_kind (kind);
- mpfr_init (tmp);
- mpfr_div (tmp, a->value.real, p->value.real, GFC_RND_MODE);
- mpfr_floor (tmp, tmp);
- mpfr_mul (tmp, tmp, p->value.real, GFC_RND_MODE);
- mpfr_sub (result->value.real, a->value.real, tmp, GFC_RND_MODE);
- mpfr_clear (tmp);
+ mpfr_fmod (result->value.real, a->value.real, p->value.real,
+ GFC_RND_MODE);
+ if ((mpfr_cmp_ui (result->value.real, 0) != 0)
+ && (mpfr_signbit (a->value.real) != mpfr_signbit (p->value.real)))
+ mpfr_add (result->value.real, result->value.real, p->value.real,
+ GFC_RND_MODE);
break;
default:
diff --git a/gcc/fortran/trans-intrinsic.c b/gcc/fortran/trans-intrinsic.c
index ac9f507..0dc3171 100644
--- a/gcc/fortran/trans-intrinsic.c
+++ b/gcc/fortran/trans-intrinsic.c
@@ -1719,21 +1719,24 @@ gfc_conv_intrinsic_cmplx (gfc_se * se, gfc_expr * expr, int both)
se->expr = fold_build2_loc (input_location, COMPLEX_EXPR, type, real, imag);
}
+
/* Remainder function MOD(A, P) = A - INT(A / P) * P
- MODULO(A, P) = A - FLOOR (A / P) * P */
-/* TODO: MOD(x, 0) */
+ MODULO(A, P) = A - FLOOR (A / P) * P
+
+ The obvious algorithms above are numerically instable for large
+ arguments, hence these intrinsics are instead implemented via calls
+ to the fmod family of functions. It is the responsibility of the
+ user to ensure that the second argument is non-zero. */
static void
gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
{
tree type;
- tree itype;
tree tmp;
tree test;
tree test2;
tree fmod;
- mpfr_t huge;
- int n, ikind;
+ tree zero;
tree args[2];
gfc_conv_intrinsic_function_args (se, expr, args, 2);
@@ -1757,16 +1760,15 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
/* Check if we have a builtin fmod. */
fmod = gfc_builtin_decl_for_float_kind (BUILT_IN_FMOD, expr->ts.kind);
- /* Use it if it exists. */
- if (fmod != NULL_TREE)
- {
- tmp = build_addr (fmod, current_function_decl);
- se->expr = build_call_array_loc (input_location,
+ /* The builtin should always be available. */
+ gcc_assert (fmod != NULL_TREE);
+
+ tmp = build_addr (fmod, current_function_decl);
+ se->expr = build_call_array_loc (input_location,
TREE_TYPE (TREE_TYPE (fmod)),
tmp, 2, args);
- if (modulo == 0)
- return;
- }
+ if (modulo == 0)
+ return;
type = TREE_TYPE (args[0]);
@@ -1780,66 +1782,23 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
test = (fmod (arg, arg2) != 0) && ((arg < 0) xor (arg2 < 0))
thereby avoiding another division and retaining the accuracy
of the builtin function. */
- if (fmod != NULL_TREE && modulo)
- {
- tree zero = gfc_build_const (type, integer_zero_node);
- tmp = gfc_evaluate_now (se->expr, &se->pre);
- test = fold_build2_loc (input_location, LT_EXPR, boolean_type_node,
- args[0], zero);
- test2 = fold_build2_loc (input_location, LT_EXPR, boolean_type_node,
- args[1], zero);
- test2 = fold_build2_loc (input_location, TRUTH_XOR_EXPR,
- boolean_type_node, test, test2);
- test = fold_build2_loc (input_location, NE_EXPR, boolean_type_node,
- tmp, zero);
- test = fold_build2_loc (input_location, TRUTH_AND_EXPR,
- boolean_type_node, test, test2);
- test = gfc_evaluate_now (test, &se->pre);
- se->expr = fold_build3_loc (input_location, COND_EXPR, type, test,
- fold_build2_loc (input_location, PLUS_EXPR,
- type, tmp, args[1]), tmp);
- return;
- }
-
- /* If we do not have a built_in fmod, the calculation is going to
- have to be done longhand. */
- tmp = fold_build2_loc (input_location, RDIV_EXPR, type, args[0], args[1]);
-
- /* Test if the value is too large to handle sensibly. */
- gfc_set_model_kind (expr->ts.kind);
- mpfr_init (huge);
- n = gfc_validate_kind (BT_INTEGER, expr->ts.kind, true);
- ikind = expr->ts.kind;
- if (n < 0)
- {
- n = gfc_validate_kind (BT_INTEGER, gfc_max_integer_kind, false);
- ikind = gfc_max_integer_kind;
- }
- mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
- test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind, 0);
+ zero = gfc_build_const (type, integer_zero_node);
+ tmp = gfc_evaluate_now (se->expr, &se->pre);
+ test = fold_build2_loc (input_location, LT_EXPR, boolean_type_node,
+ args[0], zero);
test2 = fold_build2_loc (input_location, LT_EXPR, boolean_type_node,
- tmp, test);
-
- mpfr_neg (huge, huge, GFC_RND_MODE);
- test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind, 0);
- test = fold_build2_loc (input_location, GT_EXPR, boolean_type_node, tmp,
- test);
- test2 = fold_build2_loc (input_location, TRUTH_AND_EXPR,
+ args[1], zero);
+ test2 = fold_build2_loc (input_location, TRUTH_XOR_EXPR,
boolean_type_node, test, test2);
-
- itype = gfc_get_int_type (ikind);
- if (modulo)
- tmp = build_fix_expr (&se->pre, tmp, itype, RND_FLOOR);
- else
- tmp = build_fix_expr (&se->pre, tmp, itype, RND_TRUNC);
- tmp = convert (type, tmp);
- tmp = fold_build3_loc (input_location, COND_EXPR, type, test2, tmp,
- args[0]);
- tmp = fold_build2_loc (input_location, MULT_EXPR, type, tmp, args[1]);
- se->expr = fold_build2_loc (input_location, MINUS_EXPR, type, args[0],
- tmp);
- mpfr_clear (huge);
- break;
+ test = fold_build2_loc (input_location, NE_EXPR, boolean_type_node,
+ tmp, zero);
+ test = fold_build2_loc (input_location, TRUTH_AND_EXPR,
+ boolean_type_node, test, test2);
+ test = gfc_evaluate_now (test, &se->pre);
+ se->expr = fold_build3_loc (input_location, COND_EXPR, type, test,
+ fold_build2_loc (input_location, PLUS_EXPR,
+ type, tmp, args[1]), tmp);
+ return;
default:
gcc_unreachable ();
diff --git a/libgcc/configure b/libgcc/configure
old mode 100644
new mode 100755
