quickrnd.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2024 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade, A. Buccheri, H. Menke
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program 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
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module quickrnd_oct_m
  use, intrinsic :: iso_fortran_env
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private ::                    &
    xorshift64,                 &
    to_double
  public ::                     &
    quickrnd,                   &
    shiftseed,                  &
    reservoir_sampling_aexpj,   &
    fisher_yates_shuffle
 
  ! xorshift64 needs a good seed, in particular it cannot be seeded with 0.
  ! These seeds are generated with the splitmix64 generator and have a
  ! reasonable amount of entropy.
  integer(int64), parameter, public ::                &
    splitmix64_123 = int(z'B4DC9BD462DE412B', int64), & ! splitmix64(123)
    splitmix64_321 = int(z'E95F1C4EFCF85DEE', int64)    ! splitmix64(321)
 
  interface quickrnd
    module procedure dquickrnd_single, dquickrnd_array, zquickrnd_array
  end interface quickrnd
 
contains
 
  function xorshift64(x) result(res)
    integer(int64), intent(inout) :: x
    integer(int64) :: res
 
    assert(x /= 0_int64)
 
    x = ieor(x, shiftl(x, 13))
    x = ieor(x, shiftr(x, 7))
    x = ieor(x, shiftl(x, 17))
 
    res = x
  end function xorshift64
 
  function to_double(x) result(d)
    integer(int64), intent(in) :: x
    real(real64) :: d
    real(real64), parameter :: two_pow_minus_53 = real(z'3CA0000000000000', real64)
 
    d = real(shiftr(x, 11), real64) * two_pow_minus_53
  end function to_double
 
  ! ---------------------------------------------------------
 
  subroutine dquickrnd_single(iseed, rnd)
    integer(int64), intent(inout) :: iseed
    real(real64), intent(out) :: rnd
 
    ! no PUSH_SUB, called too often
 
    rnd = to_double(xorshift64(iseed))
 
  end subroutine dquickrnd_single
 
  ! ---------------------------------------------------------
 
  subroutine dquickrnd_array(iseed, nn, rnd)
    integer(int64), intent(inout) :: iseed
    integer, intent(in) :: nn
    real(real64), intent(inout) :: rnd(:)
 
    integer :: ii
 
    push_sub(quickrnd_array)
 
    do ii = 1, nn
      rnd(ii) = to_double(xorshift64(iseed))
    end do
 
    pop_sub(quickrnd_array)
 
  end subroutine dquickrnd_array
 
  ! ---------------------------------------------------------
 
  subroutine zquickrnd_array(iseed, nn, rnd)
    integer(int64), intent(inout) :: iseed
    integer, intent(in) :: nn
    complex(real64), intent(inout) :: rnd(:)
 
    real(real64), parameter :: inv_sqrt2 = 1.0_real64 / sqrt(m_two)
    real(real64) :: re, im
    integer :: ii
 
    push_sub(quickrnd_array)
 
    do ii = 1, nn
      re = to_double(xorshift64(iseed))
      im = to_double(xorshift64(iseed))
      rnd(ii) = cmplx(re, im, real64) * inv_sqrt2
    end do
 
    pop_sub(quickrnd_array)
 
  end subroutine zquickrnd_array
 
  ! ---------------------------------------------------------
 
  subroutine shiftseed(iseed, n)
    integer(int64), intent(inout) :: iseed
    integer(int64), intent(in) :: n
 
    integer(int64) :: discard
    integer(int64) :: ii
 
    push_sub(shiftseed)
 
    do ii = 1, n
      discard = xorshift64(iseed)
    end do
 
    pop_sub(shiftseed)
 
  end subroutine shiftseed
 
 
  ! TODO(Alex) Issue #1074 . When MR !2523 is merged, scrap `xorshift_rnd_32` and `random_real_in_range`
  ! in favour of functions introduced. Use them to implement `random_number_in_range`, as described in the
  ! issue.
  subroutine xorshift_rnd_32_biased(x)
    integer(int32), intent(inout) :: x
 
    x = ieor(x, ishft(x,  13_int32))  ! x ^= x << 13
    x = ieor(x, ishft(x, -17_int32))  ! x ^= x >> 17
    x = ieor(x, ishft(x,   5_int32))  ! x ^= x << 5
 
    ! Ensure the highest bit is zero.
    ! This keeps the result in the range [0, 2^31 - 1] (positive) BUT biases against zero
    x = iand(x, z'7FFFFFFF')
 
  end subroutine xorshift_rnd_32_biased
 
 
  function random_real_in_range(x_min, x_max, seed) result(rand)
    real(real64),   intent(in   )  :: x_min, x_max
    integer(int32), intent(inout)  :: seed
    real(real64)    :: rand
 
    real(real64), parameter :: max_value_xorshift_32 = 2147483647_real64 
    real(real64) :: scale
 
    call xorshift_rnd_32_biased(seed)
    scale = real(seed, real64) / max_value_xorshift_32
    rand = scale * (x_max - x_min) + x_min
 
  end function random_real_in_range
 
 
  ! TODO Alex. Issue #1074 Review the use of this function
  function random_integer_in_range(x_min, x_max, seed) result(rand)
    integer(int64), intent(in   ) :: x_min, x_max
    integer(int64), intent(inout) :: seed
    integer(int64)                :: rand
 
    integer(int64) :: range
 
    seed = xorshift64(seed)
    ! Is this is required?
    if (seed < 0_int64) seed = -seed  ! Enforce a positive result
    range = x_max - x_min + 1_int64
    rand = x_min + mod(seed, range)
 
  end function random_integer_in_range
 
 
  subroutine fisher_yates_shuffle(m, n, seed, values)
    integer(int32), intent(in   ) :: m
    integer(int64), intent(in   ) :: n
    integer(int64), intent(inout) :: seed
    integer(int64), intent(  out) :: values(:)
 
    integer(int64), allocatable :: indices(:)
    integer(int64) :: i, j, temp
 
    push_sub(fisher_yates_shuffle)
 
    safe_allocate(indices(1:n))
    do j = 1, n
      indices(j) = j
    end do
 
    do i = 1, m
      j = random_integer_in_range(i, n, seed)
      ! Perform the swap
      temp = indices(i)
      indices(i) = indices(j)
      indices(j) = temp
      values(i) = indices(i)
    end do
 
    values(m) = indices(m)
    safe_deallocate_a(indices)
 
    pop_sub(fisher_yates_shuffle)
 
  end subroutine fisher_yates_shuffle
 
 
  subroutine reservoir_sampling_aexpj(weight, reservoir, seed_value)
    real(real64),   intent(in )           :: weight(:)
    integer(int32), intent(out)           :: reservoir(:)
    integer(int32), intent(in ), optional :: seed_value
 
    integer(int32) :: m
    integer(int32) :: n
    real(real64)   :: u_i, u_i2
    real(real64)   :: xw
    real(real64), allocatable :: key(:)
 
    integer(int32)              :: min_key_index, i, seed_size
    integer(int32), allocatable :: seed(:)
    real(real64)                :: thres_w
 
    ! Random number range must be (0, 1), as log(0) would result in a key of -inf, and log(1) = 0
    ! Set range to [a, b], chosen empirically (may not be optimal)
    real(real64), parameter :: a = 1.e-10_real64
    real(real64), parameter :: b = 1._real64 - 1.e-10_real64
 
    push_sub(reservoir_sampling_aexpj)
 
    m = size(reservoir)
    n = size(weight)
    ! Cannot ask for more numbers than are present in the population
    assert(m <= n)
 
    if (present(seed_value)) then
      safe_allocate_source(seed(1:1), seed_value)
    else
      ! Fortran intrinsic returns an array of seeds
      call random_seed(size=seed_size)
      safe_allocate(seed(1:seed_size))
      call random_seed(get=seed)
    end if
 
    ! Initialise the reservoir with the first m items of the population
    min_key_index = 1
    safe_allocate(key(1:m))
 
    do i = 1, m
      reservoir(i) = i
      u_i = random_real_in_range(a, b, seed(1))
      key(i) = u_i ** (1._real64 / weight(i))
      if (key(i) < key(min_key_index)) min_key_index = i
    enddo
 
    u_i = random_real_in_range(a, b, seed(1))
    xw = log(u_i) / log(key(min_key_index))
 
    ! Perform swaps at jump-intervals, until the population is exhausted
    i = m + 1
    do while(i <= n)
      xw = xw - weight(i)
      if (xw <= 0._real64) then
        reservoir(min_key_index) = i
        thres_w = key(min_key_index)**weight(i)
        ! U__{i2} \in (t_w, 1)
        u_i2 = random_real_in_range(thres_w + a, b, seed(1))
        key(min_key_index) = u_i2 ** (1._real64 / weight(i))
        min_key_index = minloc(key, dim=1)
        u_i = random_real_in_range(a, b, seed(1))
        xw = log(u_i) / log(key(min_key_index))
      endif
      i = i + 1
    enddo
 
    safe_deallocate_a(key)
    safe_deallocate_a(seed)
 
    pop_sub(reservoir_sampling_aexpj)
 
  end subroutine reservoir_sampling_aexpj
 
end module quickrnd_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: