[Bug libfortran/99740] New: floating point exception in rand() in gfortran

[Bug libfortran/99740] New: floating point exception in rand() in gfortran

[pvoytas at gmail dot com]

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

            Bug ID: 99740
           Summary: floating point exception in rand() in gfortran
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libfortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: pvoytas at gmail dot com
  Target Milestone: ---

Possible bug in gfortran rand() function?

Included code fails with IEEE divide by zero signalling. 
When I request a backtrace, it is from a line with "-log(rand(0))". 

Last output to the file for code included below is 0.183492497 on line 7372837.
If I just write the rand(0) to file (not taking the log), the code finishes
with no errors and the values near line 7372837 don't seem to have any weird
values that would account for the error.

Code runs with out errors either way if I use "call random_number(x)" to
generate random numbers. Since I'm porting old code, the backward compatibility
of rand(0) is needed.

(note: Code compiles and runs successfully with either g77 or gfortran on
Centos 6.x machine with distro default gcc)

gcc -v output:
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/opt/ohpc/pub/compiler/gcc/8.3.0/libexec/gcc/x86_64-pc-linux-gnu/8.3.0/lto-wrapper
Target: x86_64-pc-linux-gnu
Configured with: ../configure --disable-multilib
--enable-languages=c,c++,fortran --prefix=/opt/ohpc/pub/compiler/gcc/8.3.0
--disable-static --enable-shared
Thread model: posix
gcc version 8.3.0 (GCC)

Distro:
Oracle Linux (7.9): uname -a gives:
Linux <hostname> 5.4.17-2011.6.2.el7uek.x86_64 #2 SMP Thu Sep 3 14:09:14 PDT
2020 x86_64 x86_64 x86_64 GNU/Linux

compile command:
gfortran -ffpe-trap=zero -g rantest.f

compiles with no errors

runtime errors:
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic
operation.

Backtrace for this error:
#0  0x7fc09a5113ff in ???
#1  0x7fc09ad330bc in ???
#2  0x4008cf in rantest
        at /home/pvoytas/ps_cp/accidentals/time_series/rantest.f:10
#3  0x400968 in main
        at /home/pvoytas/ps_cp/accidentals/time_series/rantest.f:13
Floating exception (core dumped)

code:
      program rantest
      implicit none
      integer*4 i
      real*8 x

      do i=1,10000000
c       call random_number(x)
c       write(11,*)x
c       write(11,*)rand(0)
       write(11,*)-log(rand(0))
      enddo
      stop
      end

[Bug libfortran/99740] floating point exception in rand() in gfortran

[anlauf at gcc dot gnu.org]

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

anlauf at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Last reconfirmed|                            |2021-03-23
     Ever confirmed|0                           |1
             Status|UNCONFIRMED                 |WAITING

--- Comment #1 from anlauf at gcc dot gnu.org ---
rand(3f) is not a Fortran standard function but provided by your Linux system.

The man page says:

       rand and drand return values in the range 0.0 through 1.0.

So rand may return 0. and an FP exception is not surprising.

If you agree we close your report as invalid.

[Bug libfortran/99740] floating point exception in rand() in gfortran

[kargl at gcc dot gnu.org]

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

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|WAITING                     |RESOLVED
                 CC|                            |kargl at gcc dot gnu.org
           Priority|P3                          |P4
         Resolution|---                         |INVALID

--- Comment #2 from kargl at gcc dot gnu.org ---
It is not rand().   Modifying your code gives

      program rantest
      implicit none
      integer*4 i
      real*4 x
      do i = 1, 7372933
        write(11,*) -log(rand(0))
      end do
      x = rand(0)
      write(*,*) x
      end

% gfcx -o z -O -ffpe-trap=zero,invalid -g a.f && ./z
   0.00000000    

log() is from your math library, which likely documents that log(0) returns
divide-by-zero exception.  Nothing here to fix.

[Bug libfortran/99740] floating point exception in rand() in gfortran

[pvoytas at gmail dot com]

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

--- Comment #3 from Paul A. Voytas <pvoytas at gmail dot com> ---
I see what you mean--if i test for rand(0)=0.d0 I do get hist with gfortran on
the EL7 machine. 

But it seems like there must still be something different from older versions.
The info pages for rand() say "between 0 and 1" which I always took to be
exclusive of the endpoints (of course there's machine precision). On CentOS6
with g77 when I run the above code I get no errors--even with 10x more
attempts--and the test for rand(0)=0.d0 never is true. On that same CentOS6
machine with gfortran, code does show rand(0)=0.d0 cases (and the -log(rand(0))
returns +Infinity.

[Bug libfortran/99740] floating point exception in rand() in gfortran

[sgk at troutmask dot apl.washington.edu]

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

--- Comment #4 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Thu, Mar 25, 2021 at 12:52:53PM +0000, pvoytas at gmail dot com wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=99740
> 
> --- Comment #3 from Paul A. Voytas <pvoytas at gmail dot com> ---
> I see what you mean--if i test for rand(0)=0.d0 I do get hist with gfortran on
> the EL7 machine. 
> 
> But it seems like there must still be something different from older versions.
> The info pages for rand() say "between 0 and 1" which I always took to be
> exclusive of the endpoints (of course there's machine precision). On CentOS6
> with g77 when I run the above code I get no errors--even with 10x more
> attempts--and the test for rand(0)=0.d0 never is true. On that same CentOS6
> machine with gfortran, code does show rand(0)=0.d0 cases (and the -log(rand(0))
> returns +Infinity.
> 

g77 has not been supported for nearly 1.5 decades (gcc-3.4.6
released on 6 Mar 2006).  gfortran tries to provide some level
of API compatibility with g77.  I believe there has never been
a commitment to provide bit-for-bit compatibility (especially
for something like rand()).

Having now looked at g77's rand() implementation, I see that it
is a wrapper around the libc routine by the same name rand(3).
The quality of implementation of libc's rand(3) certainly 
varies across operating systems.  gfortran's rand() implements
the classic Park and Miller LCG RNG with an internal range of
[0, 2**31-1), which yields a random number in [0., 1.).

You really want to switch to using random_number().  The internal
generator used in it has an enormous period and very good
statistical quality, and random_number() is specified by the
actual Fortran standard.

[Bug libfortran/99740] floating point exception in rand() in gfortran

[pvoytas at gmail dot com]

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

--- Comment #5 from Paul A. Voytas <pvoytas at gmail dot com> ---
I understand random_number() is much better--which is great and I'll use it
going forward. I was just trying to not have to recode a lot in cases where I
used rand() because it was good enough. I certainly wouldn't expect bit
equivalent behavior, this just seemed beyond that--since it was included for
backwards compatibility but doesn't quite seem to be functionally so. Thanks
for your time looking into it.