Question about the behaviour in practice of the random_number intrinsic. The standard says pseudo-random numbers in the range 0<=r<1. I believe that in practice we always get random numbers from an equally spaced set in that range. I wonder if, in practice, I can count on it that the spacing in that set is the narrowest possible for the involved floating point type. Equivalently, may I expect (in practice) the largest possible result of random_number (r) for a given precision of the real variable r to be nearest(1.0,-1.0) for arguments of that same precision? (That is the nearest floating point number to 1.0 in the negative direction.)
I do think this is the way most pseudorandom number generators (PRNG) are implemented, but that is not specified or mandated by the fortran standard. This means that, in practice, only a tiny fraction of the possible floating point numbers are returned by random_number().
There are also random number generators that are based on physical processes, such as thermal noise. These typically involve an analog to digital conversion step, so the question of which subset of possible floating point numbers are returned depends on the details of that conversion.
Simple PRNGs work by first generating strings of random bits. Then the appropriate number of bits is extracted from that bit string to form a nonnegative integer value X. A floating point value is then computed from the ratio real(X)/real(M) where M is an upper bound with X<M. Here is a little program that shows the maximum value for 32-bit floating point.
program prng
use, intrinsic :: iso_fortran_env, only: int32, real32
integer(int32), parameter :: M = 2**24, X=M-1 ! 24 bits in a normal floating point number.
real(real32), parameter :: nearx = nearest(1.0,-1.0) ! largest number less than 1.0.
real(real32), parameter :: r = real(X) / real(M)
print '(a,b32.32/8x,a,b32.32)', 'binary: X=', X, 'M=', M
print '(8x,a,b32.32)', 'r=', r
print '(4x,a,b32.32)', 'nearx=', nearx
print *, 'X=', X, 'M=', M, 'r=', r, nearx
end program prng
$ gfortran prng.f90 && a.out
binary: X=00000000111111111111111111111111
M=00000001000000000000000000000000
r=00111111011111111111111111111111
nearx=00111111011111111111111111111111
X= 16777215 M= 16777216 r= 0.999999940 0.999999940
This is 23 stored bits plus one hidden bit, which is the way I think the simple PRNGs work. As you can see from the binary and decimal values for r, the floating point value nearest to 1.0 is generated from the division real(X)/real(M). This means for random X that all of the floating point values associated with the largest exponent for values 0.5≤r<1.0 can be generated, 1/2 of the values for the next smaller exponent, 1/4 of the values for the next smaller exponent, and so on. That is, of course, only a tiny fraction of the possible floating point values.
Some PRNGs generate several random bit patterns using various algorithms and then combine them to produce the final random number. If that all happens for the integer bit pattern, then the above comments still apply, only a tiny fraction of the floating point values will be returned. But if instead separate floating point values are generated and combined, then it is possible for larger fractions of the floating point space to be generated. There are also direct ways to generate all possible floating point values with a uniform distribution, but I’m not aware of any fortran compiler that uses that kind of algorithm. All that I know about use the above real(X)/real(M) approach.
One deficiency of the fortran PRNG intrinsic is that it only returns floating point values. It would be much more useful if it could also return random bit strings in integer variables. Those bit strings are generated internally already anyway, so it is just an issue of allowing the argument to be of type integer.
1 Like
Here is a simple test of random_number using two Fortran compilers I use.
These tests indicate:
The returned value is in the range ( 0 <= random_number(x) < 1 )
Silverfrost FTN95 is a reproduceable sequence, while
Gfortran 15.1.0 uses a different starting value ( for the compile options I used)
For my Monte-Carlo simulations I mostly used a reproduceable sequence.
implicit none
integer :: n,k,i
real*4 :: x, mi_x, ma_x
open (unit=11, file='random.log', position='append')
write (11,*) 'Vern : ',compiler_version ()
write (11,*) 'Opts : ',compiler_options ()
n = 2
mi_x = 1
ma_x = 0
do k = 1,24
do i = 1,n
call random_number (x)
if ( x == 0. ) write (11,*) 'zero detected', i
if ( x == 1. ) write (11,*) 'one detected', i
if ( x > ma_x ) ma_x = x
if ( x < mi_x ) mi_x = x
end do
write (11,fmt='(i8,2g18.8)') n,mi_x, ma_x
n = n*2
end do
end
Vern : FTN95 v9.20
Opts : 64;CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
2 0.33098495 0.67307037
4 0.89860201E-01 0.74144781
8 0.55308938E-02 0.75059688
16 0.55308938E-02 0.93005449
32 0.55308938E-02 0.98703653
64 0.55308938E-02 0.98703653
128 0.55308938E-02 0.99690789
256 0.11609793E-02 0.99690789
512 0.10038018E-02 0.99690789
1024 0.78243017E-03 0.99844843
2048 0.78243017E-03 0.99930221
4096 0.16486645E-03 0.99948233
8192 0.69141388E-04 0.99998528
16384 0.28550625E-04 0.99998528
32768 0.97155571E-05 0.99998528
65536 0.35762787E-05 0.99999893
131072 0.35762787E-05 0.99999893
262144 0.35762787E-05 0.99999952
524288 0.11920929E-06 0.99999952
1048576 0.11920929E-06 0.99999952
2097152 0.11920929E-06 0.99999958
4194304 0.11920929E-06 0.99999970
zero detected 7634751
8388608 0.0000000 0.99999988
zero detected 2906636
16777216 0.0000000 0.99999994
Vern : GCC version 15.1.0
Opts : -cpp -iprefix C:/Program Files (x86)/gcc_eq/gcc_15.1.0/bin/../lib/gcc/x86_64-w64-mingw32/15.1.0/ -U_REENTRANT -mtune=generic -march=x86-64 -fno-underscoring -fdollar-ok
2 0.48348904E-02 0.69568729
4 0.48348904E-02 0.79884946
8 0.48348904E-02 0.85110766
16 0.59008598E-04 0.92572713
32 0.59008598E-04 0.95438445
64 0.59008598E-04 0.99653119
128 0.59008598E-04 0.99690586
256 0.59008598E-04 0.99839717
512 0.59008598E-04 0.99957311
1024 0.59008598E-04 0.99957311
2048 0.59008598E-04 0.99957311
4096 0.59008598E-04 0.99957311
8192 0.55015087E-04 0.99961603
16384 0.30994415E-05 0.99991792
32768 0.30994415E-05 0.99993944
65536 0.30994415E-05 0.99999046
131072 0.16093254E-05 0.99999100
262144 0.16093254E-05 0.99999577
524288 0.16093254E-05 0.99999577
1048576 0.95367432E-06 0.99999845
2097152 0.59604645E-06 0.99999887
zero detected 3443890
4194304 0.0000000 0.99999952
8388608 0.0000000 0.99999994
zero detected 811674
zero detected 4808671
16777216 0.0000000 0.99999994
Vern : FTN95 v9.20
Opts : 64;CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
2 0.33098495 0.67307037
4 0.89860201E-01 0.74144781
8 0.55308938E-02 0.75059688
16 0.55308938E-02 0.93005449
32 0.55308938E-02 0.98703653
64 0.55308938E-02 0.98703653
128 0.55308938E-02 0.99690789
256 0.11609793E-02 0.99690789
512 0.10038018E-02 0.99690789
1024 0.78243017E-03 0.99844843
2048 0.78243017E-03 0.99930221
4096 0.16486645E-03 0.99948233
8192 0.69141388E-04 0.99998528
16384 0.28550625E-04 0.99998528
32768 0.97155571E-05 0.99998528
65536 0.35762787E-05 0.99999893
131072 0.35762787E-05 0.99999893
262144 0.35762787E-05 0.99999952
524288 0.11920929E-06 0.99999952
1048576 0.11920929E-06 0.99999952
2097152 0.11920929E-06 0.99999958
4194304 0.11920929E-06 0.99999970
zero detected 7634751
8388608 0.0000000 0.99999988
zero detected 2906636
16777216 0.0000000 0.99999994
Vern : GCC version 15.1.0
Opts : -cpp -iprefix C:/Program Files (x86)/gcc_eq/gcc_15.1.0/bin/../lib/gcc/x86_64-w64-mingw32/15.1.0/ -U_REENTRANT -mtune=generic -march=x86-64 -fno-underscoring -fdollar-ok
2 0.67090213E-01 0.50203764
4 0.40025771E-01 0.96130651
8 0.11717677E-01 0.96130651
16 0.11717677E-01 0.96130651
32 0.11717677E-01 0.99553567
64 0.11436999E-01 0.99949199
128 0.55414438E-02 0.99949199
256 0.34532547E-02 0.99949199
512 0.21323562E-02 0.99949199
1024 0.45263767E-03 0.99949199
2048 0.40495396E-03 0.99982810
4096 0.33193827E-03 0.99982810
8192 0.63836575E-04 0.99998480
16384 0.39398670E-04 0.99998719
32768 0.14841557E-04 0.99999207
65536 0.71525574E-05 0.99999207
131072 0.11920929E-06 0.99999797
262144 0.11920929E-06 0.99999797
524288 0.11920929E-06 0.99999964
zero detected 859344
1048576 0.0000000 0.99999964
2097152 0.0000000 0.99999964
4194304 0.0000000 0.99999982
8388608 0.0000000 0.99999994
zero detected 3340268
16777216 0.0000000 0.99999994
I found this fairly recent article on PRNGs that is an interesting read. Generating Random Floating-Point Numbers by Dividing Integers: A Case Study - PMC This article points out that when the integer division is used to generate the floating point values, only a small subset of the possible values are returned. One of the features discussed is that expressions like real(X)/real(M) where an actual division is executed are not necessarily the same due to rounding as expressions like real(X)*RM where RM=1.0/real(M). Of course, in fortran, the programmer cannot necessarily control that evaluation since the compiler is allowed to optimize the former to the latter. And whether a compiler does that might even depend on compiler options, the context of the expressions, or if the expression is evaluated at compile time or at run time. I added that expression to my previous fortran code, and both expressions produced the same result, but I did not experiment with compiler options or with different compilers, so maybe it is just luck when the correctly rounded result is produced.
The standard does not specify the quality level of the RANDOM_NUMBER implementation, and there are variations in algorithms across different compilers. One also needs to be aware of the resolution of the return data type as a single-precision real typically has at most 23 fraction bits, and there’s nothing to say that using higher-precision arguments will give finer resolution.
As for sequences, earlier editions were silent as to whether you’d get the same or different sequence on each run, or across images. Fortran 2018 added the new intrinsic subroutine RANDOM_INIT to deal with this.
3 Likes
I understand why this was done in legacy times, but I wonder if that should be reconsidered now. In fortran there are several issues that should be considered, including quality (e.g. as you point out, 23 bits vs. 53 bits, etc.), but also for parallel code whether the PRNGs generated on one image is distinct, or can partially or entirely overlap, with the PRNGs on other images. Currently, if a programmer wants to avoid overlap of values and still have robust random numbers, then he must avoid the fortran intrinsic entirely and use his own code. I wonder if it is time to ask a little more from the fortran standard in this respect?
The distinct sequences across images is handled by RANDOM_INIT. We had extensive discussions on these issues and the consensus was that programmers who really cared about RNG quality would use their own chosen implementation, especially as there are choices to be made about distribution. The intrinsic is “good enough” for most uses.
1 Like
As you stated originally, the standard does not specify the quality of the intrinsic. I’m suggesting that now, in these modern times, it might be time to do that. I am a fan of some of the features of the current PRNG used by gfortran. It has a very long cycle, and it has the capability to skip forward in that sequence by a large number of steps. Those two features combined allow the sequences generated on different images, even millions of images, to still be very long but with no overlap at all. That might be a good candidate algorithm for the fortran standard to specify. The committee might want to make it abstract somehow, instead of specifying the algorithm they might just specify the features (say a minimum cycle length, and the ability to skip forward by arbitrary numbers of steps), but that is what I’m thinking.
1 Like
I have always found this to be the case.
In my simulation studies, I found any limitations in random_number could be easily replaced by revising the approach. Rather that testing for rare events (which can focus only on the extremes of the random values), it can be better to estimate the time till the next event, using a Poisson distribution (which is easily sampled from the full range of random_number). My first studies were done successfully using my random number generator based on 2 integer*2 primes.
RANDOM_INIT is a good improvement, as having a repeatable generated sequence was essential for program testing.
That’s correct but it also means that of questions asked by OP,
- the largest random number should indeed be
nearest(1.0,-1.0) - the spacing should be equal to
spacing(0.75)over the entire 0.0<=x<1.0 interval
Yes, I agree, this is correct for PRNGs computed using the division expression. For real32 floating point, that spacing is 5.96e-8. That also means that none of the floating point numbers between 0.0 and 5.96e-8 will be returned by the PRNG. The binary exponent for that value is -23. The exponent for tiny(1.0) is -125, so there are 102 exponents that are basically ignored. There are (2^24-1) fractions for each of those exponents, so there are 102*(2^24-1) floating point numbers that will not be returned. So if I’m treating the hidden bit and the normal exponents correctly, there are 126*(2^24-1) total floating point numbers in the range 0.0≤r<1.0, and only 2^24-1 of them are returned by a PRNG that uses the division expression. That is only 1/126 or 0.79% of the total candidate floating point values.
But if the full set were to be considered, then there would be a large emphasis on very small values, wouldn’t there? The random_number routine is supposed to return a uniformly distributed set.
2 Likes
If you have a source of random bits, it is most convenient to return a binary floating-point number in [1.0 , 2.0), as they are all in the same binade, meaning they share sign and exponent bits and have uniformly random bits in the mantissa.
1 Like
I guess this gets back to the definition of “uniformly distributed”? It does not mean that the set of all possible floating point numbers is equally sampled – that distribution would indeed have more small numbers (<0.5) than large numbers (>0.5). Rather, it means that the numbers in [0.25,0.5) are sampled half as often as those in [0.5,1.0), and those in [0.125,0.25) are sampled half as often as those in [0.25,0.5), and so on. Each binary exponent should be sampled half as often as the next larger exponent. In principle, that can continue all the way down to the exponent for tiny(1.0), and within each exponent, the full set of fractions could be equally sampled. Such a PRNG could then return the full set of normal floating point numbers from 0.0 up to nearest(1.0,-1.0) with uniform distribution. As I said previously, I do not know of any fortran compiler that implements such a PRNG for its random_number() intrinsic. The reason, as discussed in this thread, is that the simpler division approach, which return only 1/126 of those numbers, is “good enough” for most applications.
I experimented a little with this idea of uniform distributions, and here is what I came up with.
module random_mod
use, intrinsic :: iso_fortran_env, only: int32, real32, int64, real64
interface random
module procedure random_s
module procedure random_d
end interface random
contains
function xoshiro256ss() result (res)
! return an integer with random bits.
! based on fortran stdlib xoshiro256ss() at
! https://stdlib.fortran-lang.org/sourcefile/stdlib_random.fypp.html
implicit none
integer(int64) :: res, t
integer(int64), parameter :: si = 614872703977525537_int64 ! default seed value
integer(int64), save :: st(4) = si ! internal state.
intrinsic :: shiftl, ieor, ishftc
res = ishftc(st(2) * 5, 7) * 9
! update the internal state for the next call.
t = shiftl(st(2), 17)
st(3) = ieor(st(3), st(1))
st(4) = ieor(st(4), st(2))
st(2) = ieor(st(2), st(3))
st(1) = ieor(st(1), st(4))
st(3) = ieor(st(3), t)
st(4) = ishftc(st(4), 45)
return
end function xoshiro256ss
impure elemental subroutine random_s( x )
! return a random number x with uniform distribution from the set
! of all possible floating point numbers in the range [0.0,1.0).
implicit none
real(real32), intent(out) :: x
integer(int32) :: e, ix, lz
integer, parameter :: emin = minexponent(x)
intrinsic :: shiftr, leadz, mvbits, transfer, int
! ieee real32 format is: s (31), e (30...23), f (22...0)
ix = int( shiftr( xoshiro256ss(), 41 ), int32) ! throw away the low bits for the fraction.
! compute a random integer in [0,-emin] with a geometric distribution: 1/2, 1/4, 1/8...
lz = leadz( xoshiro256ss() )
if ( lz == 64 ) lz = lz + leadz( xoshiro256ss() )
if ( lz > -emin ) then
x = 0.0_real32
else
e = 1 - emin - lz ! biased exponent in [1,1-emin].
call mvbits( e, 0, 8, ix, 23 )
x = transfer( ix, x )
endif
return
end subroutine random_s
impure elemental subroutine random_d( x )
! return a random number x with uniform distribution from the set
! of all possible floating point numbers in the range [0.0,1.0).
implicit none
real(real64), intent(out) :: x
integer(int64) :: e, ix
integer :: i, lz, lzi
integer, parameter :: emin = minexponent(x)
integer, parameter :: numw = ((-emin-1) / bit_size(ix) + 1) ! number of int64 words needed to hold -emin bits.
intrinsic :: shiftr, leadz, mvbits, transfer
! ieee real64 format is: s (63), e (62...52), f (51...0)
ix = shiftr( xoshiro256ss(), 12 ) ! throw away the low bits for the fraction.
! compute a random integer in [0,-emin] with a geometric distribution: 1/2, 1/4, 1/8...
lz = 0
do i = 1, numw
lzi = leadz( xoshiro256ss() )
lz = lz + lzi
if ( lz < 64 ) exit
enddo
if ( lz > -emin ) then
x = 0.0_real64
else
e = 1 - emin - lz ! biased exponent in [1,1-emin].
call mvbits( e, 0, 11, ix, 52 )
x = transfer( ix, x )
endif
return
end subroutine random_d
end module random_mod
program prng
use random_mod
implicit none
integer :: i, emin, emax, imax=10**8
real(real32) :: x
real(real64) :: d, sum, mean, xmin, xmax
character(*), parameter :: cfmt = '(*(g0.9:1x))'
emin = huge(emin)
emax = -huge(emax)
xmin = huge(xmin)
xmax = -1.0
sum = 0.0
do i = 1, imax
call random( x )
d = x
emin = min( emin, exponent(d) )
emax = max( emax, exponent(d) )
xmin = min( xmin, d )
xmax = max( xmax, d )
sum = sum + d
enddo
mean = sum / imax
print cfmt, 'imax=', imax, 'mean=', mean, 'xmin=', xmin, &
& 'xmax=', xmax, 'emin=', emin, 'emax=', emax, 'emin(tiny)=', minexponent(tiny(x))
emin = huge(emin)
emax = -huge(emax)
xmin = huge(xmin)
xmax = -1.0
sum = 0.0
do i = 1, imax
call random( d )
emin = min( emin, exponent(d) )
emax = max( emax, exponent(d) )
xmin = min( xmin, d )
xmax = max( xmax, d )
sum = sum + d
enddo
mean = sum / imax
print cfmt, 'imax=', imax, 'mean=', mean, 'xmin=', xmin, &
& 'xmax=', xmax, 'emin=', emin, 'emax=', emax, 'emin(tiny)=', minexponent(tiny(d))
end program prng
When executed with gfortran, here is the output:
$ gfortran prng.f90 && a.out
imax= 100000000 mean= 0.499983841 xmin= 0.151877551E-7 xmax= 0.999999940 emin= -25 emax= 0 emin(tiny)= -125
imax= 100000000 mean= 0.500022983 xmin= 0.822534641E-10 xmax= 0.999999998 emin= -33 emax= 0 emin(tiny)= -1021
I think the exponents are handled correctly for both real32 and real64, but beware of off-by-1 errors. First, notice that I’m generating int64 random numbers using the fortran stdlib routine xoshiro256ss(). I simplified it for this posting, but that is where that code came from.
Then from those random bit strings, I just use bit operators to load in the exponent and fraction bits. So in both cases, this should result in the full set of floating point values in [0.0,1.0) being sampled if called enough times. Because the real number is constructed in this way, there is absolutely no chance of a number x>=1.0 being returned or of a subnormal number being returned. There are no floating point operations, so no chance of overflows or underflows occurring. The spacing of the numbers tracks with the exponent, so the smallest spacing will be spacing(tiny(x)) when x has that same exponent. But if you look at the geometric distribution of the exponents, the chance of such small numbers is very small.
However, it occurred to me that it would be extremely unlikely for an exact 0.0 value to be returned. Certainly for the real64 case, but also even for the real32 case. If you compare that to the results for the simple division expression, where 0.0 is returned relatively often, that seems odd. I wonder which approach actually gives the best statistical results. I only test the mean here, so nothing sophisticated beyond that.
I’m generating that geometric distribution for the exponent in a brute force way with repeated calls to xoshiro256ss(). I suspect that I could do that with a single call by reusing the same random bits for both the fraction and for the exponent calculation; only rarely would the second call (or beyond) be necessary for the exponent. However, I was unsure about that reuse, so in this code I took the safe option with repeated xoshiro256ss() calls.
I use the equivalence hack here to put the bits into the real result. This could also be done with transfer(), of course, but that seems rather inelegant here for some reason. However, at least in the real32 case, there is still the big/little-endian dependence on the hardware. This code should work for little-endian addressing and for IEEE floating point formats. [edit: in the updated code, I removed the equivalence hack and used transfer(). This version should compile on all fortran compilers.]
Any comments on this code would be welcome.
edit: I did have an off-by-one error in the exponent that allowed subnormal return values. I modified the fortran code to correct this error.
1 Like
I’m not claiming to be an expert on random numbers, but I did do some work in this area a few years ago:
Specifically the way I constructed the tests may be of interest, which is to take a large number of samples from the RNG, then split the samples into bins and see whether the expected number of samples were in each bin.
3 Likes
Did you find any discrepancies, or did all implementations pass this test?
That already puts you beyond my experience in this field. Maybe you can answer this question I have about returning an exact 0.0? In the division-based codes, an exact 0.0 will be returned fairly often. If there are 24 bits in the fraction (e.g. 23 plus 1 hidden bit), then an exact 0.0 will be returned on average once every (2^24-1) calls, or roughly once every 16.7 million calls. That was observed already in the post Spacing of the results of Fortran intrinsic random_number - #3 by JohnCampbell where multiple 0.0 return values were detected in the outputs.
However, if you look at the code I posted, then an exact 0.0 will be returned only when the binary exponent is below -125. In the computed geometric distribution of those exponents, that would be expected to occur only once every 2^126 of the time, or roughly once every 8.5e+37. Practically speaking, that should be never, but I left the possibility in my code anyway.
So once every 16.7 million is a lot different than once every 8.5e+37. What consequences does that difference have for, say, monte carlo simulations? Is that difference something that is detected by the standard random number tests?
1 Like
I expanded my test program using 2^26 samples of RANDOM_NUMBER to identify :
- the extremes of returned random values, (min, 1-max)
- the uniform distribution of values using 10x10% bins (plus 2 extreme1% bins),
- the repeated special values (zero, one and first value)
I tested two compilers I use (Silverfrost FTN95 and Gfortran 15.1).
Both indicated “good” uniform distribution across 10 x 10% bins. (plus 2 x 1% special bins at the extremes) I report the variation from the expected number of values for each bin; “good” is an observation, not a statistical test, std_dev could be used?
For both compilers, multiple zero values were returned (about 5 in 2^26)
Silverfrost provides the same sequence for each run, while Gfortran by default provides a unique sequence for each run.
For one of 6 Gfortran tests, the first random value was returned 8 times in 2^25 samples which looks unusual.
These results show the two random_number intrinsics tested are very good for the monti-carlo studies I have carried out.
use iso_fortran_env
implicit none
integer :: n,k,i
real*4 :: x, mi_x, ma_x, first
integer :: num(0:11), exp(0:11),j
open (unit=11, file='random.log', position='append')
write (11,*) 'Vern : ',compiler_version ()
write (11,*) 'Opts : ',compiler_options ()
write (11,fmt='("small: ",g18.8)') epsilon(x)
call random_number (first)
n = 2
mi_x = 1
ma_x = 0
num = 0
do k = 1,25
do i = 1,n
call random_number (x)
if ( x == first ) then
write (11,*) 'first detected', i
cycle
end if
if ( x == 0. ) then
write (11,*) 'zero detected', i
cycle
end if
if ( x == 1. ) then
write (11,*) 'one detected', i
cycle
end if
if ( x > ma_x ) ma_x = x
if ( x < mi_x ) mi_x = x
j = (x*10.)+1
if ( x < 0.01 ) j = 0
if ( x > 0.99 ) j = 11
num(j) = num(j) + 1
end do
i = sum(num(0:11))/10
exp(0) = i/10 ; exp(11) = exp(0)
exp(1) = i - i/10 ; exp(10) = exp(1)
exp(2:9) = i
write ( *,11) n,mi_x, 1.-ma_x, i, (num(j)-exp(j),j=0,11)
write (11,11) n,mi_x, 1.-ma_x, i, (num(j)-exp(j),j=0,11)
11 format ( i8, 2g18.8, 13(1x,i0) )
n = n*2
end do
end
Two (of many) recent runs are:
Vern : GCC version 15.1.0
Opts : -cpp -iprefix C:/Program Files (x86)/gcc_eq/gcc_15.1.0/bin/../lib/gcc/x86_64-w64-mingw32/15.1.0/ -U_REENTRANT -mtune=generic -march=x86-64 -fno-underscoring -fdollar-ok
small: 0.11920929E-06
2 0.40233296 0.10071009 0 0 0 0 0 0 1 0 0 0 1 0 0
4 0.77126265E-01 0.10071009 0 0 1 1 1 0 2 0 0 0 1 0 0
8 0.59891820E-01 0.93671441E-01 1 0 1 1 0 -1 3 -1 0 -1 2 0 0
16 0.13292909E-01 0.23644745E-01 3 0 2 2 -1 -1 3 -3 -2 0 1 -1 0
32 0.13292909E-01 0.37516356E-02 6 0 0 4 -1 -1 2 -2 -2 0 1 -1 2
64 0.13292909E-01 0.67180395E-03 12 -1 -2 3 -3 -3 4 0 0 2 7 -3 2
128 0.13292909E-01 0.67180395E-03 25 -2 5 4 0 -7 5 -3 4 3 4 -10 1
256 0.41491389E-02 0.67180395E-03 51 -1 4 5 7 0 6 -11 7 0 -8 -8 -1
512 0.81789494E-03 0.67180395E-03 102 3 7 7 -9 4 15 -22 8 5 7 -18 -5
1024 0.81789494E-03 0.67180395E-03 204 -1 16 6 6 -4 18 -32 7 3 8 -16 -5
2048 0.18888712E-03 0.18984079E-03 409 4 4 -24 3 4 43 -50 19 -12 21 -7 -1
4096 0.18888712E-03 0.18131733E-03 819 -2 -23 -10 33 -17 37 -23 29 -29 -3 1 7
8192 0.51140785E-04 0.94771385E-04 1638 -10 -22 9 66 -10 34 -14 20 -46 -22 -8 5
16384 0.44226646E-04 0.11980534E-04 3276 7 -2 39 72 5 55 -79 -48 8 -31 -32 12
32768 0.28014183E-04 0.11980534E-04 6553 8 -25 -127 169 -28 76 -63 -81 32 -12 76 -21
65536 0.12576580E-04 0.11980534E-04 13107 46 -47 -235 289 49 -39 -187 -38 47 -20 164 -29
131072 0.12576580E-04 0.35762787E-06 26214 71 -276 -159 313 55 -95 -84 43 175 -15 48 -74
262144 0.42915344E-05 0.35762787E-06 52428 78 -26 -95 386 97 -75 24 -129 -97 -10 -75 -72
524288 0.35762787E-06 0.35762787E-06 104857 113 251 -97 348 107 -268 -53 -125 -363 145 7 -61
1048576 0.35762787E-06 0.35762787E-06 209715 19 -75 -212 614 -143 368 -108 -42 -555 221 -7 -80
zero detected 76215
2097152 0.59604645E-07 0.35762787E-06 419430 -164 -508 -840 1340 -95 -299 732 -476 -170 522 -33 -8
zero detected 2307888
4194304 0.59604645E-07 0.23841858E-06 838860 -443 -446 -552 1874 375 -226 878 -637 -463 -190 -118 -48
first detected 1718498
first detected 6104323
8388608 0.59604645E-07 0.59604645E-07 1677721 -723 -38 -25 1204 -264 -419 3046 -1481 -576 923 -1422 -225
first detected 6159289
16777216 0.59604645E-07 0.59604645E-07 3355442 -495 -1110 -1515 2274 472 -598 3814 -1603 -2110 1733 -914 57
zero detected 6091588
zero detected 11099072
zero detected 11649354
zero detected 14245640
zero detected 18043806
first detected 23185239
zero detected 26210469
zero detected 27926720
33554432 0.59604645E-07 0.59604645E-07 6710884 -1747 -2815 928 1380 262 203 1826 -1930 -1308 3476 -650 384
Vern : FTN95 v9.20
Opts : 64;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
small: 0.11920929E-06
2 0.14748472 0.66901505 0 0 0 1 0 1 0 0 0 0 0 0 0
4 0.89860201E-01 0.25855219 0 0 1 1 0 2 0 0 0 2 0 0 0
8 0.55308938E-02 0.24940312 1 1 0 0 -1 1 1 1 1 2 -1 -1 0
16 0.55308938E-02 0.69945514E-01 3 1 1 -1 -2 1 0 0 0 1 1 -2 0
32 0.55308938E-02 0.12963474E-01 6 1 0 0 -4 2 0 2 2 2 0 -3 0
64 0.55308938E-02 0.12963474E-01 12 0 -4 -2 -3 0 0 5 4 4 3 0 -1
128 0.55308938E-02 0.30921102E-02 25 0 -1 -6 -4 1 -5 4 5 9 1 0 0
256 0.11609793E-02 0.30921102E-02 51 1 0 -4 -2 -1 -7 3 9 12 -6 -3 -2
512 0.10038018E-02 0.30921102E-02 102 1 7 -2 -7 4 -11 0 9 9 1 -9 0
1024 0.78243017E-03 0.15515685E-02 204 7 9 -2 -12 17 -27 5 5 14 3 -13 0
2048 0.78243017E-03 0.69779158E-03 409 10 14 -5 -20 11 -20 -2 -17 42 19 -26 -2
4096 0.16486645E-03 0.51766634E-03 819 7 -25 -52 3 11 -31 -7 9 58 47 -12 -8
8192 0.69141388E-04 0.14722347E-04 1638 11 -10 -18 18 5 -27 -49 -15 66 24 4 -7
16384 0.28550625E-04 0.14722347E-04 3276 34 24 -53 2 -54 1 -93 11 59 23 57 -5
32768 0.97155571E-05 0.14722347E-04 6553 18 -2 -19 -14 -16 -68 -18 -22 27 40 74 4
65536 0.35762787E-05 0.10728836E-05 13107 -7 -79 -37 62 59 -105 25 -79 -111 75 175 22
131072 0.35762787E-05 0.10728836E-05 26214 -34 -99 -42 -25 -47 -111 172 -113 -128 168 263 -2
262144 0.35762787E-05 0.47683716E-06 52428 -17 -29 -21 21 -216 -20 -45 -53 -245 298 344 -11
524288 0.11920929E-06 0.47683716E-06 104857 15 -388 -65 134 -35 -209 4 -134 -344 723 173 130
1048576 0.11920929E-06 0.47683716E-06 209715 66 -230 -274 124 -187 -281 571 -210 -593 747 47 220
2097152 0.11920929E-06 0.41723251E-06 419430 380 -17 414 282 573 -1144 625 -1024 -1240 392 395 366
4194304 0.11920929E-06 0.29802322E-06 838860 407 -961 -208 903 615 -1660 269 -1169 -1098 1038 1183 687
first detected 6209467
zero detected 7634750
8388608 0.59604645E-07 0.11920929E-06 1677721 620 -1967 585 1282 1049 -1871 791 -1496 -1149 103 1140 915
zero detected 2906635
16777216 0.59604645E-07 0.59604645E-07 3355442 1105 -1997 -199 469 817 -361 -516 -176 -675 419 627 494
first detected 2294863
zero detected 4297299
zero detected 19307731
zero detected 25049961
33554432 0.59604645E-07 0.59604645E-07 6710885 689 -684 318 -2167 1101 -1505 291 1343 -1695 2435 55 -176
I replaced the intrinsic random_number() call in @JohnCampbell code with my random() call, and this is the output:
small: 0.11920929E-06
2 0.58845115 0.37459946 0 0 0 0 0 0 0 1 1 0 0 0 0
4 0.56399345E-01 0.19467908 0 0 1 0 0 0 1 1 1 1 1 0 0
8 0.56399345E-01 0.87157488E-02 1 0 1 0 -1 -1 2 1 1 0 1 -1 1
16 0.55183452E-01 0.87157488E-02 3 0 1 1 -1 -2 2 0 1 -2 0 -1 1
32 0.67169042E-02 0.60788393E-02 6 1 3 2 -1 -4 4 1 0 -5 0 -1 2
64 0.61297151E-02 0.55150986E-02 12 1 3 2 -2 -5 1 1 2 0 -1 2 2
128 0.61297151E-02 0.51218271E-03 25 0 5 6 -2 -11 -4 -2 5 -2 1 6 2
256 0.36427409E-02 0.51218271E-03 51 2 0 6 -5 -18 -5 -9 8 0 8 9 4
512 0.36427409E-02 0.51218271E-03 102 1 -5 -1 4 -26 -1 -5 7 8 5 11 4
1024 0.46532866E-03 0.35923719E-03 204 9 -12 -3 -2 -23 0 -9 11 20 6 4 5
2048 0.46532866E-03 0.18090010E-03 409 3 5 -29 1 -8 -19 8 -4 39 -14 5 17
4096 0.25507278E-03 0.18090010E-03 819 3 -14 -17 16 -29 -31 22 -15 75 -37 9 18
8192 0.18587007E-03 0.17511845E-03 1638 5 -17 45 44 -50 -34 27 17 31 -51 -17 2
16384 0.11991579E-04 0.97453594E-04 3276 -27 21 72 90 -136 3 20 30 66 -115 -21 3
32768 0.10942373E-04 0.41961670E-04 6553 -49 -20 111 104 -119 -31 13 -99 32 25 52 -15
65536 0.13819766E-05 0.89406967E-05 13107 -27 -25 160 164 -206 -177 52 -261 35 49 244 -8
131072 0.72946484E-07 0.11920929E-05 26214 35 101 220 -189 -70 -194 14 -232 44 37 258 -22
262144 0.72946484E-07 0.11920929E-05 52428 21 172 169 -325 116 -191 3 -114 45 198 -135 47
524288 0.72946484E-07 0.11920929E-05 104857 78 186 294 -606 499 -73 -147 -180 153 27 -218 -9
1048576 0.72946484E-07 0.59604645E-06 209715 92 97 514 -136 226 -348 258 -402 13 -48 -210 -56
2097152 0.72946484E-07 0.35762787E-06 419430 -119 342 624 -756 310 -1226 941 -237 -88 -20 58 173
4194304 0.11625332E-07 0.35762787E-06 838860 -99 855 677 -111 -31 -1463 333 -1034 98 109 508 164
8388608 0.11625332E-07 0.11920929E-06 1677721 31 1953 657 226 -593 -2341 -372 -1099 1764 -152 252 -322
16777216 0.11625332E-07 0.59604645E-07 3355443 -82 3442 1331 -448 -193 -2034 1210 -3774 1729 424 217 -1822
first detected 12006649
first detected 29178890
33554432 0.10724137E-07 0.59604645E-07 6710886 1009 2629 709 -1017 -1756 2727 -1494 -5605 1598 349 2266 -1415
My random() code did result in two repetitions of the first value, which I guess is to be expected with large numbers of calls, but it does not report any 0.0 values returned, which is consistent with the above discussion. I ran the code with gfortran, flang, and nagfor, and all three compilers give exactly the same output.