[Bug fortran/104535] don't use fmod?

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=104535

kargl at gcc dot gnu.org changed:

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

--- Comment #1 from kargl at gcc dot gnu.org ---
This should be closed with WONTFIX.  gcc/gfortran already provides at
least three options to get the requested behavior; namely, -ffast-math,
-Ofast, and -funsafe-math-optimization.  In addition, gfortran maps mod()
to __builtin_fmod() for REAL arguments.  It is trivial to show that gfortran
compiled Fortran code matches gcc compiled C code.

% cat u.f90
function mymod(x,y)
   real mymod
   real, intent(in), value :: x, y
   mymod = mod(x,y)
end function mymod

% cat u.c
#include <math.h>
float
mymod(float x, float y)
{
    return (fmodf(x,y));
}

% gfortran -O3 -S u.f90
% cat u.s
...
   .cfi_startproc
   jmp   fmodf
   .cfi_endproc
...

% gcc -O3 -S u.f90
% cat u.s
...
   .cfi_startproc
   jmp   fmodf
   .cfi_endproc
...

% gfortran -O3 -S -ffast-math u.f90
% cat u.s
...
   .cfi_startproc
   pushq   %rbp
   .cfi_def_cfa_offset 16
   .cfi_offset 6, -16
   movq    %rsp, %rbp
   .cfi_def_cfa_register 6
   movss   %xmm0, -4(%rbp)
   movss   %xmm1, -8(%rbp)
   flds    -8(%rbp)
   flds    -4(%rbp)
.L2:
   fprem
   fnstsw  %ax
   testb   $4, %ah
   jne     .L2
   fstp    %st(1)
   fstps   -4(%rbp)
   movss   -4(%rbp), %xmm0
   popq    %rbp
   .cfi_def_cfa 7, 8
   .ret
   .cfi_endproc
...


% gcc -O3 -S -ffast-math u.f90
   .cfi_startproc
   pushq   %rbp
   .cfi_def_cfa_offset 16
   .cfi_offset 6, -16
   movq    %rsp, %rbp
   .cfi_def_cfa_register 6
   movss   %xmm0, -4(%rbp)
   movss   %xmm1, -8(%rbp)
   flds    -8(%rbp)
   flds    -4(%rbp)
.L2:
   fprem
   fnstsw  %ax
   testb   $4, %ah
   jne     .L2
   fstp    %st(1)
   fstps   -4(%rbp)
   movss   -4(%rbp), %xmm0
   popq    %rbp
   .cfi_def_cfa 7, 8
   ret
   .cfi_endproc


Now, let's look at trans-intrinsic.cc.  One finds the following comment

/* Remainder function MOD(A, P) = A - INT(A / P) * P
                      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.  */

In infinite precision arithmetic the above equations could use
a straightforward naive in-lined implementation.  Unfortunately,
computers tend to use finite precision arithmetic, so numerical
algorithms must be appropriately considered (See Goldberg).