@brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". ************80
HAMMERSLEY computes an element of a Hammersley sequence.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer(kind=4) | :: | i | ||||
integer(kind=4) | :: | m | ||||
integer(kind=4) | :: | n | ||||
real(kind=8) | :: | r(m) |
subroutine hammersley ( i, m, n, r )
!*****************************************************************************80
!
!! HAMMERSLEY computes an element of a Hammersley sequence.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 20 August 2016
!
! Author:
!
! John Burkardt
!
! Reference:
!
! John Hammersley,
! Monte Carlo methods for solving multivariable problems,
! Proceedings of the New York Academy of Science,
! Volume 86, 1960, pages 844-874.
!
! Parameters:
!
! Input, integer ( kind = 4 ) I, the index of the element of the sequence.
! 0 <= I.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) N, the "base" for the first component.
! 1 <= N.
!
! Output, real ( kind = 8 ) R(M), the element of the sequence with index I.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) d
integer ( kind = 4 ) i
integer ( kind = 4 ) i1
integer ( kind = 4 ) j
integer ( kind = 4 ) n
! integer ( kind = 4 ) prime
real ( kind = 8 ) prime_inv(2:m)
real ( kind = 8 ) r(m)
integer ( kind = 4 ) t(2:m)
t(2:m) = abs ( i )
!
! Carry out the computation.
!
do i1 = 2, m
prime_inv(i1) = 1.0D+00 / real ( prime ( i1 - 1 ), kind = 8 )
end do
r(1) = real ( mod ( i, n + 1 ), kind = 8 ) / real ( n, kind = 8 )
r(2:m) = 0.0D+00
do while ( any ( t(2:m) /= 0 ) )
do j = 2, m
d = mod ( t(j), prime ( j - 1 ) )
r(j) = r(j) + real ( d, kind = 8 ) * prime_inv(j)
prime_inv(j) = prime_inv(j) / real ( prime ( j - 1 ), kind = 8 )
t(j) = ( t(j) / prime ( j - 1 ) )
end do
end do
return
end subroutine hammersley