hammersley_inverse Subroutine

  subroutine hammersley_inverse ( r, m, n, i )

    !*****************************************************************************80
    !
    !! HAMMERSLEY_INVERSE inverts an element of the Hammersley sequence.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    20 August 2016
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real ( kind = 8 ) R(M), the I-th element of the Hammersley sequence.
    !    0 <= R < 1.0
    !
    !    Input, integer ( kind = 4 ) M, the spatial dimension.
    !
    !    Input, integer ( kind = 4 ) N, the "base" for the first component.
    !    1 <= N.
    !
    !    Output, integer ( kind = 4 ) I, the index of the element of the sequence.
    !
    implicit none

    integer ( kind = 4 ) m

    integer ( kind = 4 ) d
    integer ( kind = 4 ) i
    integer ( kind = 4 ) n
    integer ( kind = 4 ) p
    real ( kind = 8 ) r(m)
    real ( kind = 8 ) t

    if ( any ( r(1:m) < 0.0D+00 ) .or. any ( 1.0D+00 < r(1:m) ) ) then
       write ( *, '(a)' ) ''
       write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!'
       write ( *, '(a)' ) '  0 <= R <= 1.0 is required.'
       stop 1
    end if

    if ( m < 1 ) then
       write ( *, '(a)' ) ''
       write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!'
       write ( *, '(a)' ) '  1 <= M is required.'
       stop 1
    end if

    if ( n < 1 ) then
       write ( *, '(a)' ) ''
       write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!'
       write ( *, '(a)' ) '  1 <= N is required.'
       stop 1
    end if
    !
    !  Invert using the second component only, because working with base
    !  2 is accurate.
    !
    if ( 2 <= m ) then

       i = 0
       t = r(2)
       p = 1

       do while ( t /= 0.0D+00 )
          t = t * 2.0D+00
          d = int ( t )
          i = i + d * p
          p = p * 2
          t = mod ( t, 1.0D+00 )
       end do

    else

       i = nint ( real ( n, kind = 8 ) * r(1) )

    end if

    return
  end subroutine hammersley_inverse