[Patch, fortran] PR24518 - Intrinsic MOD incorrect for large arg1/arg2 and slow.

["Paul Richard Thomas" <paul dot richard dot thomas at gmail dot com>]

:ADDPATCH fortran:

This patch concerned the accuracy of results obtained with  MOD(arg1,
arg2) and MODULO(arg1, arg2) at large values of the ratio of arg1 to
arg2; these were drifting off and exceeding the absolute value of arg2
for relatively small values of the ratio. For example:

 real*8 :: x = 10.0e9
 do i = 10, 22
   x = 10d0 * x
   print '(a,i2,a,g14.8)', "mod (10**",i,", 1.7_8) = ",mod (x, 1.7_8)
 end do
end

gives

mod (10**10, 1.7_8) =  1.3000026
mod (10**11, 1.7_8) =  1.1000261
mod (10**12, 1.7_8) = 0.80026150
mod (10**13, 1.7_8) =  1.2026138
mod (10**14, 1.7_8) = 0.12609863
mod (10**15, 1.7_8) =  1.2607422
mod (10**16, 1.7_8) =  5.8125000
mod (10**17, 1.7_8) = -50.687500
mod (10**18, 1.7_8) =  364.00000
mod (10**19, 1.7_8) = -.70000000E+20
mod (10**20, 1.7_8) = -.70000000E+21
mod (10**21, 1.7_8) = -.70000000E+22
mod (10**22, 1.7_8) = -.70000000E+23

with the patch, which calls built_in fmod,

mod (10**10, 1.7_8) =  1.3000026
mod (10**11, 1.7_8) =  1.1000261
mod (10**12, 1.7_8) = 0.80026123
mod (10**13, 1.7_8) =  1.2026123
mod (10**14, 1.7_8) = 0.12612289
mod (10**15, 1.7_8) =  1.2612289
mod (10**16, 1.7_8) = 0.71228947
mod (10**17, 1.7_8) = 0.32289470
mod (10**18, 1.7_8) =  1.5289470
mod (10**19, 1.7_8) =  1.6894697
mod (10**20, 1.7_8) =  1.5946971
mod (10**21, 1.7_8) = 0.64697063
mod (10**22, 1.7_8) = 0.86970625

is returned, which agrees with a leading brand commercial compiler... or two.

The patch is self explanatory, including the calculation of modulo.
This latter is a bit contorted in order to avoid any further
divisions; this is rewarded by the computation time only being 30%
more than that for mod.

We have not included a testcase because anything useful would fail on
platforms without the fmod built_in.  If none such exist, we would be
happy to provide a testcase!

Regtested on FC5/x86_ia64 - OK for trunk, and 4.2?

We would like to acknowledge Uros Bizjak for his work in ensuring that
the library can provide the goods for any compiler option.  Thanks
Uros!

It so happens that this version of mod and modulo execute pretty quickly too....

FX and Paul

2006-10-31  Francois-Xavier Coudert  <fxcoudert@gcc.gnu,org>
	    Paul Thomas  <pault@gcc.gnu.org>

	PR fortran/24518
	* trans-intrinsic.c (gfc_conv_intrinsic_mod): Use built_in fmod for
	both MOD and MODULO, if it is available.

Index: gcc/fortran/f95-lang.c
===================================================================
*** gcc/fortran/f95-lang.c	(révision 118179)
--- gcc/fortran/f95-lang.c	(copie de travail)
*************** gfc_init_builtin_functions (void)
*** 896,901 ****
--- 896,908 ----
  		      BUILT_IN_COPYSIGN, "copysign", true);
    gfc_define_builtin ("__builtin_copysignf", mfunc_float[1], 
  		      BUILT_IN_COPYSIGNF, "copysignf", true);
+  
+   gfc_define_builtin ("__builtin_fmodl", mfunc_longdouble[1], 
+ 		      BUILT_IN_FMODL, "fmodl", true);
+   gfc_define_builtin ("__builtin_fmod", mfunc_double[1], 
+ 		      BUILT_IN_FMOD, "fmod", true);
+   gfc_define_builtin ("__builtin_fmodf", mfunc_float[1], 
+ 		      BUILT_IN_FMODF, "fmodf", true);
  
    /* These are used to implement the ** operator.  */
    gfc_define_builtin ("__builtin_powl", mfunc_longdouble[1], 
Index: gcc/fortran/trans-intrinsic.c
===================================================================
*** gcc/fortran/trans-intrinsic.c	(révision 118179)
--- gcc/fortran/trans-intrinsic.c	(copie de travail)
*************** gfc_conv_intrinsic_mod (gfc_se * se, gfc
*** 976,989 ****
    int n, ikind;
  
    arg = gfc_conv_intrinsic_function_args (se, expr);
-   arg2 = TREE_VALUE (TREE_CHAIN (arg));
-   arg = TREE_VALUE (arg);
-   type = TREE_TYPE (arg);
  
    switch (expr->ts.type)
      {
      case BT_INTEGER:
        /* Integer case is easy, we've got a builtin op.  */
        if (modulo)
         se->expr = build2 (FLOOR_MOD_EXPR, type, arg, arg2);
        else
--- 976,990 ----
    int n, ikind;
  
    arg = gfc_conv_intrinsic_function_args (se, expr);
  
    switch (expr->ts.type)
      {
      case BT_INTEGER:
        /* Integer case is easy, we've got a builtin op.  */
+       arg2 = TREE_VALUE (TREE_CHAIN (arg));
+       arg = TREE_VALUE (arg);
+       type = TREE_TYPE (arg);
+ 
        if (modulo)
         se->expr = build2 (FLOOR_MOD_EXPR, type, arg, arg2);
        else
*************** gfc_conv_intrinsic_mod (gfc_se * se, gfc
*** 991,1001 ****
        break;
  
      case BT_REAL:
!       /* Real values we have to do the hard way.  */
        arg = gfc_evaluate_now (arg, &se->pre);
        arg2 = gfc_evaluate_now (arg2, &se->pre);
  
        tmp = build2 (RDIV_EXPR, type, arg, arg2);
        /* Test if the value is too large to handle sensibly.  */
        gfc_set_model_kind (expr->ts.kind);
        mpfr_init (huge);
--- 992,1060 ----
        break;
  
      case BT_REAL:
!       n = END_BUILTINS;
!       /* Check if we have a builtin fmod.  */
!       switch (expr->ts.kind)
! 	{
! 	case 4:
! 	  n = BUILT_IN_FMODF;
! 	  break;
! 
! 	case 8:
! 	  n = BUILT_IN_FMOD;
! 	  break;
! 
! 	case 10:
! 	case 16:
! 	  n = BUILT_IN_FMODL;
! 	  break;
! 
! 	default:
! 	  break;
! 	}
! 
!       /* Use it if it exists.  */
!       if (n != END_BUILTINS)
! 	{
! 	  tmp = built_in_decls[n];
! 	  se->expr = build_function_call_expr (tmp, arg);
! 	  if (modulo == 0)
! 	    return;
! 	}
! 
!       arg2 = TREE_VALUE (TREE_CHAIN (arg));
!       arg = TREE_VALUE (arg);
!       type = TREE_TYPE (arg);
! 
        arg = gfc_evaluate_now (arg, &se->pre);
        arg2 = gfc_evaluate_now (arg2, &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 (n != END_BUILTINS && modulo)
+ 	{
+ 	  tree zero = gfc_build_const (type, integer_zero_node);
+ 	  tmp = gfc_evaluate_now (se->expr, &se->pre);
+ 	  test = build2 (LT_EXPR, boolean_type_node, arg, zero);
+ 	  test2 = build2 (LT_EXPR, boolean_type_node, arg2, zero);
+ 	  test2 = build2 (TRUTH_XOR_EXPR, boolean_type_node, test, test2);
+ 	  test = build2 (NE_EXPR, boolean_type_node, tmp, zero);
+ 	  test = build2 (TRUTH_AND_EXPR, boolean_type_node, test, test2);
+ 	  test = gfc_evaluate_now (test, &se->pre);
+ 	  se->expr = build3 (COND_EXPR, type, test,
+ 			     build2 (PLUS_EXPR, type, tmp, arg2), tmp);
+ 	  return;
+ 	}
+ 
+       /* If we do not have a built_in fmod, the calculation is going to
+ 	 have to be done longhand.  */
        tmp = build2 (RDIV_EXPR, type, arg, arg2);
+ 
        /* Test if the value is too large to handle sensibly.  */
        gfc_set_model_kind (expr->ts.kind);
        mpfr_init (huge);

Attachment: Change.Logs
Description: Binary data