!----------------------------------------------------------------------------------------
!
! This file is part of CosmicFish.
!
! Copyright (C) 2015-2017 by the CosmicFish authors
!
! The CosmicFish code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file LICENSE at
! the top level of the CosmicFish distribution.
!
!----------------------------------------------------------------------------------------
!> @file 002_Random_utilities.f90
!! This file contains several subroutine and functions that can be used to
!! generate random numbers from different distributions.
!----------------------------------------------------------------------------------------
!> This module contains the subroutine and functions to generate random numbers from
!! different distributions.
!! This file contains part of code from the CosmoMC code with some additions and some modifications.
!! The part coming from the CosmoMC code is copyrighted by its author.
!> @author Antony Lewis (CosmoMC part), Marco Raveri
module Random_utilities
use precision
implicit none
integer :: rand_inst = 0
integer, parameter :: krand = KIND(1.d0)
integer :: Rand_Feedback = 1
contains
! ---------------------------------------------------------------------------------------------
subroutine initRandom(i, i2)
implicit none
integer, optional, intent(in) :: i
integer, optional, intent(in) :: i2
integer seed_in,kl,ij
character(len=10) :: fred
real(krand) :: klr
if (present(i)) then
seed_in = i
else
seed_in = -1
end if
if (seed_in /=-1) then
if (present(i2)) then
kl=i2
if (i2 > 30081) stop 'initRandom:second seed too large'
else
kl = 9373
end if
ij = i
else
call system_clock(count=ij)
ij = mod(ij + rand_inst*100, 31328)
call date_and_time(time=fred)
read (fred,'(e10.3)') klr
kl = mod(int(klr*1000), 30081)
end if
if (Rand_Feedback > 0 ) write(*,'(" Random seeds:",1I6,",",1I6," rand_inst:",1I4)') ij,kl,rand_inst
call rmarin(ij,kl)
end subroutine initRandom
! ---------------------------------------------------------------------------------------------
!> This function returns a random number generated from the uniform distribution in [a,b(.
!! To improve performances no check is done in order to ensure that the random number generator
!! is initialised.
real(dl) function random_uniform(a,b)
implicit none
real(dl), intent(in) :: a !< lower bound of the interval
real(dl), intent(in) :: b !< upper bound on the interval
random_uniform = (b-a)*RANMAR() + a
end function random_uniform
! ---------------------------------------------------------------------------------------------
!> This function returns a random number generated from the gaussian distribution with mean zero
!! and variance one. It works with the Box-Muller transformation.
!! To improve performances no check is done in order to ensure that the random number generator
!! is initialised.
real(dl) function gaussian_1()
implicit none
real(dl) :: R, V1, V2, FAC
integer , save :: iset = 0
real(dl), save :: gset
if ( iset==0 ) then
R = 2._dl
do while (R >= 1._dl)
V1 = 2.d0*ranmar()-1.d0
V2 = 2.d0*ranmar()-1.d0
R = V1**2+V2**2
end do
FAC = sqrt(-2.d0*log(R)/R)
GSET = V1*FAC
gaussian_1 = V2*FAC
iset = 1
else
gaussian_1 = GSET
iset = 0
endif
end function gaussian_1
! ---------------------------------------------------------------------------------------------
!> This function returns a random number generated from the gaussian distribution with mean mu
!! and variance sigma.
!! To improve performances no check is done in order to ensure that the random number generator
!! is initialised.
real(dl) function random_gaussian( mu, sigma )
implicit none
real(dl), intent(in) :: mu !< the mean of the distribution
real(dl), intent(in) :: sigma !< the variance of the distribution
random_gaussian = gaussian_1()*sigma + mu
end function random_gaussian
! This random number generator originally appeared in ''Toward a Universal
! Random Number Generator'' by George Marsaglia and Arif Zaman.
! Florida State University Report: FSU-SCRI-87-50 (1987)
!
! It was later modified by F. James and published in ''A Review of Pseudo-
! random Number Generators''
!
! THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
! (However, a newly discovered technique can yield
! a period of 10^600. But that is still in the development stage.)
!
! It passes ALL of the tests for random number generators and has a period
! of 2^144, is completely portable (gives bit identical results on all
! machines with at least 24-bit mantissas in the floating point
! representation).
!
! The algorithm is a combination of a Fibonacci sequence (with lags of 97
! and 33, and operation "subtraction plus one, modulo one") and an
! "arithmetic sequence" (using subtraction).
!
! On a Vax 11/780, this random number generator can produce a number in
! 13 microseconds.
!========================================================================
!
! PROGRAM TstRAN
! INTEGER IJ, KL, I
! Thee are the seeds needed to produce the test case results
! IJ = 1802
! KL = 9373
!
!
! Do the initialization
! call rmarin(ij,kl)
!
! Generate 20000 random numbers
! do 10 I = 1, 20000
! x = RANMAR()
!10 continue
!
! If the random number generator is working properly, the next six random
! numbers should be:
! 6533892.0 14220222.0 7275067.0
! 6172232.0 8354498.0 10633180.0
!
!
!
! write(6,20) (4096.0*4096.0*RANMAR(), I=1,6)
!20 format (3f12.1)
! end
!
subroutine RMARIN(IJ,KL)
! This is the initialization routine for the random number generator RANMAR()
! NOTE: The seed variables can have values between: 0 <= IJ <= 31328
! 0 <= KL <= 30081
!The random number sequences created by these two seeds are of sufficient
! length to complete an entire calculation with. For example, if sveral
! different groups are working on different parts of the same calculation,
! each group could be assigned its own IJ seed. This would leave each group
! with 30000 choices for the second seed. That is to say, this random
! number generator can create 900 million different subsequences -- with
! each subsequence having a length of approximately 10^30.
!
! Use IJ = 1802 & KL = 9373 to test the random number generator. The
! subroutine RANMAR should be used to generate 20000 random numbers.
! Then display the next six random numbers generated multiplied by 4096*4096
! If the random number generator is working properly, the random numbers
! should be:
! 6533892.0 14220222.0 7275067.0
! 6172232.0 8354498.0 10633180.0
double precision U(97), C, CD, CM, S, T
integer I97, J97,i,j,k,l,m
integer ij,kl
integer ii,jj
! INTEGER IRM(103)
common /RASET1/ U, C, CD, CM, I97, J97
if( IJ < 0 .or. IJ > 31328 .or. &
KL < 0 .or. KL > 30081 ) then
print '(A)', ' The first random number seed must have a value between 0 and 31328'
print '(A)',' The second seed must have a value between 0 and 30081'
stop
endif
I = mod(IJ/177, 177) + 2
J = mod(IJ , 177) + 2
K = mod(KL/169, 178) + 1
L = mod(KL, 169)
do II = 1, 97
S = 0.0
T = 0.5
do JJ = 1, 24
M = mod(mod(I*J, 179)*K, 179)
I = J
J = K
K = M
L = mod(53*L+1, 169)
if (mod(L*M, 64) >= 32) then
S = S + T
endif
T = 0.5 * T
3 continue
end do
U(II) = S
2 continue
end do
C = 362436.0 / 16777216.0
CD = 7654321.0 / 16777216.0
CM = 16777213.0 /16777216.0
I97 = 97
J97 = 33
end subroutine RMARIN
! ---------------------------------------------------------------------------------------------
!> This function returns a random number generated from the uniform distribution in [0,1(
!! This is the random number generator proposed by George Marsaglia in
!! Florida State University Report: FSU-SCRI-87-50
!! It was slightly modified by F. James to produce an array of pseudorandom
!! numbers.
double precision function RANMAR()
double precision U(97), C, CD, CM
integer I97, J97
double precision uni
common /RASET1/ U, C, CD, CM, I97, J97
! INTEGER IVEC
UNI = U(I97) - U(J97)
if( UNI < 0.0 ) UNI = UNI + 1.0
U(I97) = UNI
I97 = I97 - 1
if(I97 == 0) I97 = 97
J97 = J97 - 1
if(J97 == 0) J97 = 97
C = C - CD
if( C < 0.d0 ) C = C + CM
UNI = UNI - C
if( UNI < 0.d0 ) UNI = UNI + 1.0
RANMAR = UNI
end function RANMAR
! ---------------------------------------------------------------------------------------------
end module Random_utilities
!----------------------------------------------------------------------------------------