/*****************************************************************************
*                                                                            *
*   BINCOEFF.SPL Copyright (C) 2014 DSP Development Corporation              *
*                            All Rights Reserved                             *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the binomial coefficient                         *
*                                                                            *
*   Revisions:   22 Jul 2014  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_BINCOEFF

    BINCOEFF

    Purpose: Calculates the generalized binomial coefficient.

    Syntax:  BINCOEFF(n, k)

              n - A scalar or series.

              k - A scalar or series.

    Returns: A scalar or series. For positive integer values, the number 
             of combinations of n items taken in unordered groups of k.

    Example:
             bincoeff(5, 3)

             returns 10

    Example:
             bincoeff(5.3, 3)

             returns 12.5345

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

             c == {0, 0, 1, 4, 10}

    Example:
             c = bincoeff(3..7, 1..5)

             c == {3, 6, 10, 15, 21}

    Example:
             c = bincoeff(3..0.5..7, 1..0.5..5)

             c == {3, 4.375, 6, 7.875, 10, 12.375, 15, 17.875, 21}

    Remarks:
             BINCOEFF returns the generalized binomial coefficient. The
             inputs can be any series or scalar values.

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

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

             The result is accurate to 15 digits. A warning is
             displayed if the computation results in an overflow.

             See NCHOOSEK to return the actual combinations of series
             elements for a given input series.

    See Also:
             Factorial
             Gamma
             Nchoosek
             Pascal
#endif


/* binomial coefficient, n choose k */
bincoeff(n, k)
{
        local c, overflow;

        if (argc < 2)
        {
                if (argc < 1)
                {
                        n = 1;
                }

                k = 1;
        }

        /* core computation */
        (c, overflow) = bincoeff_core(n, k);

        if (outargc > 1)
        {
                return(c, any(overflow));
        }
        else
        {
                if (any(overflow))
                {
                        message("Warning", "bincoeff - possible overflow, result only accurate to 15 digits", 8);
                }

                return(c);
        }
}


/* binomial coeffcient core routine */
ITERATE bincoeff_core(n, k)
{
        local c, ni, ki, idx, nlz, ngz, ibin, intn, intk, overflow = 0;

        /* sort out sizes */
        if (isarray(n) && not(isarray(k)))
        {
                k = ones(length(n), 1, deltax(n)) * k;
        }
        else if (isarray(k) && not(isarray(n)))
        {
                n = ones(length(k), 1, deltax(k)) * n;
        }

        if (length(n) != length(k))
        {
                error("bincoeff - inputs must be of conforming size");
        }

        /* initialize to zeros of same length */
        c = n * 0;

        /* integer values */
        intn = floor(n) == real(n);
        intk = floor(k) == real(k);

        /* n < 0 */
        nlz = n < 0;

        if (any(nlz))
        {
                /* n < 0, non-integer n */
                idx = find(nlz && not(intn));
                ni  = n[idx];
                ki  = k[idx];

                c[idx] = bincoeff_bico(ni, ki);

                /* n < 0, integer n, integer k, k >= 0 */
                idx = find(nlz && intn && intk && k >= 0);
                ni  = n[idx];
                ki  = k[idx];

                c[idx] = (-1)^(ki) * bincoeff_bico(-ni + ki - 1, ki);

                /* n < 0, k <= n, integer n, integer k */
                idx = find(nlz && intn && intk && k <= n);
                ni  = n[idx];
                ki  = k[idx];

                c[idx] = (-1)^(ni - ki) * bincoeff_bico(-ki - 1, ni - ki);

                /* n < 0, integer n, non-integer k */
                idx    = find(nlz && intn && not(intk));
                c[idx] = inf;
        }

        ngz = n >= 0;

        if (any(ngz))
        {
                /* use integer binomial here */
                ibin = intn && intk && (k >= 0) && (k <= n);

                /* n > 0, k > 0 && k < n && real integer */
                idx = find(ngz && ibin);
                ni  = n[idx];
                ki  = k[idx];

                /* positive integer binomial coefficient */
                c[idx] = binomial(ni, ki);

                /* n > 0 and k anything */
                idx = find(ngz && not(ibin));
                ni  = n[idx];
                ki  = k[idx];

                c[idx] = bincoeff_bico(ni, ki);
        }

        /* nan */
        idx    = find(isnan(n) || isnan(k));
        c[idx] = nan;

        /* fast check if all values real */
        if (isreal(n) && isreal(k))
        {
                c = real(c);
        }
        else
        {
                /* some values might be real */
                idx    = find(imag(n) == 0 && imag(k) == 0);
                c[idx] = real(c);
        }

        /* handle integer results */
        idx    = find(intn && intk);
        c[idx] = bincoeff_round(c[idx]);

        /* check for overflow */
        if (any(c[idx] > 1e15))
        {
                overflow = 1;
        }

        return(c, overflow);
}


/* binomial coefficient using gammaln */
bincoeff_bico(n, k)
{
        return(exp(gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1)));
}


/* rounding that handles large numbers */
bincoeff_round(x)
{
        return(floor(x) == x ? x : floor(x + 0.5));
}