poly.src

来自「没有说明」· SRC 代码 · 共 281 行

/*
** poly.src - polynomial functions
** (C) Copyright 1988-1998 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**  Format                  Purpose                          Line
** =================================================================
** c = POLYCHAR(x);         characteristic polynomial          27
** y = POLYEVAL(x,c);       evaluate polynomial                53
** c = POLYMAKE(r);         coefficients of polynomial        133
** c = POLYMULT(c1,c2);     multiply polynomials              170
** r = POLYROOT(c);         roots of polynomial               201
*/

/*
**> polychar
**
**  Purpose:    To compute the characteristic polynomial of a square
**              matrix.
**
**  Format:     c = polychar(x);
**
**  Input:      x    NxN matrix.
**
**  Output:     c    N+1x1 vector of coefficients of the Nth order
**                   characteristic polynomial of x:
**
**                   p(z)=c[1,1]*z^n + c[2,1]*z^(n-1) + ... +
**                   c[n,1]*z + c[n+1,1];
**
**  Remarks:    The coefficient of z^n is set to unity (c[1,1]=1).
**
**  See Also:   polymake, polymult, polyroot, polyeval
*/


proc polychar(x);
    retp( polymake( eig( x ) ) );
endp;

/*
**> polyeval
**
**  Purpose:    To evaluate polynomials. Can either be 1 or more
**              scalar polynomials, or a single matrix polynomial.
**
**  Format:     y = polyeval(x,c);
**
**  Input:      x    1xK or NxN; that is, x can either represent K
**                   separate scalar values at which to evaluate the
**                   (scalar) polynomial(s), or it can represent a
**                   single NxN matrix.
**
**              c    P+1xK or P+1x1 matrix of coefficients of
**                   polynomials to evaluate. If x is 1xK, then c
**                   must be P+1xK. If x is NxN, c must be P+1x1.
**                   That is, if x is a matrix, it can only be
**                   evaluated at a single set of coefficients.
**
**  Output:     y    Kx1 vector (if c is P+1xK) or NxN matrix (if c
**                   is P+1x1 and x is NxN):
**
**                   y = ( c[1,.].*x^p + c[2,.].*x^(p-1)  + ... +
**                   c[p+1,.] )';
**
**  Remarks:    In both the scalar and the matrix case, Horner's
**              rule is used to do the evaluation. In the scalar
**              case, the function recsercp is called (this
**              implements an elaboration of Horner's rule).
**
**  Example:    x = 2; let c = 1 1 0 1 1;
**              y = polyeval(x,c);
**
**              The result is 27. Note that this is the decimal
**              value of the binary number 11011.
**
**                   y = polyeval(x,1|zeros(n,1));
**
**              This will raise the matrix x to the nth power (e.g:
**              x*x*x*x*...*x).
**
**  See Also:   polymake, polychar, polymult, polyroot
*/

proc polyeval(x,c);
    local c1,y,rx,cx,rc,cc,i,p;

    rx = rows(x);
    cx = cols(x);
    rc = rows(c);
    cc = cols(c);

    if rx == 1;     /* scalar polynomial(s) */
        c1 = c[1,.];
        if not(c1 /= 0);    /* some 0's in first row of c */
            c1 = missrv(miss(c1,0),1);      /* replace 0's with 1's  */
        endif;
        p = rows(c) - 1;
        y = recsercp(x,trimr(c./c1,1,0));
        retp( conj((c1.*y[p,.] - x^p.*(c[1,.] .== 0))') );

    elseif rx > 1 and rx == cx and cc == 1;         /* single matrix
                                                    :: polynomial
                                                    */

        y = c[1];
        i = 2;
        do until i > rc;
            y = y * x + diagrv(zeros(rc,rc),c[i]);
            i = i + 1;
        endo;
        retp( y );

    else;
        errorlog "ERROR: Input matrices are not the right size in POLYEVAL.";
        end;
    endif;

endp;

/*
**> polymake
**
**  Purpose:    To compute the coefficients of a polynomial, given
**              the roots of the polynomial. (Restricted to real
**              roots).
**
**  Format:     c = polymake(r);
**
**  Input:      r    Nx1 vector containing roots of the desired
**                   polynomial.
**
**  Output:     c    N+1x1 vector containing the coefficients of the
**                   Nth order polynomial with roots r:
**
**                   p(z)=c[1]*z^n + c[2]*z^(n-1) + ... +
**                     c[n]*z + c[n+1];
**
**  Remarks:    The coefficient of z^n is set to unity (c[1]=1).
**
**  See Also:   polychar, polymult, polyroot, polyeval
*/

proc polymake(r);
    local n,c,j;
    n = rows(r);
    c = zeros(1,n+1);
    j = 1;
    c[1] = 1;
    do until j > n;
        c = c - r[j] * shiftr(c,1,0);
  /*      c[2:j+1] = c[2:j+1] - r[j]*c[1:j]; */
        j = j+1;
    endo;
    retp(c');
endp;

/*
**> polymult
**
**  Purpose:    To multiply two polynomials together.
**
**  Format:     c = polymult(c1,c2);
**
**  Input:      c1   d1+1x1 vector containing the coefficients of
**                   the first polynomial.
**
**              c2   d2+1x1 vector containing the coefficients of
**                   the second polynomial.
**
**  Output:     c    d1+d2x1 vector containing the coefficients of
**                   the product of the two polynomials.
**
**  Remarks:    If the degree of c1 is d1 (eg, if d1=3, then the
**              polynomial corresponding to c1 is cubic), then there
**              must be d1+1 elements in c1 (e.g. 4 elements for a
**              cubic). Thus, for instance the coefficients for  the
**              polynomial 5*x.^3 + 6*x + 3 would be: c1=5|0|6|3.
**              (Note that zeros must be explicitly given if the
**              are powers of x missing.)
**
**  See Also:   polymake, polychar, polyroot, polyeval
*/

proc polymult(c1,c2);
    retp( conv(c1,c2,0,0) );
endp;

/*
**> polyroot
**
**  Purpose:    To compute the roots of a polynomial whose coefficients are c.
**
**  Format:     r = polyroot(c);
**
**  Input:      c      N+1x1 vector of coefficients of an Nth order
**                     polynomial:
**
**                     p(z)=c[1]*z^n + c[2]*z^(n-1) + ... + c[n]*z + c[n+1];
**
**                     Zero leading terms will be stripped from c.  When that
**                     occurs the order of the r will be the order of the
**                     polynomial after the leading zeros have been stripped.
**
**              c[1]   need not be normalized to unity.
**
**  Output:     r      Nx2 vector containing the roots of c.
**
**  See Also:   polymake, polychar, polymult, polyeval
*/

#ifcplx

    proc polyroot(c);

#else
    #ifdef OLD_POLYROOT

    proc polyroot(c);

    #else

    proc (2)=polyroot(c);

    #endif
#endif

    local n,a;

#ifreal
    #ifdef OLD_POLYROOT

    local  rl,im;

    #endif
#endif

       do while c[1] == 0;
         if rows(c) <= 1;
               retp(error(1));
         else;
               c=c[2:rows(c)];
         endif;
       endo;

       n = rows(c)-1;
       if n == 1;
              retp( -c[2]/c[1] );
       else;
              a = ( -trimr(c,1,0)'/c[1,1] )| eye(n-1)~zeros(n-1,1);
#ifcplx

              retp( eig(a) );

#else
    #ifdef OLD_POLYROOT

              { rl,im } = cmsplit(eig(a));
              retp( rl~im );

    #else

             retp( cmsplit(eig(a)) );

    #endif
#endif

       endif;
endp;