[Bug fortran/49010] New: Result of MOD and MODULO intrinsic has wrong sign

[Bug fortran/49010] New: Result of MOD and MODULO intrinsic has wrong sign

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

           Summary: Result of MOD and MODULO intrinsic has wrong sign
           Product: gcc
           Version: 4.7.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: thenlich@users.sourceforge.net


If the first argument to the MOD and MODULO intrinsic is negative, and the
magnitude of the result is zero, the function returns a negative zero.

Fortran 2003/2008:
MOD: The value of the result is A − INT (A/P) * P.
MODULO: The value of the result is A − FLOOR (A / P) * P.

The values returned by these intrinsics should be the same as the expressions
by which the functions are defined; which is zero, without a negative sign.

Example:
program test_mod
    real :: a, p

    a = -4.0
    p = 2.0
    print *, mod(a, p), a - int(a / p) * p ! -0 0 expected 0 0
    print *, modulo(a, p), a - floor(a / p) * p ! -0 0 expected 0 0
end program test_mod

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

Tobias Burnus <burnus at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu.org

--- Comment #1 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-05-16 10:18:20 UTC ---
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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[kargl at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org

--- Comment #2 from kargl at gcc dot gnu.org 2011-05-16 20:44:08 UTC ---
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).

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #3 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2011-05-16 21:17:44 UTC ---
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

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #4 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2011-05-16 21:43:57 UTC ---
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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

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

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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

Janne Blomqvist <jb at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jb at gcc dot gnu.org

--- Comment #6 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-05-17 09:32:18 UTC ---
I suppose we could still use __builtin_fmod if -fno-signed-zeros is in effect.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #7 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-05-17 11:57:31 UTC ---
I suppose we could still use __builtin_fmod if we reset the sign bit if the
result is -0.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #8 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-05-17 14:02:07 UTC ---
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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #9 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2011-05-17 14:02:11 UTC ---
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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #10 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2011-05-17 14:50:52 UTC ---
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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #11 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-05-17 16:18:41 UTC ---
(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.

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

--- Comment #12 from Janne Blomqvist <jb at gcc dot gnu.org> 2012-05-05 07:59:28 UTC ---
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

[Bug fortran/49010] Result of MOD and MODULO intrinsic has wrong sign

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49010

Janne Blomqvist <jb at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |FIXED

--- Comment #13 from Janne Blomqvist <jb at gcc dot gnu.org> 2012-05-05 08:11:35 UTC ---
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.