lundi 31 août 2015

Generate *truly* random numbers fortran 90

I am relatively new to Fortran 90 and despite having scoured a lot of manuals, forums, tutorials, university pages, and books, I have not managed to modify the F90 routine below using the system clock so that the random numbers generated (in [0,1]) are independent from one run to another.

! Simple random number generator in the range [0,1]
! ranset_(iseed) initializes with iseed
! ranf_() returns next random number




      real function ranf_()
       implicit none
       real rand_
!       ranf_ = rand_(0)
        call random_number(ranf_)
       return
      end


      subroutine ranset_(iseed)
      implicit none
      real rand_,ranf_
      integer iseed, i, m, nsteps
!      i = rand_(1) ! reinitialize (reset)
      nsteps = iseed*10000
      do i = 1,nsteps
        m = ranf_()
!       m = rand_(0)
      end do
      return
      end





      real function rand_(iseed)
      implicit none
      integer iseed
      integer ia1, ia0, ia1ma0, ic, ix1, ix0, iy0, iy1
      save ia1, ia0, ia1ma0, ic, ix1, ix0
      data ix1, ix0, ia1, ia0, ia1ma0, ic/0,0,1536,1029,507,1731/
      if (iseed.ne.0) then
        ia1 = 1536
        ia0 = 1029
        ia1ma0 = 507
        ic = 1731
        ix1 = 0
        ix0 = 0
        rand_ = 0
      else
       iy0 = ia0*ix0
       iy1 = ia1*ix1 + ia1ma0*(ix0-ix1) + iy0
       iy0 = iy0 + ic
       ix0 = mod (iy0, 2048)
       iy1 = iy1 + (iy0-ix0)/2048
       ix1 = mod (iy1, 2048)
       rand_ = ix1*2048 + ix0
       rand_ = rand_ / 4194304.
      end if
      return
      end

The code is not that easy to fully understand for a beginner like me. Any insights/suggestions/solutions on how to properly do this using the system clock?




Aucun commentaire:

Enregistrer un commentaire