1 /*
2 * $Id: zroots.i,v 1.1 2005-09-18 22:06:13 dhmunro Exp $
3 * Laguerre's method for finding roots of complex polynomial.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func zroots(a, imsort=)
12 /* DOCUMENT zroots(a)
13 returns the vector of (complex) roots of the (complex)
14 polynomial: Sum(a(i)*x^(i-1)) using Laguerre's method.
15 The roots are sorted in order of increasing real parts. Call with
16 non-zero imsort keyword to sort into increasing imaginary parts.
17 See Numerical Recipes (Press, et. al., Cambridge Univ. Press 1988),
18 section 9.5.
19 */
20 {
21 m= numberof(a)-1; /* degree of poly */
22 EPS= 1.e-6;
23 EPS*= 2.*EPS;
24 is_complex= (structof(a)==complex);
25 if (!is_complex) a= a+0i;
26 ad= a; /* copy of coeffs for successive deflation */
27 roots= array(0i, m);
28
29 for (j=m ; j>=1 ; j--) {
30 /* Loop over each root to be found */
31 x= laguerre(ad,0i); /* start at zero to favor smallest rem root */
32 if (abs(x.im)<=EPS*abs(x.re)) x= x.re+0i;
33 roots(j)= x;
34 b= ad(j+1); /* forward deflation */
35 for (jj=j ; jj>=1 ; jj--) {
36 c= ad(jj);
37 ad(jj)= b;
38 b= x*b+c;
39 }
40 ad= ad(1:j); /* poly has one lower degree */
41 }
42
43 /* polish */
44 for (j=1 ; j<=m ; j++) roots(j)= laguerre(a,roots(j));
45
46 /* sort roots */
47 if (numberof(roots) == 1) return roots;
48 re= roots.re;
49 im= roots.im;
50 if (!is_complex &&
51 allof(abs(im)<=EPS*abs(re))) return re(sort(re));
52
53 /* sorting on multiple keys is difficult because all fast sorting
54 algorithms (such as the quicksort used by sort) randomize the
55 order of equal elements in the list to be sorted -- the following
56 is a modification of the general deterministic sorting algorithm
57 implemented in msort.i, modified to ensure that conjugate roots
58 are recognized even when their real parts are not quite equal */
59 if (imsort) {
60 rank= re;
61 re= im;
62 im= rank;
63 }
64 list= rank= sort(re);
65 re= re(list);
66 rank(list)= (re(dif)>EPS*abs(re(zcen)))(cum); /* see msort.i */
67 list= rerank= sort(im);
68 rerank(list)= indgen(m);
69 return roots(sort(rank+double(rerank)/m));
70 }
71
laguerre(a,x)72 func laguerre(a,x)
73 /* DOCUMENT laguerre(a,x)
74 Given the coefficients a(1:m+1) of the m'th degree
75 complex polynomial Sum(a(i)*x^(i-1)) and a guess x, returns a root.
76 See Numerical Recipes (Press, et. al., Cambridge Univ. Press 1988),
77 section 9.5.
78 */
79 {
80 EPSS=3.e-15;
81 MR=8;
82 MT=10;
83 MAXIT=MT*MR;
84 frac=[.5,.25,.75,.13,.38,.62,.88,1.]; /* Fractions to break limit cycles */
85
86 m= numberof(a)-1; /* degree of poly */
87 a= a+0i; /* make complex */
88 if (m<2) return -a(1)/a(2);
89
90 for (iter=1 ; iter<=MAXIT ; iter++) {
91 b= a(0); /* coefft of max degree */
92 d= f= 0i;
93 err= abs(b);
94 abx=abs(x);
95 for (j=m ; j>=1 ; j--) { /* compute poly and 2 derivs */
96 f= f*x+d;
97 d= d*x+b;
98 b= b*x+a(j);
99 err= abs(b)+abx*err;
100 }
101 err*= EPSS; /* estimated round-off error in poly */
102 if (abs(b)<=err) return x; /* we're on the root already */
103
104 g= d/b;
105 g2= g*g;
106 h= g2-2.*f/b;
107 sq= sqrt((m-1)*(m*h-g2));
108 gp= g+sq;
109 gm= g-sq;
110 abp= abs(gp);
111 abm= abs(gm);
112 if (abp<abm) gp= gm;
113
114 if (gp) dx= m/gp;
115 else dx= exp(log(1+abx)+1.i*iter);
116
117 x1= x-dx;
118 if (abs(x-x1)<EPSS) return x1; /* converged to machine accuracy */
119
120 if (iter%MT) x= x1; /* next iteration */
121 else x= x-dx*frac(iter/MT); /* break limit cycle? */
122 }
123 write, "WARNING: too many iterations in laguerre";
124 return x;
125 }
126
127
128
129