Re: [Patch, fortran] PR 49010/24518 MOD/MODULO fixes

Re: [Patch, fortran] PR 49010/24518 MOD/MODULO fixes

[Tobias Burnus]

Janne Blomqvist wrote:
> the attached patch implements a few fixes and cleanups for the MOD and
> MODULO intrinsics.

> 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.

> +@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.

I wonder whether one should extend the note, stating that that using
"fmod" might lead to a `wrongly' signed 0. Something like: "; depending
on its implementation, this might lead to a diffent sign for 0 results
- compared to the result of the obvious algorithm."

Since you modify intrinsic.texi, could you also do:
- Add a cross @ref between mod and modulo.
- Show in the example (as comments) also the result of the mod/module
  operations. As many confuse mod (=remainder) and modulo, it makes
  sense to help them by showing the result of examples.*

Please also update the Copyright year for simplify.c.

Otherwise, the patch looks OK.

(Except for the mode change of libgcc/configure - which requires a build
maintainer approval. If you want to make it executable, do the same also
for libitm/configure - and do so in a separate patch.)

Tobias


* Regarding "mod" vs. "module" see also
https://en.wikipedia.org/wiki/Remainder#The_case_of_general_integers
and the right column at
https://en.wikipedia.org/wiki/Modulo_operation

Re: [Patch, fortran] PR 49010/24518 MOD/MODULO fixes

[Janne Blomqvist]

On Wed, Apr 4, 2012 at 12:11, Tobias Burnus
<tobias.burnus@physik.fu-berlin.de> wrote:
> Janne Blomqvist wrote:
>> the attached patch implements a few fixes and cleanups for the MOD and
>> MODULO intrinsics.
>
>> 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.
>
>> +@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.
>
> I wonder whether one should extend the note, stating that that using
> "fmod" might lead to a `wrongly' signed 0. Something like: "; depending
> on its implementation, this might lead to a diffent sign for 0 results
> - compared to the result of the obvious algorithm."

Yes, but if one describes the behavior for one special case, IMHO one
should also specify it for other special cases. E.g. from the glibc
fmod manpage:

...the returned value has the same sign as x and a magnitude less than
the magnitude of y.

       If x or y is a NaN, a NaN is returned.

       If x is an infinity, a domain error occurs, and a NaN is returned.

       If y is zero, a domain error occurs, and a NaN is returned.

       If x is +0 (-0), and y is not zero, +0 (-0) is returned.

For compile-time evaluation, MPFR fmod conforms to the specification
for fmod in C99 Annex F, and for runtime so does glibc fmod (see
above). But for other targets, runtime might be different. But I
couldn't come up with a formulation that I was happy with; not that
I'm happy about my text in the patch either. Suggestions welcome..

> Since you modify intrinsic.texi, could you also do:
> - Add a cross @ref between mod and modulo.

Ok.

> - Show in the example (as comments) also the result of the mod/module
>  operations. As many confuse mod (=remainder) and modulo, it makes
>  sense to help them by showing the result of examples.*

Well, Fortran mod != C99 remainder().

So for a general remainder operation the result is "x-n*y" (evaluated
with infinite precision), where n is x/y rounded to an integer in some
way. Then we have at least the following cases

- C fmod and Fortran MOD: n is rounded towards zero

- Fortran MODULO: n is rounded towards -Infinity

- C99 and IEEE754 remainder: n is rounded to the nearest integer.
Incidentally, this has the nice property that abs(result) <= abs(y/2).

> Please also update the Copyright year for simplify.c.

Ok.

> Otherwise, the patch looks OK.
>
> (Except for the mode change of libgcc/configure - which requires a build
> maintainer approval. If you want to make it executable, do the same also
> for libitm/configure - and do so in a separate patch.)

Ah, that was some autogenerated stuff that was accidentally included
in the patch. Please disregard it.


-- 
Janne Blomqvist

Re: [Patch, fortran] PR 49010/24518 MOD/MODULO fixes

[Janne Blomqvist]

PING**2

On Wed, Mar 21, 2012 at 23:45, Janne Blomqvist
<blomqvist.janne@gmail.com> wrote:
> PING
>
> On Wed, Mar 14, 2012 at 01:03, Janne Blomqvist
> <blomqvist.janne@gmail.com> wrote:
>> 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  <jb@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
>
>
>
> --
> Janne Blomqvist



-- 
Janne Blomqvist

Re: [Patch, fortran] PR 49010/24518 MOD/MODULO fixes

[Janne Blomqvist]

PING

On Wed, Mar 14, 2012 at 01:03, Janne Blomqvist
<blomqvist.janne@gmail.com> wrote:
> 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  <jb@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



-- 
Janne Blomqvist

[Patch, fortran] PR 49010/24518 MOD/MODULO fixes

[Janne Blomqvist]

[-- Attachment #1: Type: text/plain, Size: 1761 bytes --]

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  <jb@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

[-- Attachment #2: mod.diff --]
[-- Type: text/x-patch, Size: 8749 bytes --]

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