!     Last change:  DOS  31 Jul 2000    7:54 pm
!                                       *** copyright 2000 ***
! *** filename halton.f95               *** John F. Monahan **
!                                       **********************
      program halton
!  routine for generating Halton sequences
      implicit none
      real x
!  declare the variable list long enough so that the number of points
!   sought does not exceed base**length
      integer, parameter :: length = 10
      integer, dimension(length) :: list
      integer, dimension(5) :: baslst
      integer i,nn,base,nnzero,k,j
!
  21  format(//i6,' points from Van de Corput sequence for base',i4)
  22  format(i6,f12.8,9x,10i3)
!                   interface block
      interface
      subroutine ncrmnt(list,base,nnzero)
      integer, dimension(:) :: list
      integer, INTENT(IN) :: base
      integer, INTENT(IN OUT) :: nnzero
      end subroutine ncrmnt
      real function expand(list,base,nnzero)
      integer, dimension(:) :: list
      integer, INTENT(IN) :: base
      integer, INTENT(IN OUT) :: nnzero
      end function expand
      end interface
!
      data baslst/ 2, 3, 5, 7, 10 /
!
      open( unit=6, file='halton.out' )
!                   how many points to get?
      nn = 10
!                   what bases?
      do k = 1,5
      base = baslst(k)
!                   write header
      write(6,21) nn,base
!                   initialize digit counter and index vector to zero
      nnzero = 0
      list = 0                ! initialize whole vector list to zero
!                   do nn of them
      do i = 1,nn
!                   get the next number
      call ncrmnt(list,base,nnzero)
!                   get the value
      x = expand(list,base,nnzero)
      write(6,22) i,x,(list(j),j=1,nnzero)
      end do  ! loop on i
!                   done with this base
      end do  ! loop on k
      stop
      end program halton
      
      subroutine ncrmnt(list,base,nnzero)
!  increments by one a number expressed as an integer in base "base"
!   with each digit residing in "list"
!   arguments
!  list     integer vector of length at least nnzero holding digit
!            of base "base" expansion -- dummy dimensioned as length 1
!  base     integer base
!  nnzero   integer giving the number of nonzero digits
!
!  if list initialized to all zeros, then after n calls to nrcrmnt
!                                 nnzero
!                              n = SUM list(j)*base**(j)
!                                  j=1
!
      implicit none
      integer, dimension(:) :: list
      integer, INTENT(IN) :: base
      integer, INTENT(IN OUT) :: nnzero
      integer j
!                   if nnzero is goofy, then do nothing
      if( nnzero .lt. 0 ) return
!                   increment first digit
      j = 1
      do            ! repeat loop without condition
      if( j .gt. nnzero ) nnzero = j 
      list(j) = list(j) + 1
!                   if no carry, then done
      if( list(j) .lt. base ) exit
!                   carrying
      list(j) = 0
      j = j + 1
      end do        ! end loop
      return
      end subroutine ncrmnt

      real function expand(list,base,nnzero)
!  computes the real number corresponding to its base "base"
!   expansion to nnzero digits listed in list
!
!                                 nnzero
!   x = expand(list,base,nnzero) = SUM list(j)*base**(-j)
!                                  j=1
!   arguments
!  list     integer vector of length at least nnzero holding digit
!            of base "base" expansion -- dummy dimensioned as length 1
!  base     integer base
!  nnzero   integer giving the number of nonzero digits
!
! **  note that all digits are to the right of the radix point **
! **  so that expand(list,base,nnzero) is between 0 and 1      **
!
      implicit none
      integer, dimension(:) :: list
      integer, INTENT(IN) :: base
      integer, INTENT(IN OUT) :: nnzero
      integer j,jj
      real flbase
!                   get floating point value of base
      flbase = real(base)
!                   initialize value to zero
      expand = 0.
!                   increment j through number of digits given
      do j = 1,nnzero
!                   do the least significant first
      jj = nnzero + 1 - j
      expand = list(jj) + expand/flbase
      end do  ! loop on j
!                   last divide
      expand = expand / base
      return
      end function expand
! *** end of filename halton.f95        *********************