Dear all,
A quick question, based on your experience, which random number generator is the best?
Could you also include the source code or github/gitlab repo?
I mean combination of speed and robustness. It should be able to generator uniform and gaussian random number and array, to say the least.
Thank you very much in advance!
Discussion of what random number generators should be in stdlib, with links to some codes is here. I donāt know if decisions have been made.
I have used random.f90 of Alan Miller to generate normal deviates. That code uses the built-in random_number subroutine for uniform deviates. I have not compared it with other routines. Miller has several codes for uniform deviates. His rand3.f90 code using an algorithm of Knuth may need a fix, because Knuth wrote there was a problem in one of his RNGs.
One could have several versions of a module with a subroutine that generates uniform deviates and verify that simulation results are not sensitive to the RNG chosen.
Note also that there is a vast body of literature on the generation of random numbers. One name in particular is useful: George Marsaglia. Google for his name and you ought to get a fair amount of publications to chew on 
That said, the random number generators in most Fortran compilers are doing a decent enough job, unless you want really large numbers of random numbers or have cryptographic needs.
But in all such matters: there is no best generator, certainly not if you do not specify what your requirements are.
Compiling and running
ziggurat.f90 George Marsagliaās functions for generating random samples from the uniform, normal and exponential distributions. Translated from C
of Alan Miller with gfortran -O3 gives output
Time to generate 10000000 uniform random nos. = 0.02sec.
Time to generate 10000000 normal random nos. = 0.06sec.
last normal deviate = -1.7869397538387441
Time to generate 10000000 exponential random nos. = 0.05sec.
Std. devn. of chi-squared with 999 d.of.f. = 44.70
Values of chi-squared below should be within about 90. of 999.
Uniform distribution, chi-squared = 1009.86 with 999 deg. of freedom
Normal distribution, chi-squared = 1076.53 with 999 deg. of freedom
Exponential distribution, chi-squared = 958.54 with 999 deg. of freedom
Calling Millerās random.f90, the time taken is much higher,
Time to generate 10000000 normal random nos. = 0.42sec.
for the calling program
program xxrandom_normal_time
use random, only: dp, random_normal
implicit none
real(kind=dp) :: t1,t2,x
integer, parameter :: n = 10**7
integer :: i
call cpu_time(t1)
do i=1,n
x = random_normal()
end do
call cpu_time(t2)
write (*,"(a, i10, a, f9.2, a)") " Time to generate ",n," normal random nos. = ",t2-t1,"sec."
print*,"x=",x
end program xxrandom_normal_time
I am running Windows 10 for Processor Intel(R) Core(TM) i3-8350K CPU @ 4.00GHz, 4008 Mhz, 4 Core(s), 4 Logical Processor(s)
See here for a list and discussion of proposals for stdlib: Pseudorandom number generator Ā· Issue #135 Ā· fortran-lang/stdlib Ā· GitHub, @Beliavsky already posted it above.
Please check also stdlib specs, especially this page mentioning a first implementation for pseudo-random number generators.
Also, there are several PRs in stdlib github project waiting for reviewers:
I agree with @Arjen that compiler random_number() generators are good enough for most applications. The Gaussian random deviate generator is not in the Fortran standard. I have an implementation of it here (which is inspired by the recipe of Press et al. 1990). For Multivariate Normal Distribution, this implementation is likely among the fastest algorithms (much better than MATLABās implementation, for example) as it works directly with the Cholesky factor of the covariance matrix of the distribution instead of the covariance matrix. The performance degrades considerably if the covariance matrix is used for random deviate generation. The Cholesky factorization can be readily computed as exmplified here.
I have noticed that the intrinsic random_number() sometimes generates nonsensical values when the seed of the generator is set to some specific values. I do not know any specific rules for what the seed can be, but there does not seem to exist any rules or restrictions in the standard either.
In such cases, when you burn off the first say N=1000 values (which takes 0.04s with gfortran on my PC), are the following values ok?
Exactly, that is what I realized to be the case. But that seems concerning to me, although I have not done anything to fix it. Based on my old internal notes, this seems to be a cross-compiler issue which means it is likely algorithmic rather than a compiler bug (seen in both GFortran and Intel, for example).
I agree.
In my shallow experience, the intrinsic random_number() seems not only gives mediocre performance, and the thing you guys said. But also, perhaps its period is not long enough (perhaps the one in gfortran has a longer period than intelās intrinsic one).
By no means I am saying they are bad.
Just it seems the performance of the intrinsic rng is slightly slower than the one I am using.
The period of say, intelās random_number is not perhaps very long, it is 10^18 or so, see below
The MT rng seems have much longer period,
By the way,
How about the random number below, I always use it instead of the intrinsic one. You could check it with the intrinsic one.
The code below is a minimal working f90 file as an example, can just compile and run.
module random2
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
contains
subroutine ran1(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48_8)-1_8
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult1 = 44485709377909_8
integer(kind=i8), parameter :: m11 = iand(mult1,mask24)
integer(kind=i8), parameter :: m12 = iand(ishft(mult1,-24),mask24)
integer(kind=i8), parameter :: iadd1 = 96309754297_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m12+is2*m11,mask24),24)+is1*m11+iadd1,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran1
subroutine ran2(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48)-1
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult2 = 34522712143931_8
integer(kind=i8), parameter :: m21 = iand(mult2,mask24)
integer(kind=i8), parameter :: m22 = iand(ishft(mult2,-24),mask24)
integer(kind=i8), parameter :: iadd2 = 55789347517_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m22+is2*m21,mask24),24)+is1*m21+iadd2,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran2
end module random2
module random
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
integer(kind=i8), private, parameter :: mask48 = ishft(1_8,48)-1
integer(kind=i8), private, save :: irn = 1_8,irnsave
contains
function randn(n)
! return an array of random variates (0,1)
use random2
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randn
integer(kind=i8) :: i
do i=1,n
call ran1(randn(i),irn)
enddo
return
end function randn
function gaussian(n)
! return an array of gaussian random variates with variance 1
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: gaussian
real(kind=r8), dimension(2*((n+1)/2)) :: rn
real(kind=r8) :: rn1,rn2,arg,x1,x2,x3
real(kind=r8), parameter :: pi=4.0d0*atan(1.0_r8)
integer(kind=i8) :: i
rn=randn(size(rn,kind=i8))
do i=1,n/2
rn1=rn(2*i-1)
rn2=rn(2*i)
arg=2.0_r8*pi*rn1
x1=sin(arg)
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(2*i-1)=x1*x3
gaussian(2*i)=x2*x3
enddo
if (mod(n,2).ne.0) then
rn1=rn(n)
rn2=rn(n+1)
arg=2.0_r8*pi*rn1
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(n)=x2*x3
endif
return
end function gaussian
subroutine setrn(irnin)
! set the seed
integer(kind=i8) :: irnin
integer(kind=i8) :: n
integer(kind=i8), allocatable :: iseed(:)
irn=iand(irnin,mask48)
call savern(irn)
call RANDOM_SEED(SIZE=n)
allocate(iseed(n))
iseed=int(irnin)
call RANDOM_SEED(PUT=iseed) ! set the internal ran function's seed.
return
end subroutine setrn
subroutine savern(irnin)
! save the seed
integer(kind=i8) :: irnin
irnsave=irnin
return
end subroutine savern
subroutine showirn(irnout)
! show the seed
integer(kind=i8) :: irnout
irnout=irnsave
return
end subroutine showirn
function randomn(n) ! the default ramdom_number subroutine.
! return an array of random variates (0,1)
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randomn
call RANDOM_NUMBER(randomn)
return
end function randomn
subroutine random1(x) ! generate random number x using the default ramdom_number subroutine.
real(kind=r8) :: x
call RANDOM_NUMBER(x)
return
end subroutine random1
end module random
module constants
implicit none
integer, public, parameter :: i4=selected_int_kind(9)
integer, public, parameter :: i8=selected_int_kind(15)
integer, public, parameter :: r8=selected_real_kind(15,9)
real(kind=r8), public, parameter :: zero=0.0_r8,one=1.0_r8,two=2.0_r8,three=3.0_r8,four=4.0_r8 &
,five=5.0_r8,six=6.0_r8,seven=7.0_r8,eight=8.0_r8,nine=9.0_r8 &
,ten=10.0_r8,tenth=.1_r8,half=.5_r8,third=1.0_r8/3.0_r8,sixth=1.0_r8/6.0_r8 &
,pi=4.0_r8*atan(1.0_r8) &
,normpdf_factor=1.0_r8/sqrt(8.0_r8*atan(1.0_r8)) ! 1/sqrt(2 pi)
complex(kind=r8), public, parameter :: czero=(0.0_r8,0.0_r8),ci=(0.0_r8,1.0_r8),cone=(1.0_r8,0.0_r8)
end module constants
program main
use random
use random2
use constants
implicit none
integer(kind=i8) :: seed
real (kind=r8) :: x1(1),x100(100)
! use the RNG illustration
seed = 123456789
call setrn(int(seed,8))
! generate one uniform randum number
x1=randn(1)
write (6,*) 'the one random number is', x1(1)
! generate 100 uniform random number array
x100=randn(100)
write (6,*) 'the 100 randum numbers are ', x100
! generate one gaussian randum number
x1=gaussian(1)
write (6,*) ' the one Gaussian random number is', x1(1)
! generate 100 gaussian random number array
x100=gaussian(100)
write (6,*) 'the 100 Gaussian randum numbers are ', x100
stop ('Program end normally.')
end program main
Again, I am not expert in RNG, if you have RNG better than the one I am using now, please advise. Thank you very much indeed!
Thanks for the code. Intel Fortran compiles and runs it with the default options, but
gfortran -std=f2018 says
c:\fortran\test>gfortran -std=f2018 random.f90
random.f90:106:25:
106 | call RANDOM_SEED(SIZE=n)
| 1
Error: 'size' argument of 'random_seed' intrinsic at (1) must be of kind 4
random.f90:109:24:
109 | call RANDOM_SEED(PUT=iseed) ! set the internal ran function's seed.
| 1
Error: 'put' argument of 'random_seed' intrinsic at (1) must be of kind 4
random.f90:88:13:
88 | if (mod(n,2).ne.0) then
| 1
Error: GNU Extension: Different type kinds at (1)
random.f90:148:2:
148 | ,five=5.0_r8,six=6.0_r8,seven=7.0_r8,eight=8.0_r8,nine=9.0_r8 &
| 1
Warning: Nonconforming tab character at (1) [-Wtabs]
random.f90:149:2:
149 | ,ten=10.0_r8,tenth=.1_r8,half=.5_r8,third=1.0_r8/3.0_r8,sixth=1.0_r8/6.0_r8 &
| 1
Warning: Nonconforming tab character at (1) [-Wtabs]
random.f90:150:2:
150 | ,pi=4.0_r8*atan(1.0_r8) &
| 1
Warning: Nonconforming tab character at (1) [-Wtabs]
random.f90:151:2:
151 | ,normpdf_factor=1.0_r8/sqrt(8.0_r8*atan(1.0_r8)) ! 1/sqrt(2 pi)
| 1
Warning: Nonconforming tab character at (1) [-Wtabs]
random.f90:156:9:
156 | use random
| 1
Fatal Error: Cannot open module file 'random.mod' for reading at (1): No such file or directory
compilation terminated.
Line 88 can be changed to if (mod(n,2_i8).ne.0) then .
Do you have references for the algorithms?
Iād also mention the SPRNG library, if you need the random numbers at a large scale, e.g. for Monte Carlo simulations:
http://www.sprng.org/SPRNGmain.html
One issue that I found with the default random_number in GFortran is that it was hard to get exactly the same ārandom numbersā across platforms, for testing purposes. Also I needed something that is very fast. So I ended up writing this:
subroutine lcg_int32(x)
! Linear congruential generator
! https://en.wikipedia.org/wiki/Linear_congruential_generator
integer(int32), intent(out) :: x
integer(int32), save :: s = 26250493
s = modulo(s * 48271, 2147483647)
x = s
end subroutine
subroutine lcg_realsp(x)
real(sp), intent(out) :: x
integer(int32) :: s
call lcg_int32(s)
x = s / 2147483647._sp
end subroutine
subroutine lcg_realdp(x)
real(dp), intent(out) :: x
integer(int32) :: s
call lcg_int32(s)
x = s / 2147483647._dp
end subroutine
Itās probably not very good, but it worked for my testing purposes.
I think at some point I translated the C++ LCG code from Numerical Recipes 3rd ed. and used that for my testing purposes.
MODULE matgen
IMPLICIT NONE
INTEGER, PARAMETER :: dp = KIND(1.D0)
INTEGER, PARAMETER :: i64 = SELECTED_INT_KIND(18)
CONTAINS
REAL(dp) FUNCTION myran(seed)
! This is a pseudorandom number generator adapted from Numerical Recipes 3e
! but with different numbers from the l'Ecuyer tables
INTEGER(i64), INTENT(IN) :: seed
! "Good" numbers from Numerical Recipes
INTEGER, PARAMETER :: A2_1 = 20
INTEGER, PARAMETER :: A2_2 = 41
INTEGER, PARAMETER :: A2_3 = 8
INTEGER(i64), PARAMETER :: D5 = 702098784532940405_i64
REAL(dp), PARAMETER :: bit = 5.42101086242752217D-20
! States
INTEGER(i64), SAVE :: v, w
! Init with seed
v = seed * D5
! Shift and XOR
w = IEOR(w, ISHFT(v, A2_1) )
w = IEOR(w, ISHFT(w, -A2_2) )
w = IEOR(w, ISHFT(w, A2_3) )
! Return double
myran = 0.5D0 + w * bit
END FUNCTION myran
SUBROUTINE genvec(N, V)
! Arguments
INTEGER, INTENT(IN) :: N
REAL(dp), DIMENSION(N), INTENT(OUT) :: V
! Local variables
INTEGER :: i
REAL(dp) :: r, s
! Calculate analytical sum to scale the vector
r = REAL(N,dp)
s = SQRT( r * (r + 1.D0) * (2.D0*r + 1.D0) / 6.0D0 ) ! Analytical sum
DO i = 1, N
V(i) = REAL(i,dp) / s;
END DO
RETURN
END SUBROUTINE genvec
SUBROUTINE genmat(M, N, A)
! Arguments
INTEGER, INTENT(IN) :: M, N
REAL(dp), DIMENSION(M,N), INTENT(OUT) :: A
! Local variables
INTEGER :: i, j
INTEGER(i64), PARAMETER :: seed = 8652445042_i64
DO j = 1, N
DO i = 1, M
A(i,j) = myran(seed)
END DO
END DO
RETURN
END SUBROUTINE genmat
END MODULE matgen
Thanks @certik . I added a few lines to your code to make it a stand-alone module:
module lcg_random_mod
use, intrinsic :: iso_fortran_env, only: int32, real32, real64
implicit none
integer, parameter :: sp = real32, dp = real64
contains
subroutine lcg_int32(x)
! Linear congruential generator
! https://en.wikipedia.org/wiki/Linear_congruential_generator
integer(int32), intent(out) :: x
integer(int32), save :: s = 26250493
s = modulo(s * 48271, 2147483647)
x = s
end subroutine lcg_int32
subroutine lcg_realsp(x)
real(sp), intent(out) :: x
integer(int32) :: s
call lcg_int32(s)
x = s / 2147483647._sp
end subroutine lcg_realsp
subroutine lcg_realdp(x)
real(dp), intent(out) :: x
integer(int32) :: s
call lcg_int32(s)
x = s / 2147483647._dp
end subroutine lcg_realdp
end module lcg_random_mod
Have you tested this against the intrinsic random_number() for performance or quality?
NAG Fortran Compiler employs the MT with the advice to users shown here:
https://www.nag.com/nagware/np/r70_doc/manual/compiler_2_13.html
Sorry. Please see the code below.
Just added some _i8 so gfortran is happier.
gfortran -O3 allinone2.f90 -o test
then
./test
should work.
It is RNG from my PhD advisor, probably box muller. It has been used in many nuclear quantum Monte Carlo code. If you have better ones please advise.
Thank you so much!
module random2
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
contains
subroutine ran1(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48_8)-1_8
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult1 = 44485709377909_8
integer(kind=i8), parameter :: m11 = iand(mult1,mask24)
integer(kind=i8), parameter :: m12 = iand(ishft(mult1,-24),mask24)
integer(kind=i8), parameter :: iadd1 = 96309754297_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m12+is2*m11,mask24),24)+is1*m11+iadd1,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran1
subroutine ran2(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48)-1
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult2 = 34522712143931_8
integer(kind=i8), parameter :: m21 = iand(mult2,mask24)
integer(kind=i8), parameter :: m22 = iand(ishft(mult2,-24),mask24)
integer(kind=i8), parameter :: iadd2 = 55789347517_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m22+is2*m21,mask24),24)+is1*m21+iadd2,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran2
end module random2
module random
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
integer(kind=i8), private, parameter :: mask48 = ishft(1_8,48)-1
integer(kind=i8), private, save :: irn = 1_8,irnsave
contains
function randn(n)
! return an array of random variates (0,1)
use random2
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randn
integer(kind=i8) :: i
do i=1,n
call ran1(randn(i),irn)
enddo
return
end function randn
function gaussian(n)
! return an array of gaussian random variates with variance 1
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: gaussian
real(kind=r8), dimension(2*((n+1)/2)) :: rn
real(kind=r8) :: rn1,rn2,arg,x1,x2,x3
real(kind=r8), parameter :: pi=4.0d0*atan(1.0_r8)
integer(kind=i8) :: i
rn=randn(size(rn,kind=i8))
do i=1,n/2
rn1=rn(2*i-1)
rn2=rn(2*i)
arg=2.0_r8*pi*rn1
x1=sin(arg)
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(2*i-1)=x1*x3
gaussian(2*i)=x2*x3
enddo
if (mod(n,2).ne.0) then
rn1=rn(n)
rn2=rn(n+1)
arg=2.0_r8*pi*rn1
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(n)=x2*x3
endif
return
end function gaussian
subroutine setrn(irnin)
! set the seed
integer(kind=i8) :: irnin
integer(kind=i8) :: n
integer(kind=i8), allocatable :: iseed(:)
irn=iand(irnin,mask48)
call savern(irn)
!call RANDOM_SEED(SIZE=n)
allocate(iseed(n))
iseed=int(irnin)
!call RANDOM_SEED(PUT=iseed) ! set the internal ran function's seed.
return
end subroutine setrn
subroutine savern(irnin)
! save the seed
integer(kind=i8) :: irnin
irnsave=irnin
return
end subroutine savern
subroutine showirn(irnout)
! show the seed
integer(kind=i8) :: irnout
irnout=irnsave
return
end subroutine showirn
function randomn(n) ! the default ramdom_number subroutine.
! return an array of random variates (0,1)
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randomn
call RANDOM_NUMBER(randomn)
return
end function randomn
subroutine random1(x) ! generate random number x using the default ramdom_number subroutine.
real(kind=r8) :: x
call RANDOM_NUMBER(x)
return
end subroutine random1
end module random
module constants
implicit none
integer, public, parameter :: i4=selected_int_kind(9)
integer, public, parameter :: i8=selected_int_kind(15)
integer, public, parameter :: r8=selected_real_kind(15,9)
real(kind=r8), public, parameter :: zero=0.0_r8,one=1.0_r8,two=2.0_r8,three=3.0_r8,four=4.0_r8 &
,five=5.0_r8,six=6.0_r8,seven=7.0_r8,eight=8.0_r8,nine=9.0_r8 &
,ten=10.0_r8,tenth=.1_r8,half=.5_r8,third=1.0_r8/3.0_r8,sixth=1.0_r8/6.0_r8 &
,pi=4.0_r8*atan(1.0_r8) &
,normpdf_factor=1.0_r8/sqrt(8.0_r8*atan(1.0_r8)) ! 1/sqrt(2 pi)
complex(kind=r8), public, parameter :: czero=(0.0_r8,0.0_r8),ci=(0.0_r8,1.0_r8),cone=(1.0_r8,0.0_r8)
end module constants
program main
use random
use random2
use constants
implicit none
integer(kind=i8) :: seed
real (kind=r8) :: x1(1),x100(100)
! use the RNG illustration
seed = 123456789
call setrn(int(seed,8))
! generate one uniform randum number
x1=randn(1_8)
write (6,*) 'the one random number is', x1(1)
! generate 100 uniform random number array
x100=randn(100_i8)
write (6,*) 'the 100 randum numbers are ', x100
! generate one gaussian randum number
x1=gaussian(1_i8)
write (6,*) ' the one Gaussian random number is', x1(1)
! generate 100 gaussian random number array
x100=gaussian(100_i8)
write (6,*) 'the 100 Gaussian randum numbers are ', x100
stop ('Program end normally.')
end program mainThe intrinsic random_number is never used, so that two lines with random_numbers or seeds I commented it.
I set it just for fun, in case I do not use my version, I can use the internal one with the seed I set. But I never used intrinsic one.
See the code here, it should work just fine with gfortran.