home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / poly < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.7 KB  |  78 lines
  1. //----------------------------------------------------------------------
  2. //
  3. // poly
  4. //
  5. // Syntax: c=poly(x)
  6. //
  7. // This routine evaluates the characteristic polynomial contained
  8. // in x. If x is an N by N matrix, then poly(x) will return a row
  9. // vector of length N+1 whose elements are the coefficients of the
  10. // characteristic polynomial (det(lambda*eye(N,N)-a)).
  11. //
  12. // If x is a vector, then poly(x) will return a vector whose elements
  13. // are the coefficients of the polynomial whose roots are the elements
  14. // of x.
  15. //
  16. // For example,
  17. //
  18. // > rfile poly
  19. // > A=[0,4,5,6];
  20. // > poly(A)
  21. //         1        -15         74       -120          0  
  22. //
  23. //
  24. // > rfile poly
  25. // > A=[1,2,3;4,5,6;7,8,0];
  26. // > poly(A)
  27. //         1         -6        -72        -27  
  28. //
  29. // In general, the routines poly and roots are the inverse of each
  30. // other (But only for vector inputs).
  31. //
  32. // Ref: MATLAB Reference Guide, The Mathworks Inc. August 1992.
  33. //
  34. // Copyright (C), by Jeffrey B. Layton, 1994-95
  35. // Version JBL 950105
  36. //----------------------------------------------------------------------
  37.  
  38. poly = function(x)
  39. {
  40.    local(n,z,i,c)
  41.  
  42. // Check for the empty case
  43.    if (max(size(x)) == 0) {
  44.        c=1;
  45.        return c
  46.    }
  47.  
  48. // find eigenvalues of x
  49.    if (x.nc == x.nr) {
  50.        z=eig(x).val;
  51.    else if (x.nr != 1) {
  52.        z=x.';
  53.    else
  54.        z=x;
  55.    }}
  56.  
  57. // strip away infinities and NaN
  58.    z=z[finite(z).*[1:z.nc]];
  59.  
  60.    n=max([z.nr,z.nc]);
  61.  
  62. // Initialize c (characteristic polynomial)
  63.    c=zeros(1,n+1);
  64.    c[1]=1;
  65.  
  66. // Expand recursion formula
  67.    for (i in 1:n) {
  68.        c[2:(i+1)]=c[2:(i+1)]-z[i].*c[1:i];
  69.    }
  70.  
  71. // If x was real, then make c real.
  72.    if (x.type == "real") {
  73.        c=real(c);
  74.    }
  75.  
  76.    return c
  77. };
  78.