function combinations(n,k) {
	a = J(comb(n,k),k,.)
	u = 1..k
	for (i=1; 1; i++) {
		a[i,.] = u
		for (j=k; j>0; j--) {
			if (u[j]-j<n-k) break
		}
		if (j<1) return(a)
		u[j..k] = u[j]+1..u[j]+1+k-j
	}
}

combinations(5,3)