What random number generator is the best/fastest or do you use?

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 main