I obtain as result for the test program: -0  0;  -0  0

If "a" and "p" are PARAMETERs, the result is:  0  0;  0  0.

This appears to be caused by the use of __builtin_fmod in
trans-intrinsic.c (gfc_conv_intrinsic_mod).  By hacking
the code to disallow the use of the builtin, one uses
the fallback code (which implements the simply code).

program foo
   real a, p, m1, m2
   a = -4.
   p = 2.
   m1 = mod(a, p)
   m2 = a - int(a / p) * p
   print *, m1, m2
end program foo
troutmask:sgk[239] gfc4x -o z t.f90 && ./z
   0.00000000       0.00000000    

Thus, one gets the right answer.  Also, note

#include <stdio.h>
#include <math.h>
int
main(void)
{
  float a, p, m1, m2;
  a = -4.f;
  p = 2.f;
  m1 = fmodf(a, p);
  m2 = a - (int)(a / p) * p;
  printf("%f %f\n", m1, m2);
  return 0;
}
troutmask:sgk[241] cc -o z t.c -lm && ./z
-0.000000 0.000000

which probably means that we either want to not
use the __builtin_fmod() (which is probably some
optimized bit twiddling routine).

There is an additional problem with MOD(A,P) and MODULO(A,P).

In F95, one finds "P = 0, the result is processor dependent."

In F2003 and F2008, one finds "P shall not be zero."

Consider the code,

program foo
   real, parameter :: b = 0.
   real a, p, m1, m2
   a = 2.
   p = 0.
   m1 = mod(a, p)
   m2 = mod(a, b)
   print *, m1, m2
end program foo

with the __builtin_fmod function we get

troutmask:sgk[215] gfc4x -o z t.f90 && ./z
              NaN              NaN

(ie., there is no error or warning that P=0).

If we simply eliminate the __builtin_fmod function,
we get 

troutmask:sgk[212] gfc4x -o z t.f90 && ./z
   2.00000000       2.00000000    

which is worse than the NaN. :(
__ 
Steve

On Mon, May 16, 2011 at 09:31:57PM +0000, sgk at troutmask dot apl.washington.edu wrote:
> In F95, one finds "P = 0, the result is processor dependent."
> 
> In F2003 and F2008, one finds "P shall not be zero."
> 
> Consider the code,
> 
> program foo
>    real, parameter :: b = 0.
>    real a, p, m1, m2
>    a = 2.
>    p = 0.
>    m1 = mod(a, p)
>    m2 = mod(a, b)

for the above case, I forgot to also make 'a' a parameter.
so simplification did not occur.  both mod and modulo is
an error if both 'a' and 'b' are constants and b = 0.
So, we are only missing a runtime error, which should 
probably be triggered with -fno-fast-math.

The fmod behaviour is correct for x < 0 according to N1548:

double fmod(double x, double y);
float fmodf(float x, float y);

The fmod functions return the value x−ny, for some integer n such that, if y is nonzero, the result has the same sign as x and magnitude less than the magnitude of y. If y is zero, whether  a  domain  error  occurs  or  the fmod functions  return  zero  is  implementation-defined.

I suppose we could still use __builtin_fmod if -fno-signed-zeros is in effect.

I suppose we could still use __builtin_fmod if we reset the sign bit if the result is -0.

So does the fallback path actually ever get used? AFAICS the builtins are always available, and if the builtin results in a call to fmod{f,,l,Q} we have fallback implementations in c99_functions.c for non-C99 functions.  See PR 29810.

A problem with a naive implementation of the algorithm specified by the standard is that it doesn't work for large arguments, which was what prompted the usage of builtin_fmod in the first place. See PR 24518.

Taken together, it seems the proper approach would be to just remove the fallback path (not really related to this PR per se, just as a general janitorial patch), and fix the result if it's negative zero. E.g. something like the following pseudocode for MOD(a,p):

if (!options.fast_math && fabs(p) == 0.0)
  generate_error(...) // TBH, aborting the program here seems quite drastic.. [1]
res = __builtin_fmod (a, p)
if (options.signed_zeros)
{
  if (res == -0.0)
    res = 0.0
}

[1] While ISO C leaves it implementation-defined what happens here, POSIX specifies that "If y is zero, a domain error shall occur, and either a NaN (if supported), or an implementation-defined value shall be returned.".  Similarly, MS fmod returns NaN (according to MSDN, I haven't tested). So while not strictly conforming to the Fortran spec, the POSIX/MS approach seems more sensible, and is IMHO a better choice than aborting the program.

On Tue, May 17, 2011 at 06:05:50AM +0000, thenlich at users dot sourceforge.net wrote:
> --- Comment #5 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-05-17 05:51:56 UTC ---
> The fmod behaviour is correct for x < 0 according to N1548:
> 

Yes, I know.

On Tue, May 17, 2011 at 02:17:22PM +0000, jb at gcc dot gnu.org wrote:
>
> So does the fallback path actually ever get used? AFAICS the builtins are
> always available, and if the builtin results in a call to fmod{f,,l,Q} we have
> fallback implementations in c99_functions.c for non-C99 functions.  See PR
> 29810.
> 

I don't know all the targets on which gfortran can run.  Of course,
if all target have the builtin, then the fallback code won't be used.
We can garbage collect the code if you're confident that it is unneeded.

> Taken together, it seems the proper approach would be to just remove the
> fallback path (not really related to this PR per se, just as a general
> janitorial patch), and fix the result if it's negative zero. E.g. something
> like the following pseudocode for MOD(a,p):
> 
> if (!options.fast_math && fabs(p) == 0.0)
>   generate_error(...) // TBH, aborting the program here

The restriction that "P shall not be zero" is on the user
not the compiler, and the compiler is not required to 
diagnosis that problem.  I think we can simply issue a
warning and let the builtin set res = NaN. 

> res = __builtin_fmod (a, p)
> if (options.signed_zeros)
> {
>   if (res == -0.0)
>     res = 0.0
> }

I don't have n1256.pdf handy (draft of C standard), but
IIRC, this check is expensive because zero and negative zero
compare equal.  One needs to explicitly check the sign
bit or simply do

   if (res == 0.0) res = abs(res)

to clear the sign.

(In reply to comment #10)
> On Tue, May 17, 2011 at 02:17:22PM +0000, jb at gcc dot gnu.org wrote:
> >
> > So does the fallback path actually ever get used? AFAICS the builtins are
> > always available, and if the builtin results in a call to fmod{f,,l,Q} we have
> > fallback implementations in c99_functions.c for non-C99 functions.  See PR
> > 29810.
> > 
> 
> I don't know all the targets on which gfortran can run.  Of course,
> if all target have the builtin, then the fallback code won't be used.
> We can garbage collect the code if you're confident that it is unneeded.

I'm fairly confident. AFAICS the idea is that the builtins are always available, but if the compiler cannot optimize it in some meaningful way (e.g. constant folding) then a call to the corresponding library function is emitted, and this is the situation where things can go wrong as not all targets provide a C99 libm. However, as we provide fallbacks in c99_functions.c this case is covered.

This is, in short, the story behind PR 29810; when we started using BUILT_IN_FMOD{F,,L} the float and long double builtins were available also on targets where fmod{f,l} were not.

> > Taken together, it seems the proper approach would be to just remove the
> > fallback path (not really related to this PR per se, just as a general
> > janitorial patch), and fix the result if it's negative zero. E.g. something
> > like the following pseudocode for MOD(a,p):
> > 
> > if (!options.fast_math && fabs(p) == 0.0)
> >   generate_error(...) // TBH, aborting the program here
> 
> The restriction that "P shall not be zero" is on the user
> not the compiler, and the compiler is not required to 
> diagnosis that problem.

Ah, good point.

>  I think we can simply issue a
> warning and let the builtin set res = NaN. 

Well, I was thinking of not testing whether p == 0 at all; As the NaN will presumably propagate the user will find out that something went wrong quickly enough (if the user bothers to check the output, that is! :) ), and finding the place where the NaN occurred is as simple as recompiling with -ffpe-trap=invalid.

> > res = __builtin_fmod (a, p)
> > if (options.signed_zeros)
> > {
> >   if (res == -0.0)
> >     res = 0.0
> > }
> 
> I don't have n1256.pdf handy (draft of C standard), but
> IIRC, this check is expensive because zero and negative zero
> compare equal.  One needs to explicitly check the sign
> bit or simply do
> 
>    if (res == 0.0) res = abs(res)
> 
> to clear the sign.

Ah, so it seems.

Author: jb
Date: Sat May  5 07:59:22 2012
New Revision: 187191

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=187191
Log:
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.


Added:
    trunk/gcc/testsuite/gfortran.dg/mod_large_1.f90
    trunk/gcc/testsuite/gfortran.dg/mod_sign0_1.f90
Modified:
    trunk/gcc/fortran/ChangeLog
    trunk/gcc/fortran/intrinsic.texi
    trunk/gcc/fortran/simplify.c
    trunk/gcc/fortran/trans-intrinsic.c
    trunk/gcc/testsuite/ChangeLog

Closing as fixed. The new behavior should now be identical for constant and non-constant arguments (assuming that at runtime the target inline expands fmod (at least x86(-64) does this) or the library fmod conforms to C99 Annex F).

Wrt. the sign of the result, we now provide the following behavior (again, assuming that also the runtime behavior of fmod conforms to C99 Annex F) which also applies when the result is (signed) zero:

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.

Note that this is not the same as calculating the result according to the formula in the standard using the IEEE 754 rules for signed zero arithmetic, but rather makes sure that the sign behavior is consistent for zero and non-zero results.