hammersley Subroutine

  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