BEGIN_PROVIDER [ double precision, halton_seed, (3) ]
  implicit none
  BEGIN_DOC
  ! Seed of the quasi Monte Carlo sequence generator
  END_DOC
  
  halton_seed(:) = 1.d0
END_PROVIDER

BEGIN_PROVIDER [ double precision, halton_base, (3) ]
  implicit none
  BEGIN_DOC
  ! Base of the quasi Monte Carlo sequence generator
  END_DOC
  
  halton_base(1) = 2.d0
  halton_base(2) = 3.d0
  halton_base(3) = 5.d0
END_PROVIDER

double precision function halton_ranf(axis)
  implicit none
  integer, intent(in)            :: axis
  
!call random_number(halton_ranf)
!return
  ASSERT (axis > 0)
  ASSERT (axis <= 3)
  
  double precision :: l_seed
  l_seed = halton_seed(axis)
  halton_seed(axis) += 1_8
  
  halton_ranf = 0.d0
  
  double precision               :: base_inv, base_inv0
  base_inv  = 1.d0/halton_base(axis)
  base_inv0 = 1.d0/halton_base(axis)
  
  do while (l_seed > 0.d0)
    halton_ranf = halton_ranf + dble(mod(int(l_seed),int(halton_base(axis)))) * base_inv
    base_inv *= base_inv0
    l_seed *= base_inv0
  enddo
end