/*****************************************************************************
*                                                                            *
*   NCHOOSEK.SPL Copyright (C) 2014 DSP Development Corporation              *
*                            All Rights Reserved                             *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Computes binomial coefficient or returns all combinations   *
*                                                                            *
*   Revisions:    8 Jun 2014  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_NCHOOSEK

    NCHOOSEK

    Purpose: Computes the binomial coefficient or returns all the
             unordered combinations of the elements of the input series.

    Syntax:  NCHOOSEK(n, k)

              n - A scalar or series. For a scalar, the binomial coefficient
                  is returned. For a series, all the unordered combinations
                  of the series elements are returned for k a positive
                  integer.

              k - A scalar or series. If n is a series and k a positive
                  integer less than length(n), all the unordered combinations
                  of the series elements are returned.

    Returns: A scalar or series, the number of combinations of n items
             taken in unordreed groups of k, or if n is a series and k
             a positive integer, the actual unordered combinations of
             the series elements where each row represents a group of k
             elements.

    Example:
             nchoosek(5, 3)

             returns 10

    Example:
             c = nchoosek(1..5, 3)

             c == {{1, 2, 3},
                   {1, 2, 4},
                   {1, 2, 5},
                   {1, 3, 4},
                   {1, 3, 5},
                   {1, 4, 5},
                   {2, 3, 4},
                   {2, 3, 5},
                   {2, 4, 5},
                   {3, 4, 5}}

    Example:
             nchoosek(5.3, 3)

             returns 12.5345

    Example:
             nchoosek(3..5, 1..3)

             returns {3, 6, 10}, the same as:

             {nchoosek(3, 1), nchoosek(4, 2), nchoosek(5, 3)}

    Remarks:
             To return the actual combinations, N must be a series and
             K must be an integer where 0 <= k <= length(n).

             For scalar inputs N and K can be any value.

             For integer values where 0 <= k <= n:

                                n         n!
             nchoosek(n, k) = (   ) = ----------- 
                                k     k! (n - k)!


             See BINCOEFF to return the generalized binomial coefficient
             for any scalar or series inputs.

    See Also:
             Bincoeff
             Factorial
             Pascal
#endif


/* from n items, choose groups of k */
ITERATE nchoosek(v, k)
{
        local c, i, m, n, p, r, cp, rp, nk;

        if (argc < 2)
        {
                error("nchoosek - 2 input values required");
        }

        n = length(v);

        if (n == 1 || length(k) > 1)
        {
                /* binomial coefficient for scalar or single value series */
                return(bincoeff(v, k));
        }
        else if (floor(k) != real(k) || k < 0)
        {
                error("nchoosek - k must be a non-negative integer for n a series");
        }

        /* ensure integer k */
        k = castint(k);

        /* v is a series and k a positive integer */
        if (k == 0 || k > n)
        {
                p = {};
        }
        else if (k == 1)
        {
                /* n x 1 */
                p = v;
        }
        else if (k == 2)
        {
                r = repelem(v, n-1..-1..1);
                p = serextract(v, 2..n, n-1..-1..1);
                p = ravel(r, p);
        }
        else if (k == n)
        {
                /* 1 x n */
                p = v';
        }
        else
        {
                /* iterate combinations back to front */
                nk = 1..n-k+1;
                rp = {};
                p  = v[k..n];

                c = m = nk;

                loop (i = 2..k)
                {
                        cp = concat(rp, serextract(p, m, {-1}));

                        c  = 1 + length(p) - c;
                        r  = repelem(v[k-i+1..n-i+1], nk, c);
                        p  = ravel(r, cp);

                        c  = partsum(c);
                        m  = 1 + c[1..n-k];
                        c  = {1, m};

                        /* for speed, reference array, do not copy */
                        rp = refseries(p);
                }
        }

        /* set to table */
        setplotstyle(p, 4);

        return(p);
}