mata
real matrix queens_rec(real scalar n) {
	real rowvector a, u, v
	real matrix m
	
	a = 1..n
	u = J(1, 2*n-1, 1)
	v = J(1, 2*n-1, 1)
	m = J(0, n, .)
	queens_aux(n, 1, a, u, v, m)
	return(m)
}

void queens_aux(real scalar n, real scalar i, real rowvector a,
	real rowvector u, real rowvector v, real matrix m) {
	real scalar j, k
	
	if (i > n) {
		m = m\a
	} else {
		for (j = i; j <= n; j++) {
			k = a[j]
			p = i-k+n
			q = i+k-1
			if (u[p] & v[q]) {
				u[p] = v[q] = 0
				a[j] = a[i]
				a[i] = k
				queens_aux(n, i+1, a, u, v, m)
				u[p] = v[q] = 1
				a[i] = a[j]
				a[j] = k
			}
		}
	}
}
end