/*****************************************************************************
*                                                                            *
*   PASCAL.SPL   Copyright (C) 2014 DSP Development Corporation              *
*                            All Rights Reserved                             *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates Pascal's triangle                                *
*                                                                            *
*   Revisions:    3 Jul 2014  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_PASCAL

    PASCAL

    Purpose: Calculates Pascal's triangle.

    Syntax:  PASCAL(n, order)

                  n -  An integer.

              order - Optional. An integer, the matrix order.

                       0: lower triangle of the standard Pascal
                          matrix (default)

                       1: lower triangle Cholesky factor of the Pascal
                          matrix.

                       2: the lower triangle Cholesky factor of the transposed
                          and permuted Pascal matrix.

    Returns: An array.

    Example:
             pascal(4)

             returns  {{1,  1,  1,   1},
                       {1,  2,  3,   4},
                       {1,  3,  6,  10},
                       {1,  4, 10,  20}}

    Example:
             W1: pascal(4, 1)
             W2: W1^^2
             W3: eye(4)
             W4: inv(W1)

             W1 == {{1,  0,  0,  0}
                    {1, -1,  0,  0}
                    {1, -2,  1,  0}
                    {1, -3,  3, -1}}

             W2 == {{1,  0,  0,  0}
                    {0,  1,  0,  0}
                    {0,  0,  1,  0}
                    {0,  0,  0,  1}}

             W3 == {{1,  0,  0,  0}
                    {0,  1,  0,  0}
                    {0,  0,  1,  0}
                    {0,  0,  0,  1}}

             W4 == {{1,  0,  0,  0}
                    {1, -1,  0,  0}
                    {1, -2,  1,  0}
                    {1, -3,  3, -1}}

             For order 1, the result is equal to the square root of the
             identity matrix. The result is also equal to the inverse itself.

    Example:
             W1: pascal(4, 2)
             W2: W1^^3

             W1 == {{-1, -1, -1, -1}
                    { 3,  2,  1,  0}
                    {-3, -1,  0,  0}
                    { 1,  0,  0,  0}}

             W2 == {{ 1,  0,  0,  0}
                    { 0,  1,  0,  0}
                    { 0,  0,  1,  0}
                    { 0,  0,  0,  1}}

             For order 2, the result is equal to the cube root of the
             identity matrix.

    Remarks: 
             PASCAL returns the lower triangle of Pascal's matrix.

             For ORDER == 1, the resulting array is the lower triangle
             of the Cholesky factor of the Pascal matrix. The result is
             its own inverse and is equal to the square root of the
             identity matrix, thus:

               pascal(n, 1)^^2 == eye(n)

             For ORDER == 2, the result is the lower triangle of a
             the Cholesky factor transposed and permuted such that the
             result is equal to the cube root of the identity matrix.
             Thus:

               pascal(n, 2)^^3 == eye(n)

    See Also:
             Bincoeff
             Cholesky
             Nchoosek
#endif


/* Pascal's Triangle */
pascal(n, order)
{
        local p;

        if (argc < 2)
        {
                if (argc < 1)
                {
                        error("pascal - input must be integer values");
                }

                order = 0;
        }

        n = int(n);

        if (n < 0)
        {
                p = {};
        }
        else if (order == 0)
        {
                /* lower triangle */
                p = ravel(serextract(0..(2*n-2), 1..n, {n}), n);
                p = bincoeff(p, 0..n-1);
        }
        else if (order > 0)
        {
                /* Cholesky factor */
                p = rev(ravel(serextract(-n..(n-2), 1..n, {n}), n));
                p = bincoeff(p, rep((0..n-1)', n));

                if (order > 1)
                {
                        /* transpose permuted */
                        p = (-1)^(n + 1) * rev(p)';
                }
        }

        return(p);
}