program halton_prb
!
!*******************************************************************************
!
!! HALTON_PRB calls a set of problems for HALTON.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HALTON_PRB'
  write ( *, '(a)' ) '  A set of tests for HALTON.'
 
  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10

  call test11
  call test12

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'HALTON_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests HALTON_NUMBER.
!
  implicit none
!
  integer i
  real r
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  HALTON_NUMBER computes the elements of a Halton'
  write ( *, '(a)' ) '  sequence.'
  write ( *, '(a)' ) '  Each call produces the next value.  By default,'
  write ( *, '(a)' ) '  the base is 2, and the sequence starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we call HALTON_NUMBER several times.'

  seed = 1
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    call halton_number ( r )
    write ( *, '(i6,g14.6)' ) i, r
  end do

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests HALTON_NUMBER_SEQUENCE.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer i
  real r(n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  HALTON_NUMBER_SEQUENCE computes several elements of '
  write ( *, '(a)' ) '  a Halton sequence on a single call.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  By default, the base is 2, and the sequence starts '
  write ( *, '(a)' ) '  at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we call HALTON_NUMBER_SEQUENCE once.'

  seed = 1
  call halton_seed_set ( seed )

  call halton_number_sequence ( n, r )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i6,g14.6)' ) i, r(i)
  end do

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests HALTON_NUMBER_SEQUENCE.
!! TEST03 tests HALTON_SEED_SET.
!! TEST03 tests HALTON_SEED_GET.
!
  implicit none
!
  integer, parameter :: nmax = 10
!
  integer i
  integer n
  real r(nmax)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  HALTON_SEED_SET specifies the next element of the'
  write ( *, '(a)' ) '    Halton sequence to compute.'
  write ( *, '(a)' ) '  HALTON_SEED_GET reports the next element of the'
  write ( *, '(a)' ) '    Halton sequence that will be computed.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  By default, the sequence starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we demonstrate computing elements'
  write ( *, '(a)' ) '  affects the seed, and how resetting the seed determines'
  write ( *, '(a)' ) '  the next element computed.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We start at element 1 and compute 10 elements.'
  write ( *, '(a)' ) ' '

  seed = 1
  call halton_seed_set ( seed )

  n = 10
  call halton_number_sequence ( n, r )

  do i = 1, n
    write ( *, '(i6,g14.6)' ) i, r(i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We jump back to element 6 and compute 10 elements.'
  write ( *, '(a)' ) ' '

  seed = 6
  call halton_seed_set ( seed )

  n = 10
  call halton_number_sequence ( n, r )

  do i = 1, n
    write ( *, '(i6,g14.6)' ) i + seed - 1, r(i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We restart at element 0 and compute 6 elements.'
  write ( *, '(a)' ) ' '

  seed = 0
  call halton_seed_set ( seed )

  n = 6
  call halton_number_sequence ( n, r )

  do i = 1, n
    write ( *, '(i6,g14.6)' ) i + seed - 1, r(i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We jump to element 100 and compute 5 elements.'
  write ( *, '(a)' ) ' '

  seed = 100
  call halton_seed_set ( seed )

  n = 5
  call halton_number_sequence ( n, r )

  do i = 1, n
    write ( *, '(i6,g14.6)' ) i + seed - 1, r(i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests HALTON_NUMBER.
!! TEST04 tests HALTON_NUMBER_BASE_GET.
!! TEST04 tests HALTON_NUMBER_BASE_SET.
!
  implicit none
!
  integer base
  integer i
  real r
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  HALTON_NUMBER_BASE_GET gets the current Halton base.'
  write ( *, '(a)' ) '  HALTON_NUMBER_BASE_SET sets the current Halton base.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The Halton base should be a prime, but this is'
  write ( *, '(a)' ) '  not checked.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we compute Halton numbers with the'
  write ( *, '(a)' ) '  default base, then change the base, reset the seed,'
  write ( *, '(a)' ) '  and recompute the numbers.'

  seed = 1
  call halton_seed_set ( seed )

  call halton_number_base_get ( base )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  HALTON_NUMBER_BASE_GET: Current base is ', base
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    call halton_number ( r )
    write ( *, '(i6,g14.6)' ) i, r
  end do

  base = 3
  call halton_number_base_set ( base )
  seed = 1
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Reset base to ', base
  write ( *, '(a,i12)' ) '  Reset seed to ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    call halton_number ( r )
    write ( *, '(i6,g14.6)' ) i, r
  end do

  base = 4
  call halton_number_base_set ( base )
  seed = 1
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Reset base to ', base
  write ( *, '(a,i12)' ) '  Reset seed to ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    call halton_number ( r )
    write ( *, '(i6,g14.6)' ) i, r
  end do

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests HALTON_VECTOR.
!
  implicit none
!
  integer, parameter :: ndim = 4
!
  integer i
  real r(ndim)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  HALTON_VECTOR computes the elements of a vector '
  write ( *, '(a)' ) '  Halton sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Each call produces the next value.  By default,'
  write ( *, '(a)' ) '  the bases are the first NDIM primes, and the sequence '
  write ( *, '(a)' ) '  starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we call HALTON_VECTOR several times,'
  write ( *, '(a)' ) '  with the default bases.'

  seed = 1
  call halton_seed_set ( seed )

  call halton_vector_ndim_set ( ndim )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    call halton_vector ( ndim, r )
    write ( *, '(i6,4g14.6)' ) i, r(1:ndim)
  end do

  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests HALTON_VECTOR_SEQUENCE.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: ndim = 4
!
  integer i
  real r(ndim,n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  HALTON_VECTOR_SEQUENCE computes the next N elements' 
  write ( *, '(a)' ) '  of a vector Halton sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Each call produces the next value.  By default,'
  write ( *, '(a)' ) '  the bases are the first NDIM primes, and the sequence '
  write ( *, '(a)' ) '  starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we demonstrate how one call can compute'
  write ( *, '(a)' ) '  many successive vector elements of the sequence.'

  seed = 1
  call halton_seed_set ( seed )

  call halton_vector_ndim_set ( ndim )

  call halton_vector_sequence ( ndim, n, r )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(i6,4g14.6)' ) i, r(1:ndim,i)
  end do

  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests HALTON_SEED_GET.
!! TEST07 tests HALTON_SEED_SET.
!! TEST07 tests HALTON_VECTOR_SEQUENCE.
!
  implicit none
!
  integer, parameter :: ndim = 4
  integer, parameter :: nmax = 10
!
  integer i
  integer n
  real r(ndim,nmax)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  HALTON_SEED_SET specifies the next element of the'
  write ( *, '(a)' ) '    Halton sequence to compute.'
  write ( *, '(a)' ) '  HALTON_SEED_GET reports the next element of the'
  write ( *, '(a)' ) '    Halton sequence that will be computed.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  By default, the sequence starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we show how the seed changes as we'
  write ( *, '(a)' ) '  compute elements, and how changing the seed changes'
  write ( *, '(a)' ) '  the next element computed.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We start at element 1 and compute 10 elements.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  We use a vector Halton sequence, with NDIM = ', ndim

  seed = 1
  call halton_seed_set ( seed )

  call halton_vector_ndim_set ( ndim )

  n = 10
  call halton_vector_sequence ( ndim, n, r )

  do i = 1, n
    write ( *, '(i6,4g14.6)' ) i, r(1:ndim,i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We jump back to element 6 and compute 10 elements.'
  write ( *, '(a)' ) ' '

  seed = 6
  call halton_seed_set ( seed )

  n = 10
  call halton_vector_sequence ( ndim, n, r )

  do i = 1, n
    write ( *, '(i6,4g14.6)' ) i + seed - 1, r(1:ndim,i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We restart at element 0 and compute 6 elements.'
  write ( *, '(a)' ) ' '

  seed = 0
  call halton_seed_set ( seed )

  n = 6
  call halton_vector_sequence ( ndim, n, r )

  do i = 1, n
    write ( *, '(i6,4g14.6)' ) i + seed - 1, r(1:ndim,i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We jump to element 100 and compute 5 elements.'
  write ( *, '(a)' ) ' '

  seed = 100
  call halton_seed_set ( seed )

  n = 5
  call halton_vector_sequence ( ndim, n, r )

  do i = 1, n
    write ( *, '(i6,4g14.6)' ) i + seed - 1, r(1:ndim,i)
  end do

  call halton_seed_get ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i12)' ) '  HALTON_SEED_GET reports current seed is ', seed

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests HALTON_VECTOR_BASE_GET.
!! TEST08 tests HALTON_VECTOR_BASE_SET.
!! TEST08 tests HALTON_VECTOR_SEQUENCE.
!
  implicit none
!
  integer, parameter :: n = 10
  integer, parameter :: ndim = 4
!
  integer base(ndim)
  integer i
  integer prime
  real r(ndim,n)
  integer seed
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  HALTON_VECTOR_BASE_GET gets the current bases.'
  write ( *, '(a)' ) '  HALTON_VECTOR_BASE_SET sets the current bases.'
  write ( *, '(a)' ) '  HALTON_VECTOR_SEQUENCE computes the next N elements'
  write ( *, '(a)' ) '  of a vector Halton sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Each call produces the next value.  By default,'
  write ( *, '(a)' ) '  the bases are the first NDIM primes, and the sequence '
  write ( *, '(a)' ) '  starts at element 1.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we compute the first 10 elements of the'
  write ( *, '(a)' ) '  default sequence, then change bases, reset the seed'
  write ( *, '(a)' ) '  to 1, and recompute the first 10 elements.'

  seed = 1
  call halton_seed_set ( seed )

  call halton_vector_ndim_set ( ndim )

  call halton_vector_sequence ( ndim, n, r )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(i6,4g14.6)' ) i, r(1:ndim,i)
  end do

  call halton_vector_base_get ( ndim, base )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Current Halton bases:'
  write ( *, '(a)' ) ' '
  write ( *, '(5i6)' ) base(1:ndim)

  do i = 1, ndim
    base(i) = prime(2*i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  New Halton bases:'
  write ( *, '(a)' ) ' '
  write ( *, '(5i6)' ) base(1:ndim)

  call halton_vector_base_set ( ndim, base )

  seed = 1
  call halton_seed_set ( seed )

  call halton_vector_sequence ( ndim, n, r )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Element   Halton'
  write ( *, '(a)' ) ' '
  do i = 1, 10
    write ( *, '(i6,4g14.6)' ) i, r(1:ndim,i)
  end do

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests SPHERE_UNIT_HALTON_2D.
!
  integer, parameter :: n = 2
!
  real average(n)
  real average_dot
  real dot_average
  integer i
  integer j
  integer, parameter :: n_sample = 1000
  integer seed
  real v(n)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  For the unit sphere in 2 dimensions (the circle):'
  write ( *, '(a)' ) '  SPHERE_UNIT_HALTON_2D samples;'

  seed = 0
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A few sample values:'
  write ( *, '(a)' ) ' '

  do i = 1, 5
    call sphere_unit_halton_2d ( x )
    write ( *, '(2f8.4)' ) x(1:n)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of sample points = ', n_sample

  seed = 0
  call halton_seed_set ( seed )

  average(1:n) = 0.0E+00

  do i = 1, n_sample
    call sphere_unit_halton_2d ( x )
    average(1:n) = average(1:n) + x(1:n)
  end do

  average(1:n) = average(1:n) / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the points, which should get a value'
  write ( *, '(a)' ) '  close to zero, and closer as N increases.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2f8.4)' ) '  Average:        ', average(1:n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now choose a random direction, sample the same'
  write ( *, '(a)' ) '  number of points, and compute the dot product with'
  write ( *, '(a)' ) '  the direction.'
  write ( *, '(a)' ) '  Take the absolute value of each dot product '
  write ( *, '(a)' ) '  and sum and average.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We expect a value near 2 / PI = 0.6366...'

  do j = 1, 5

    call halton_seed_randomize ( seed )

    call sphere_unit_halton_2d ( v )

    seed = 0
    call halton_seed_set ( seed )

    dot_average = 0.0E+00

    do i = 1, n_sample
      call sphere_unit_halton_2d ( x )
      dot_average = dot_average + abs ( dot_product ( x(1:n), v(1:n) ) )
    end do

    dot_average = dot_average / real ( n_sample )

    write ( *, '(a)' ) ' '
    write ( *, '(a,2f8.4)' ) '  V:                ', v(1:n)
    write ( *, '(a, f8.4)' ) '  Average |(XdotV)| ', dot_average

  end do

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests BALL_UNIT_HALTON_2D.
!
  integer, parameter :: n = 2
!
  real atan4
  real average(n)
  real average_r
  real average_theta
  integer i
  integer j
  integer, parameter :: n_sample = 1000
  real r_pi
  integer seed
  real theta
  real v(n)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  For the unit ball in 2 dimensions (the disk):'
  write ( *, '(a)' ) '  BALL_UNIT_HALTON_2D samples;'

  seed = 0
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A few sample values:'
  write ( *, '(a)' ) ' '

  do i = 1, 5
    call ball_unit_halton_2d ( x )
    write ( *, '(2f8.4)' ) x(1:n)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of sample points = ', n_sample

  seed = 0
  call halton_seed_set ( seed )

  average(1:n) = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_2d ( x )
    average(1:n) = average(1:n) + x(1:n)
  end do

  average(1:n) = average(1:n) / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the points, which should get a value'
  write ( *, '(a)' ) '  close to zero, and closer as N increases.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2f8.4)' ) '  Average:        ', average(1:n)

  seed = 0
  call halton_seed_set ( seed )

  average_r = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_2d ( x )
    r = sqrt ( sum ( x(1:n)**2 ) )
    average_r = average_r + r
  end do

  average_r = average_r / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the distance of the points from' 
  write ( *, '(a,f8.4)' ) '  the center, which should be 1/2**(1/n) = ', &
    0.5E+00**( 1.0E+00 / real ( n ) )
  write ( *, '(a)' ) ' '
  write ( *, '(a,f8.4)' ) '  Average:        ', average_r

  seed = 0
  call halton_seed_set ( seed )

  average_theta = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_2d ( x )
    theta = atan4 ( x(2), x(1) )
    average_theta = average_theta + theta
  end do

  average_theta = average_theta / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the angle THETA,'
  write ( *, '(a)' ) '  which should be PI.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,f8.4)' ) '  Average:        ', average_theta

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests SPHERE_UNIT_HALTON_3D.
!
  integer, parameter :: n = 3
!
  real average(n)
  real dot_average
  integer i
  integer j
  integer, parameter :: n_sample = 1000
  integer seed
  real v(n)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  For the unit sphere in 3 dimensions:'
  write ( *, '(a)' ) '  SPHERE_UNIT_HALTON_3D samples;'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A few sample values:'
  write ( *, '(a)' ) ' '

  seed = 123456789
  call random_seed ( seed )

  do i = 1, 5
    call sphere_unit_halton_3d ( x )
    write ( *, '(3f8.4)' ) x(1:n)
  end do

  average(1:n) = 0.0E+00

  do i = 1, n_sample
    call sphere_unit_halton_3d ( x )
    average(1:n) = average(1:n) + x(1:n)
  end do

  average(1:n) = average(1:n) / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the points, which should get a value'
  write ( *, '(a)' ) '  close to zero, and closer as N increases.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,3f8.4)' ) '  Average:        ', average(1:n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now choose a random direction, sample the same'
  write ( *, '(a)' ) '  number of points, and compute the dot product with'
  write ( *, '(a)' ) '  the direction.'
  write ( *, '(a)' ) '  Take the absolute value of each dot product '
  write ( *, '(a)' ) '  and sum and average.'

  do j = 1, 5

    call halton_seed_randomize ( seed )

    call sphere_unit_halton_3d ( v )

    seed = 0
    call halton_seed_set ( seed )

    dot_average = 0.0E+00

    do i = 1, n_sample
      call sphere_unit_halton_3d ( x )
      dot_average = dot_average + abs ( dot_product ( x(1:n), v(1:n) ) )
    end do

    dot_average = dot_average / real ( n_sample )

    write ( *, '(a)' ) ' '
    write ( *, '(a,3f8.4)' ) '  V:                ', v(1:n)
    write ( *, '(a, f8.4)' ) '  Average |(XdotV)| ', dot_average

  end do

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests BALL_UNIT_HALTON_3D.
!
  integer, parameter :: n = 3
!
  real atan4
  real average(n)
  real average_phi
  real average_r
  real average_theta
  integer i
  integer j
  integer, parameter :: n_sample = 1000
  real r_pi
  integer seed
  real theta
  real v(n)
  real x(n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  For the unit ball in 3 dimensions:'
  write ( *, '(a)' ) '  BALL_UNIT_HALTON_3D samples;'

  seed = 0
  call halton_seed_set ( seed )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A few sample values:'
  write ( *, '(a)' ) ' '

  do i = 1, 5
    call ball_unit_halton_3d ( x )
    write ( *, '(3f8.4)' ) x(1:n)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of sample points = ', n_sample

  seed = 0
  call halton_seed_set ( seed )

  average(1:n) = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_3d ( x )
    average(1:n) = average(1:n) + x(1:n)
  end do

  average(1:n) = average(1:n) / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the points, which should get a value'
  write ( *, '(a)' ) '  close to zero, and closer as N increases.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,3f8.4)' ) '  Average:        ', average(1:n)

  seed = 0
  call halton_seed_set ( seed )

  average_r = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_3d ( x )
    r = sqrt ( sum ( x(1:n)**2 ) )
    average_r = average_r + r
  end do

  average_r = average_r / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the distance of the points from' 
  write ( *, '(a,f8.4)' ) '  the center, which should be the '
  write ( *, '(a,f8.4)' ) '  1/2**(1/n) = ', &
    0.5E+00**( 1.0E+00 / real ( n ) )
  write ( *, '(a)' ) ' '
  write ( *, '(a,f8.4)' ) '  Average:        ', average_r

  seed = 0
  call halton_seed_set ( seed )

  average_theta = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_3d ( x )
    theta = atan4 ( x(2), x(1) )
    average_theta = average_theta + theta
  end do

  average_theta = average_theta / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the angle THETA,'
  write ( *, '(a)' ) '  which should be PI.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,f8.4)' ) '  Average:        ', average_theta

  seed = 0
  call halton_seed_set ( seed )

  average_phi = 0.0E+00

  do i = 1, n_sample
    call ball_unit_halton_3d ( x )
    r = sqrt ( sum ( x(1:n)**2 ) )
    phi = acos ( x(3) / r )
    average_phi = average_phi + phi
  end do

  average_phi = average_phi / real ( n_sample )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now average the angle PHI,'
  write ( *, '(a)' ) '  which should be PI/2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,f8.4)' ) '  Average:        ', average_phi

  return
end