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

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.

3 Likes

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 :slight_smile:
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.

4 Likes

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)

4 Likes

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.

1 Like

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:

1 Like

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.

2 Likes

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?

2 Likes

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).

1 Like

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

1 Like

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.

3 Likes

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

2 Likes

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
3 Likes

Have you tested this against the intrinsic random_number() for performance or quality?

1 Like

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

2 Likes

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

The 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.