random_numbers.f90

Go to the documentation of this file.

!> \file random_numbers.f90  Procedures to generate random numbers
 
 
!  Copyright (c) 2002-2025  Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
!   
!  This file is part of the libSUFR package, 
!  see: http://libsufr.sourceforge.net/
!   
!  This is free software: you can redistribute it and/or modify it under the terms of the European Union
!  Public Licence 1.2 (EUPL 1.2).  This software is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
!  PURPOSE.  See the EU Public Licence for more details.  You should have received a copy of the European
!  Union Public Licence along with this code.  If not, see <https://www.eupl.eu/1.2/en/>.
!  
!  
 
 
 
 
!***********************************************************************************************************************************
!> \brief  Procedures to generate random numbers
 
module sufr_random_numbers
  implicit none
  save
  
contains
  
  
  !*********************************************************************************************************************************
  !> \brief  Use the system clock to get a random initialisation seed (i.e., negative integer) for a random-numbed generator.
  !!
  !! \param degree         Degree of randomness: 0-completely (same result for a ms), 1-same result during a clock hour, 
  !!                       2-same result during a calendar day (int)
  !! \param date_array     Provide a date/time array as provided by date_and_time() to generate a non-random seed, e.g. for
  !!                       reproduction purposes ([year,month,day, (tz), hour,minute,second,millisecond]; optional)
  !!
  !! \retval get_ran_seed  Randon-number seed:  -1e6 < seed < 0 (int)
  
  function get_ran_seed(degree, date_array)
    implicit none
    integer, intent(in) :: degree
    integer, intent(in), optional :: date_array(8)
    integer :: get_ran_seed
    
    integer :: seed,dt(8)
    character :: tmpstr*(10)
    
    
    call date_and_time(tmpstr,tmpstr,tmpstr,dt)  ! dt: 1-year, 2-month, 3-day, 5-hour, 6-minute, 7-second, 8-millisecond
    if(present(date_array)) dt = date_array      ! Use a user-provided date/time, rather than the system clock
    
    select case (degree)
    case(1)
       seed = dt(1)*1000 + dt(2)*10000 + dt(3)*10101 + dt(5)  ! Constant result during a clock hour
    case(2)
       seed = dt(1)*1000 + dt(2)*10000 + dt(3)*10101          ! Constant result during a calendar day
    case default
       seed = dt(6)*1010 + dt(7)*10101 + dt(8)*1001           ! Different result every millisecond
    end select
    
    get_ran_seed = -abs(mod(seed+1,999999))                   ! Return a negative number, -1e6 < get_ran_seed < 0
    
  function get_ran_seed(degree, date_array) …
  end function get_ran_seed
  !*********************************************************************************************************************************
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Generate a pseudo-random number from a uniform distribution 0 < r < 1.
  !!
  !! \param seed    The seed to generate the random number from.  
  !!                Set seed<0 to initialise the generator; seed is updated between calls (int).
  !! \retval ran_unif  The random number, uniformely generated between 0 < r < 1  (double).
  !! 
  !! - Use two L'Ecuyer generators, period is \f$ \sim 10^{18}\f$
  !! - tab is a Bays-Durham shuffle table of length Ntab
  !! \see Numerical Recipes in Fortran 77, Sect.7.1.
  
  function ran_unif(seed)
    use sufr_kinds, only: double, dbl
    
    implicit none
    integer, intent(inout) :: seed
    real(double) :: ran_unif
    
    integer, parameter :: im1=2147483563, ia1=40014, iq1=53668, ir1=12211 
    integer, parameter :: im2=2147483399, ia2=40692, iq2=52774, ir2= 3791
    integer, parameter :: ntab=32, im1m1=im1-1, ndtab=67108862  ! ndtab=1+im1m1/Ntab
    
    ! rnmx should be the largest number < 1 and != 1:
    real(double), parameter :: am1  = 1.0_dbl/im1
    real(double), parameter :: eps  = epsilon(0.0_dbl)
    real(double), parameter :: rnmx = 1.0_dbl - eps
    
    integer, save :: seed2=123456789, tab(ntab+8)=0, iy=0  ! +8 needed to avoid gfortran-8 out-of-bounds warnings
    integer :: j,k
    
    if(seed.le.0) then                                  ! 'Initialise' generator
       seed = max(-seed,1)                              ! Don't allow seed=0
       seed2 = seed
       do j = ntab+8,1,-1                               ! Shuffle the table, don't save the first 8 iterations
          k = seed/iq1
          seed = ia1*(seed-k*iq1) - k*ir1
          if(seed.lt.0) seed = seed + im1
          if(j.le.ntab) tab(j) = seed
       end do
       iy = tab(1)
    end if
    
    ! Produce the random number 1:
    k = seed/iq1
    seed = ia1*(seed-k*iq1) - k*ir1                     ! Use Schrage's method to compute mod(). Update seed for next draw
    if(seed.lt.0) seed = seed + im1
    
    ! Produce the random number 2:
    k = seed2/iq2
    seed2 = ia2*(seed2-k*iq2) - k*ir2                   ! Use Schrage's method to compute mod(). Update seed for next draw
    if(seed2.lt.0) seed2 = seed2 + im2
    
    j = 1 + iy/ndtab                                    ! Result: 1 <= j <= Ntab
    iy = tab(j) - seed2                                 ! tab contains information about seed
    tab(j) = seed
    if(iy.lt.1) iy = iy + im1m1
    
    ran_unif = min(am1*iy,rnmx)                         ! Make sure r<1
    
  function ran_unif(seed) …
  end function ran_unif
  !*********************************************************************************************************************************
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Generate a pseudo-random number from a Gaussian distribution with mean 0 and standard deviation 1
  !!
  !! \param seed        The seed to generate the random number from
  !! \retval ran_gauss  The random number, using a Gaussian distribution
  !! 
  !! - Two pseudo-random numbers are drawn from a uniform distribution, using ran_unif().
  !! - These are converted into two pseudo-random numbers from a Gaussian distribution using the Box-Muller transform
  !! - Odd-numbered calls return the first, even-numbered calls the second number
  !!   - using the polar form is ~38% faster than using the basic form
  !!
  !! \see e.g.:
  !! - Numerical Recipes in Fortran 77, Sect.7.2.
  !! - http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
  
  function ran_gauss(seed)
    use sufr_kinds, only: double
    
    implicit none
    integer, intent(inout) :: seed
    real(double) :: ran_gauss, ru1,ru2, rusq, fac, rg1,rg2
    logical :: saved = .false.
    save :: rg2, saved
    
    if(.not.saved) then  ! Compute
       rusq = -1.d0
       do while(rusq.le.0.d0 .or. rusq.ge.1.d0)
          ru1 = ran_unif(seed)*2 - 1.d0
          ru2 = ran_unif(seed)*2 - 1.d0
          
          rusq = ru1**2 + ru2**2
       end do
       
       fac = sqrt(-2*log(rusq)/rusq)
       rg1 = fac * ru1
       rg2 = fac * ru2
       
       ran_gauss = rg1
       saved = .true.
    else
       ran_gauss = rg2
       saved = .false.
    end if
    
  function ran_gauss(seed) …
  end function ran_gauss
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Generate a pseudo-random number from a uniform distribution 0 < r < 1.
  !!
  !! \param  seed    The seed to generate the random number from.  
  !!                 Set seed<0 to initialise the generator; seed is updated between calls (int).
  !! \param  str     The random string  (I/O: output str has length of input str, filled with random characters)
  
  subroutine ran_str(seed, str)
    implicit none
    integer, intent(inout) :: seed
    character, intent(inout) :: str*(*)
    
    integer, parameter :: ascii_start=33, ascii_end=126  ! First and last ASCII characters to use
    integer :: dascii, ich, chari
    
    
    dascii = ascii_end - ascii_start + 1
    do ich=1,len(str)
       chari = nint(ascii_start + ran_unif(seed)*dascii)
       str(ich:ich) = char(chari)
    end do
    
  subroutine ran_str(seed, str) …
  end subroutine ran_str
  !*********************************************************************************************************************************
  
end module sufr_random_numbers
!***********************************************************************************************************************************