function x = halton(n)
% Halton Quasi-Random numbers
%
% to initialize
%    halton(n)
% to generate next vector x of length n
%    x = halton
%
% only for demonstration purposes (very slow, and other Quasi-Random number
% generators like Sobol's have better properties)
global dig pri nn ndig

if nargin == 1
  dig = zeros(n,40); % assuming fewer than 2^40 random numbers are needed
  pri = primes(100); % assuming n is less than 26
  nn = n;
  ndig = zeros(1,n);
else
  x = zeros(1,nn);
  for i=1:nn
    k = 1;
    while 1
      if dig(i,k) < pri(i)-1
        dig(i,k) = dig(i,k)+1;
        if k > ndig(i)
          ndig(i) = k;
        end
        break;
      else
        dig(i,k) = 0;
        k = k+1;
      end
    end
    q = 1;
    for j=1:ndig(i)
      q = q/pri(i);
      x(i) = x(i) + dig(i,j)*q;
    end
  end
end