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

["sgk at troutmask dot apl.washington.edu" <gcc-bugzilla@gcc.gnu.org>]

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.