/* IO(1) */

matrix combinations( vector _v, int _n)
{
/* generate combinatorial matrix of m things taken n at a time, m >= n */

  matrix c, cnext;
  vector vnext;
  int i, k, R;
  int m;

  if( narg != 2 || (m = cols(_v)) < 1 || _n < 1 || m < _n) {
    print( "combinations: bad arguments");
    return 0;
  }

  if( _n == 1)	/* R combinations are just _v' itself */
    return _v';

  if( _n == m)	/* 1 combination is just _v itself */
    return _v;

  R = prod(_n+1:m) / prod(1:m-_n); /* number of combinations == rows(c) */

  c = zeros( R, _n);
  vnext = _v;
  i = 1;

  while( i <= R) {
    cnext = combinations( vnext(2:m), _n-1);
    k = rows(cnext);
    c(i:i+k-1) = [reshape(vnext(1),k,1), cnext];
    vnext = vnext(2:m);
    --m;
    i += k;
  }

  return c;
}
