A Good Fortran Random Number Generator

Question

The GNA of FORTRAN (Rand()) seems to be bad, this because it is proved to be worse than very simple random number generators. For example, in my simulations the GNA below

SUBROUTINE GNA(iiseed)
USE Variaveis
parameter (ia=843314861,ib=453816693,m=1073741824, r231=1./2147483648.)
INTEGER :: iiseed

iiseed = ib + ia*iiseed
if (iiseed.lt.0) iiseed = (iiseed+m) + m
RndNum = iiseed*r231

END SUBROUTINE GNA

has shown itself better than the Rand() of FORTRAN, and this GNA is very simple. Someone would indicate me a good random number generator (GNA) in FORTRAN. Something that combines computational time efficiency and randomness of the generated numbers.

Answer 1

There are several algorithms for number generation pseudo random and the choice of the ideal algorithm depends very much on the objective for which the numbers will be generated.

The most widely used algorithm is the Mersenne Twister - Wikipedia which, according to the article, has, among other advantages, a long period (2^19937 − 1), so the "randomness" of the numbers is better.

If you prefer a faster algorithm, you can use the Xorshift - Wikipedia (examples of implementation are in C language), or some of its variations.

Another option to use is to previously obtain a list of random numbers in a service such as the Random.ORG and use these numbers in your program by reading them from a file.

There are other services like Random.org (and other algorithms as well), which you can refer to on Wikipedia - List of Random Number Generators.

Here are some sites (in English) that already provide the source code for Prngs in FORTRAN:

Alan Miller’s Fortran Software

Multiple stream Mersenne Twister PRNG

Mersenne Twister Home Page

Source Codes in Fortran90 and F90_RANDOM

Random number Generation

Answer 2

Hello, all right? So I use the Marsaglia. It has a very large period and very small running time. They’ve run a lot of tests on him for correlations, and he’s within range. For used, in your program or routine that will use it, just put it in the USE random_utilities header and, before starting, call the subroutine CALL initrandom(i,it). If you need to use the same Seed for comparison, you can call rmarin(ij, kl) where 0 <= IJ <= 31328 & 0 <= KL <= 30081.

MODULE varitype1
    INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
    INTEGER, PARAMETER :: k15 = SELECTED_INT_KIND(15)
    INTEGER, PARAMETER :: DP = KIND(1.0D0)
END MODULE
MODULE raset
    USE varitype1
    IMPLICIT NONE
    REAL(DP) :: U(1:97), cc, cd, cm
    INTEGER(k15) ::  i97, j97 
END MODULE
MODULE random_utilities
    USE varitype1
    IMPLICIT NONE
    INTEGER(I4B) :: rand_inst = 0
    INTEGER(I4B) :: rand_feedback = 1

CONTAINS

SUBROUTINE initrandom(i,i2)
    USE varitype1
    IMPLICIT NONE
    INTEGER(k15),OPTIONAL,INTENT(IN) :: i,i2
    INTEGER(k15) :: seed_in,kl,ij
    CHARACTER(LEN=10) :: fred
    REAL(DP) :: klr

    IF (PRESENT(i)) THEN
        seed_in = i
    ELSE
        seed_in = -1
    END IF
    IF (seed_in /= -1) THEN
        IF (PRESENT(i2)) THEN
            kl = i2
            IF (i2 > 30081) STOP 'i2 > 30081'
        ELSE
            kl = 9373
        END IF
            ij = i
    ELSE
        CALL system_clock(count=ij)
        ij = MOD(ij+rand_inst*100,31328)
        CALL date_and_time(time = fred)
        READ (fred,'(e10.3)') klr
        kl = MOD(INT(klr*1000),30081)
    END IF
        CALL rmarin(ij,kl)
END SUBROUTINE initrandom

REAL(DP) FUNCTION random_real(a,b)
    USE varitype1
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: a,b

    random_real = (b-a)*ranmar() + a

END FUNCTION random_real

INTEGER(I4B) FUNCTION random_integer(o1,o2)
    USE varitype1
    IMPLICIT NONE
    INTEGER(I4B), INTENT(IN) :: o1,o2

    random_integer = INT(o2*ranmar()) + o1

END FUNCTION random_integer

SUBROUTINE rmarin(ij, kl)
    !                           0 <= IJ <= 31328
    !                           0 <= KL <= 30081
    USE varitype1
    USE raset
    IMPLICIT NONE
    INTEGER(k15), INTENT(IN) :: ij, kl
    INTEGER(K15) :: i, j, k, l, ii, jj, m
    REAL(DP) :: s, t

    IF( ij < 0  .OR.  ij > 31328  .OR. kl < 0  .OR.  kl > 30081 ) THEN
      WRITE(*, '(A)') ' The first random number seed must have a value ',  &
                   'between 0 AND 31328'
      WRITE(*, '(A)') ' The second seed must have a value between 0 and 30081'
      STOP
    END IF

    i = MOD(ij/177, 177) + 2
    j = MOD(ij, 177) + 2
    k = MOD(kl/169, 178) + 1
    l = MOD(kl, 169)
    DO ii = 1, 97
      s = 0.0D0
      t = 0.5D0
      DO jj = 1, 24
        m = MOD(MOD(i*j, 179)*k, 179)
        i = j
        j = k
        k = m
        l = MOD(53*l + 1, 169)
        IF (MOD(l*m, 64) >= 32) THEN
          s = s + t
        END IF
        t = 0.5D0*t
      END DO
      u(ii) = s
    END DO
    cc = 362436.0D0/16777216.0D0
    cd = 7654321.0D0/16777216.0D0
    cm = 16777213.0D0/16777216.0D0
    i97 = 97
    j97 = 33

    RETURN
END SUBROUTINE rmarin

FUNCTION ranmar() RESULT(fn_val)
    USE varitype1
    USE raset
    IMPLICIT NONE
    REAL(DP) :: fn_val
    ! Local variable
    REAL(DP) :: uni

    uni = u(i97) - u(j97)
    IF( uni < 0.0D0 ) uni = uni + 1.0D0
    u(i97) = uni
    i97 = i97 - 1
    IF(i97 == 0) i97 = 97
    j97 = j97 - 1
    IF(j97 == 0) j97 = 97
    cc = cc - cd
    IF( cc < 0.0D0 ) cc = cc + cm
    uni = uni - cc
    IF( uni < 0.0D0 ) uni = uni + 1.0D0
    !IF( uni == 0.0D0 ) uni = 2.0D-38
    fn_val = uni

  RETURN
END FUNCTION ranmar

END MODULE random_utilities