/* IO(1) */

matrix C( 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;

  m = cols(v);

  if( narg != 2 || m < 1 || n < 1 || m < n) {
    print( "C: 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 = C( 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;
}
