mercredi 1 juillet 2015

Slow random seed generator--why?

I've been playing around with various random seed generators. Here is a simple one:

subroutine init_random_seed()
  integer :: i, n, clock
  integer, dimension(:), allocatable :: seed
  call random_seed(size = n)
  allocate(seed(n))
  call system_clock(count=clock)
  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
  call random_seed(put = seed)
  deallocate(seed)
end

...and a more robust one:

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))
  OPEN(newunit=un, file='/dev/urandom', access='stream', status='old', action='read', form='unformatted', iostat=istat)
  IF (istat == 0) THEN
    READ(un) seed
    CLOSE(un)
  ELSE
    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)
  DEALLOCATE(seed)
CONTAINS
  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

The second one generates higher-quality random numbers, but is comparatively slow. Does anyone see why?




Aucun commentaire:

Enregistrer un commentaire