File: poly.m

package info (click to toggle)
freemat 4.0-5
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd, wheezy
  • size: 174,736 kB
  • ctags: 67,053
  • sloc: cpp: 351,060; ansic: 255,892; sh: 40,590; makefile: 4,323; perl: 4,058; asm: 3,313; pascal: 2,718; fortran: 1,722; ada: 1,681; ml: 1,360; cs: 879; csh: 795; python: 430; sed: 162; lisp: 160; awk: 5
file content (49 lines) | stat: -rw-r--r-- 1,257 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
% POLY POLY Convert Roots To Polynomial Coefficients
% 
% Usage
% 
% This function calculates the polynomial coefficients for given roots
% 
%     p = poly(r)
% 
% when r is a vector, is a vector whose elements are the coefficients 
% of the polynomial whose roots are the elements of r.  Alternately,
% you can provide a matrix
% 
%     p = poly(A)
% 
% when A is an N x N square matrix, is a row vector with 
% N+1 elements which are the coefficients of the
% characteristic polynomial, det(lambda*eye(size(A))-A).
% 
% Contributed by Paulo Xavier Candeias under GPL.
function p = poly(r)
   if (nargin < 1) | (nargout > 1)
      error('wrong use (see help poly)')
   end
   
   % Identify arguments
   if issquare(r)
      r = eig(r);
   else
      r = r(:).';
   end
   
   % Strip out infinities
   r = r( ~isinf(r) );
   
   % Expand recursion formula
   n = length(r);
   p = [1,zeros(1,n)];
   for v = 1:n
      p(2:v+1) = p(2:v+1)-r(v)*p(1:v);
   end

   % The result should be real if the roots are complex conjugates.
   r_neg_comp_roots = r(imag(r)<0);
   r_pos_comp_roots = r(imag(r)>0);
   if (numel(r_neg_comp_roots) == numel(r_pos_comp_roots))
     if sort(r(imag(r)>0))-sort(conj(r(imag(r)<0)))<teps(r)
        p = real(p);
     end
   end