CODEWITHSUNDEEP
randomfortran
castled-rook
castled-rook

Reputation: 13

Random number generator produces same sequence even though it's seeded

I'm relatively new to Fortran and am trying to understand the RANDOM_NUMBER and RANDOM_SEED subroutines. The following code continually produces the same sequence of random numbers despite the fact that I seed the generator outside the DO LOOP at the top of the program.

  1 PROGRAM TEST
  2 
  3         IMPLICIT NONE
  4         
  5         INTEGER :: I, OUTPUT
  6         REAL :: R
  7 
  8         CALL RANDOM_SEED()
  9 
 10         DO I=1, 10
 11                 CALL RANDOM_NUMBER(R)
 12                 OUTPUT = I*R
 13                 PRINT *,'Random number ', I, ' = ',  OUTPUT
 14         END DO
 15 
 16 END PROGRAM TEST

Here is the output when I run the code

 Random number            1  =            0
 Random number            2  =            1
 Random number            3  =            2
 Random number            4  =            2
 Random number            5  =            1
 Random number            6  =            2
 Random number            7  =            0
 Random number            8  =            0
 Random number            9  =            3
 Random number           10  =            3

I get this exact sequence every time I run the code. I've even tried recompiling to see if the generator would reseed at compile time.

Upvotes: 1

Views: 1004

Answers (1)

Jodrell
Jodrell

Reputation: 35716

I'm no expert but looking at some documentation, I think you need to call RANDOM_SEED with a number, otherwise it initializes the random number generator to a default state.

The "default state" is implementation dependent and will be a fixed value on some platforms.

If you need to write portable code you should do something like this,

 subroutine init_random_seed()
    use iso_fortran_env, only: int64
    implicit none
    integer, allocatable :: seed(:)
    integer :: i, n, un, istat, dt(8), pid
    integer(int64) :: t

    call random_seed(size = n)
    allocate(seed(n))
    ! First try if the OS provides a random number generator
    open(newunit=un, file="/dev/urandom", access="stream", &
         form="unformatted", action="read", status="old", iostat=istat)
    if (istat == 0) then
       read(un) seed
       close(un)
    else
       ! Fallback to XOR:ing the current time and pid. The PID is
       ! useful in case one launches multiple instances of the same
       ! program in parallel.
       call system_clock(t)
       if (t == 0) then
          call date_and_time(values=dt)
          t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
               + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
               + dt(3) * 24_int64 * 60 * 60 * 1000 &
               + dt(5) * 60 * 60 * 1000 &
               + dt(6) * 60 * 1000 + dt(7) * 1000 &
               + dt(8)
       end if
       pid = getpid()
       t = ieor(t, int(pid, kind(t)))
       do i = 1, n
          seed(i) = lcg(t)
       end do
    end if
    call random_seed(put=seed)
  contains
    ! This simple PRNG might not be good enough for real work, but is
    ! sufficient for seeding a better PRNG.
    function lcg(s)
      integer :: lcg
      integer(int64) :: s
      if (s == 0) then
         s = 104729
      else
         s = mod(s, 4294967296_int64)
      end if
      s = mod(s * 279470273_int64, 4294967291_int64)
      lcg = int(mod(s, int(huge(0), int64)), kind(0))
    end function lcg
  end subroutine init_random_seed

Upvotes: 1

Related Questions