*! version 1.0.0  12jul2005 

mata


real matrix halton(real scalar n, d, |real scalar burn, real scalar method)
{
	real matrix x

	x = J(n,d,.)
	if (burn >= .) burn = 0
	if (method >= .) method = 1
	_halton(x, burn, method)

	return(x)
}

void _halton(real matrix x, |real scalar burn, real scalar method)
{
	real scalar n, d, d1, l, j, k, dn, xj, eps
	real rowvector xk
	real rowvector g_pr

	g_pr = (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71)

	n = rows(x)
	d = cols(x)
	if (d > length(g_pr)) {	
		errprintf("requested dimension of %g; maximum dimension is %g\n", 
			d, length(g_pr))
		exit(198)
	}

	if (burn < .) burn = abs(trunc(burn))
	else burn = 0

	if (method >= .) method = 1
	else if (method!=1 && method!=2) {
		errprintf("requested method %g; method must be 1=halton or 2=hammersly\n",
			method)
		exit(198)
	}

	eps = 100*st_numscalar("c(epsdouble)")
	dn = xj = 0
	d1 = d 
	xk = J(1,d,0)
	/* burn in the sequence */
	if (method == 2) {
		/* Hammersley */
		dn = 1/n
		xj = halton_mod1((2*burn-1)/(2*n))
		xk[1] = xj

		--d1;
	}
	x = J(n,d,0)

	for (k=1; k<=d1; k++) {
		pr = g_pr[k]
		pj = 1
		xj = 0
		/* burn in the sequence */
		j = burn
		while (j>0) {
			pj = pj/pr
			b = mod(j,pr)
			xj = xj + pj*b
			j = (j-b)/pr
		}
		if (method == 2) {
			xk[k+1] = xj
		}
		else {
			xk[k] = xj
		}
	}
	for (l=1; l<=n; l++) {
		for (j=1; j<=d; j++) {
			if (method == 2) {
				/* Hammersley */
				if (j == 1) {
					xj = halton_mod1(xk[j]+dn)
					pr = 0
				}
				else {
					pr = g_pr[j-1]
				}
			}
			else {
				pr = g_pr[j]
			}
			if (pr > 0) {
				pj = 1
				xj = 1 - xk[j] - eps	
				while (pj >= xj) {
					pj = pj/pr
				}
				xj = xk[j]+(pr+1)*pj-1
			}
			xk[j] = xj
		}
		x[l,.] = xk
	}
}

real scalar function halton_mod1(real scalar a)
{
	return(a-floor(a))
}

end