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)