PR 49010,24518 MOD/MODULO fixes.

gcc/fortran:

2012-05-05  Janne Blomqvist  <jb@gcc.gnu.org>

	PR fortran/49010
	PR fortran/24518
	* intrinsic.texi (MOD, MODULO): Mention sign and magnitude of result.
	* simplify.c (gfc_simplify_mod): Use mpfr_fmod.
	(gfc_simplify_modulo): Likewise, use copysign to fix the result if
	zero.
	* trans-intrinsic.c (gfc_conv_intrinsic_mod): Remove fallback as
	builtin_fmod is always available. For modulo, call copysign to fix
	the result when signed zeros are enabled.


testsuite:

2012-05-05  Janne Blomqvist  <jb@gcc.gnu.org>

	PR fortran/49010
	PR fortran/24518
	* gfortran.dg/mod_sign0_1.f90: New test.
	* gfortran.dg/mod_large_1.f90: New test.

From-SVN: r187191

gcc/fortran/ChangeLog

 2012-05-05  Janne Blomqvist  <jb@gcc.gnu.org>
 
+	PR fortran/49010
+	PR fortran/24518
+	* intrinsic.texi (MOD, MODULO): Mention sign and magnitude of result.
+	* simplify.c (gfc_simplify_mod): Use mpfr_fmod.
+	(gfc_simplify_modulo): Likewise, use copysign to fix the result if
+	zero.
+	* trans-intrinsic.c (gfc_conv_intrinsic_mod): Remove fallback as
+	builtin_fmod is always available. For modulo, call copysign to fix
+	the result when signed zeros are enabled.
+
+2012-05-05  Janne Blomqvist  <jb@gcc.gnu.org>
+
         * gfortran.texi (GFORTRAN_TMPDIR): Rename to TMPDIR, explain
         algorithm for choosing temp directory.
 

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
@@ -9005,14 +9004,16 @@ Elemental function
 
 @item @emph{Arguments}:
 @multitable @columnfractions .15 .70
-@item @var{A} @tab Shall be a scalar of type @code{INTEGER} or @code{REAL}
-@item @var{P} @tab Shall be a scalar of the same type as @var{A} and not
-equal to zero
+@item @var{A} @tab Shall be a scalar of type @code{INTEGER} or @code{REAL}.
+@item @var{P} @tab Shall be a scalar of the same type and kind as @var{A} 
+and not 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 type
+and kind of the return value is the same as that of the arguments. The
+returned value has the same sign as A and a magnitude less than the
+magnitude of P.
 
 @item @emph{Example}:
 @smallexample
@@ -9041,6 +9042,10 @@ end program test_mod
 @item @code{AMOD(A,P)} @tab @code{REAL(4) A,P} @tab @code{REAL(4)} @tab Fortran 95 and later
 @item @code{DMOD(A,P)} @tab @code{REAL(8) A,P} @tab @code{REAL(8)} @tab Fortran 95 and later
 @end multitable
+
+@item @emph{See also}:
+@ref{MODULO}
+
 @end table
 
 
@@ -9066,8 +9071,9 @@ Elemental function
 
 @item @emph{Arguments}:
 @multitable @columnfractions .15 .70
-@item @var{A} @tab Shall be a scalar of type @code{INTEGER} or @code{REAL}
-@item @var{P} @tab Shall be a scalar of the same type and kind as @var{A}
+@item @var{A} @tab Shall be a scalar of type @code{INTEGER} or @code{REAL}.
+@item @var{P} @tab Shall be a scalar of the same type and kind as @var{A}. 
+It shall not be zero.
 @end multitable
 
 @item @emph{Return value}:
@@ -9080,7 +9086,8 @@ The type and kind of the result are those of the arguments.
 @item If @var{A} and @var{P} are of type @code{REAL}:
 @code{MODULO(A,P)} has the value of @code{A - FLOOR (A / P) * P}.
 @end table
-In all cases, if @var{P} is zero the result is processor-dependent.
+The returned value has the same sign as P and a magnitude less than
+the magnitude of P.
 
 @item @emph{Example}:
 @smallexample
@@ -9096,6 +9103,9 @@ program test_modulo
 end program
 @end smallexample
 
+@item @emph{See also}:
+@ref{MOD}
+
 @end table
 
 

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,17 @@ 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)
+	  {
+	    if (mpfr_signbit (a->value.real) != mpfr_signbit (p->value.real))
+	      mpfr_add (result->value.real, result->value.real, p->value.real,
+			GFC_RND_MODE);
+	  }
+	else
+	  mpfr_copysign (result->value.real, result->value.real, 
+			 p->value.real, GFC_RND_MODE);
 	break;
 
       default:

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]);
 
@@ -1774,16 +1776,31 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
       args[1] = gfc_evaluate_now (args[1], &se->pre);
 
       /* Definition:
-	 modulo = arg - floor (arg/arg2) * arg2, so
-		= test ? fmod (arg, arg2) : fmod (arg, arg2) + arg2, 
-	 where
-	  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)
+	 modulo = arg - floor (arg/arg2) * arg2
+
+	 In order to calculate the result accurately, we use the fmod
+	 function as follows.
+	 
+	 res = fmod (arg, arg2);
+	 if (res)
+	   {
+	     if ((arg < 0) xor (arg2 < 0))
+	       res += arg2;
+	   }
+	 else
+	   res = copysign (0., arg2);
+
+	 => As two nested ternary exprs:
+
+	 res = res ? (((arg < 0) xor (arg2 < 0)) ? res + arg2 : res) 
+	       : copysign (0., arg2);
+
+      */
+
+      zero = gfc_build_const (type, integer_zero_node);
+      tmp = gfc_evaluate_now (se->expr, &se->pre);
+      if (!flag_signed_zeros)
 	{
-	  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,
@@ -1796,50 +1813,35 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
 				  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;
+				      fold_build2_loc (input_location, 
+						       PLUS_EXPR,
+						       type, tmp, args[1]), 
+				      tmp);
 	}
-
-      /* 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)
+      else
 	{
-	  n = gfc_validate_kind (BT_INTEGER, gfc_max_integer_kind, false);
-	  ikind = gfc_max_integer_kind;
+	  tree expr1, copysign, cscall;
+	  copysign = gfc_builtin_decl_for_float_kind (BUILT_IN_COPYSIGN, 
+						      expr->ts.kind);
+	  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);
+	  expr1 = fold_build3_loc (input_location, COND_EXPR, type, test2,
+				   fold_build2_loc (input_location, 
+						    PLUS_EXPR,
+						    type, tmp, args[1]), 
+				   tmp);
+	  test = fold_build2_loc (input_location, NE_EXPR, boolean_type_node,
+				  tmp, zero);
+	  cscall = build_call_expr_loc (input_location, copysign, 2, zero, 
+					args[1]);
+	  se->expr = fold_build3_loc (input_location, COND_EXPR, type, test,
+				      expr1, cscall);
 	}
-      mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
-      test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind, 0);
-      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,
-			       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;
+      return;
 
     default:
       gcc_unreachable ();

gcc/testsuite/ChangeLog

+2012-05-05  Janne Blomqvist  <jb@gcc.gnu.org>
+
+	PR fortran/49010
+	PR fortran/24518
+	* gfortran.dg/mod_sign0_1.f90: New test.
+	* gfortran.dg/mod_large_1.f90: New test.
+
 2012-05-04  Tobias Burnus  <burnus@net-b.de>
 
 	PR fortran/53175

gcc/testsuite/gfortran.dg/mod_large_1.f90

+! { dg-do run }
+! PR fortran/24518 
+! MOD/MODULO of large arguments.
+! The naive algorithm goes pear-shaped for large arguments, instead
+! use fmod.
+! Here we test only with constant arguments (evaluated with
+! mpfr_fmod), as we don't want to cause failures on targets with a
+! crappy libm.
+program mod_large_1
+  implicit none
+  real :: r1
+  r1 = mod (1e22, 1.7)
+  if (abs(r1 - 0.995928764) > 1e-5) call abort
+  r1 = modulo (1e22, -1.7)
+  if (abs(r1 + 0.704071283) > 1e-5) call abort
+end program mod_large_1

gcc/testsuite/gfortran.dg/mod_sign0_1.f90

+! { dg-do run }
+! PR fortran/49010 
+! MOD/MODULO sign of zero.
+
+! We wish to provide the following guarantees:
+
+! MOD(A, P): The result has the sign of A and a magnitude less than
+! that of P.  
+
+! MODULO(A, P): The result has the sign of P and a magnitude less than
+! that of P.
+
+! Here we test only with constant arguments (evaluated with
+! mpfr_fmod), as we don't want to cause failures on targets with a
+! crappy libm. But, a target where fmod follows C99 Annex F is
+! fine. Also, targets where GCC inline expands fmod (such as x86(-64))
+! are also fine.
+program mod_sign0_1
+  implicit none
+  real :: r, t
+
+  r = mod (4., 2.)
+  t = sign (1., r)
+  if (t < 0.) call abort
+
+  r = modulo (4., 2.)
+  t = sign (1., r)
+  if (t < 0.) call abort
+
+  r = mod (-4., 2.)
+  t = sign (1., r)
+  if (t > 0.) call abort
+
+  r = modulo (-4., 2.)
+  t = sign (1., r)
+  if (t < 0.) call abort
+
+  r = mod (4., -2.)
+  t = sign (1., r)
+  if (t < 0.) call abort
+
+  r = modulo (4., -2.)
+  t = sign (1., r)
+  if (t > 0.) call abort
+
+  r = mod (-4., -2.)
+  t = sign (1., r)
+  if (t > 0.) call abort
+
+  r = modulo (-4., -2.)
+  t = sign (1., r)
+  if (t > 0.) call abort
+
+end program mod_sign0_1