/*****************************************************************************
*                                                                            *
*   POLY.SPL       Copyright (C) 1997-2015 DSP Development Corporation       *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Carl Trapanni                                               *
*                                                                            *
*   Synopsis:    Calculates coefficents of characteristic polynomial         *
*                                                                            *
*   Revisions:   18 Sep 1997  CCT  Creation                                  *
*                 1 Oct 1997  CCT  complex input                             *
*                17 Feb 2005  RRR  use convolution for polynomial multiply   *
*                                                                            *
*****************************************************************************/


#if @HELP_POLY

    POLY

    Purpose:
             Calculates coefficients of the characteristic polynomial

    Syntax:
              POLY(x)

              x - A series or matrix. A matrix where the eigenvalues are
                  the roots of the returned polynomial coefficients or
                  a series of the roots of the returned polynomial
                  coefficients.

    Returns:
              A series, the polynomial coefficients in descending order..

    Example:
              W1: {{1, 2, 3},
                   {3, 4, 5},
                   {5, 6, 0}}

              W2: poly(W1)
              W3: polyroot(W2, 1)
              W4: eigval(W1)


              W2 contains the series {1, -5, -47, -14}, the coefficients in
              descending order of the characteristic polynomial of the matrix
              in W1. The roots of this polynomial are the eigenvalues of the
              input matrix.

              W3 == {9.894, -4.585, -0.309}
              W4 == {9.894, -0.309, -4.585}


              If the input is a series, the series represents the roots of
              the polynomial coefficients returned by POLY. For example,
              consider the polynomial:

              P(x) = (x + 2) * (x + 1) = x^2 + 3*x + 2

              poly({-2, -1}) == {1, 3, 2}

              polyroot({1, 3, 2}, 1) == {-2, -1}

    Remarks:
              As indicated above, if the input is a matrix, POLY returns the
              coefficients of the characteristic polynomial in descending
              order.

              If the input is a series, the input represents the roots of a
              polynomial and POLY returns the coefficients of this polynomial
              in descending order.

    See Also:
              POLYFIT
              POLYDER
              POLYROOT
              POLYVAL
              ROOTS

#endif


/* find characteristic polynomial or polynomial from roots */
poly(x)
{
        local a, k, m, n, v, y, cj1, cj2;

        (m, n) = size(x);

        /* empty case first since EIGVAL won't accept empty input */
        if (m == 0)
        {
                return({1});
        }

        else if (n == 1)
        {
                /* scalar or use faster reference */
                v = (m == 1) ? {x} : refseries(x);
        }

        /* convert single row into column */
        else if ((m == 1) && (n > 1))
        {
                v = x';
        }

        /* characteristic polynomial if square x */
        else if (m == n)
        {
                /* eigenvalues are the roots of the characteristic polynomial */
                v = eigval(x);
                
                if (length(v) == 0)
                {
                        return({1});
                }
        }
        else
        {
                error("poly - series or square matrix required");
        }

        /* strip infinites */
        v = delete(v, isinf(v));

        if (length(v) == 0)
        {
                /* empty or all inf */
                return({1});
        }

        /* create array where each column is (1 - ROOT_k) */
        a        = ones(2, length(v));
        a[2, ..] = -v';

        y = {1};

        /* use convolution to multiply root terms */
        loop (k = 1..numcols(a))
        {
                y = conv(y, col(a, int(k)));
        }

        if (not(iscomplex(x)))
        {
                /* force real if matrix was real */
                y = real(y);
        }
        else
        {
                /* check for all complex conjugate roots */
                cj1 = sort(v[find(imag(v) > 0)]);
                cj2 = sort(conj(v[find(imag(v) < 0))]);

                if (length(cj1) == length(cj2))
                {
                        if (all(cj1 == cj2))
                        {
                                y = real(y);
                        }
                }
        }

        return(y);
}