/******************************************
**    nCr(n,r)
**
**    Written by Tsz-Cheong Lai, July 2002
**
**    This code is provided gratis without any guarantees or
**    warrantees. Proprietary modifications of this code are
**    not permitted.
**
**    Function: compute all nCr combinations
**    Inputs: n -- the total number of items
**            r -- the number of items to choose from the n
**
**    Output: nCr x r matrix
**            where each row is a combination
**
**    e.g.       nCr(3,2) = 1  2
**                          1  3
**                          2  3
**
**    Auxillary function: _serial_sum
******************************************/

PROC nCr(n,r);
LOCAL errout, errmsg, start_n, mtrx, idx, rows_to_add, idx_n;

IF n <= 0;
        GOTO errout("* ERROR: n = " $+ FTOS(n,"*.*lf",CEIL(LOG(n))+1,2) $+ " is negative.");
ELSEIF r <= 0;
        GOTO errout("* ERROR: r = " $+ FTOS(r,"*.*lf",CEIL(LOG(r))+1,2) $+ " is negative.");
ELSEIF n /= TRUNC(n);
        GOTO errout("* ERROR: n = " $+ FTOS(n,"*.*lf",CEIL(LOG(n))+1,2) $+ " is not an integer.");
ELSEIF r /= TRUNC(r);
        GOTO errout("* ERROR: r = " $+ FTOS(r,"*.*lf",CEIL(LOG(r))+1,2) $+ " is not an integer.");
ELSEIF r > n;
        GOTO errout("* ERROR: r = " $+ FTOS(r,"*.*lf",CEIL(LOG(r))+1,2) $+ " is larger than n = " $+ FTOS(n,"*.*lf",CEIL(LOG(n))+1,2) $+ ".");
ELSEIF coreleft/(8*r*n!/(r!*(n-r)!)) < 2;
        GOTO errout("* ERROR: Insufficient workspace for n = " $+ FTOS(n,"*.*lf",CEIL(LOG(n))+1,2) $+ " and r = " $+ FTOS(r,"*.*lf",CEIL(LOG(r))+1,2) $+ ".");
ELSEIF r == 1;
        RETP(SEQA(1,1,n));
ELSEIF r == n;
        RETP(SEQA(1,1,n)');
ENDIF;

/* Create aCb where a = n - r + 2 and b = 2 */
start_n = n - r + 2;
mtrx = ONES(start_n - 1,1) ~ SEQA(2,1,start_n - 1);
idx = 1;
DO UNTIL idx == start_n - 1;
        mtrx = mtrx | (mtrx[1:start_n-1-idx,.] + idx);
        idx = idx + 1;
ENDO;

/* The computation is already done if r = 2 */
IF r == 2;
        RETP(mtrx);
ENDIF;

/* Move from aCb to (a+1)C(b+1) until we reach nCr */
rows_to_add = SEQA(1,1,start_n - 2);
idx_n = start_n;
DO UNTIL idx_n == n;
        /* The matrix aCb is flipped over to form the upper-right part of
        ** the new matrix (a+1)C(b+1)
        ** which contains all the combinations with the first item chosen */
        mtrx = (idx_n + 1 - mtrx);
        mtrx = rev(rev(mtrx')');
        mtrx = ONES(ROWS(mtrx), 1) ~ (1 + mtrx);

        /* The remaining combinations with the second item chosen are added,
        ** then the third, so on and so forth */
        rows_to_add = _serial_sum(rows_to_add);
        idx = 0;
        DO UNTIL idx == ROWS(rows_to_add);
                mtrx = mtrx | (mtrx[1:rows_to_add[ROWS(rows_to_add)-idx],.] + idx + 1);
                idx = idx + 1;
        ENDO;

        idx_n = idx_n + 1;
ENDO;
RETP(mtrx);

errout:
POP errmsg;
ERRORLOG "* Procedure \34nCr\34 is in use.";
ERRORLOG errmsg;
STOP;
ENDP;

PROC _serial_sum(vctr);
LOCAL idx;
idx = 2;
DO UNTIL idx > ROWS(vctr);
        vctr[idx] = vctr[idx-1] + vctr[idx];
        idx = idx + 1;
ENDO;
RETP(vctr);
ENDP;