1 /* === S Y N F I G ========================================================= */
2 /*! \file polynomial_root.cpp
3 ** \brief Template File
4 **
5 ** $Id$
6 **
7 ** \legal
8 ** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
9 **
10 ** This package is free software; you can redistribute it and/or
11 ** modify it under the terms of the GNU General Public License as
12 ** published by the Free Software Foundation; either version 2 of
13 ** the License, or (at your option) any later version.
14 **
15 ** This package is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 ** General Public License for more details.
19 ** \endlegal
20 */
21 /* ========================================================================= */
22
23 /* === H E A D E R S ======================================================= */
24
25 #ifdef USING_PCH
26 # include "pch.h"
27 #else
28 #ifdef HAVE_CONFIG_H
29 # include <config.h>
30 #endif
31
32 #include "polynomial_root.h"
33 #include <complex>
34
35 #endif
36
37 /* === U S I N G =========================================================== */
38
39 using namespace std;
40 //using namespace etl;
41 //using namespace synfig;
42
43 /* === M A C R O S ========================================================= */
44
45 /* === G L O B A L S ======================================================= */
46 typedef complex<float> Complex;
47
48 /* === P R O C E D U R E S ================================================= */
49
50 /* === M E T H O D S ======================================================= */
51
52 #define EPSS 1.0e-7
53 #define MR 8
54 #define MT 10
55 #define MAXIT (MT*MR)
56
57 /*EPSS is the estimated fractional roundoff error. We try to break (rare) limit
58 cycles with MR different fractional values, once every MT steps, for MAXIT total allowed iterations.
59 */
60
61 /* Explanation:
62
63 A polynomial can be represented like so:
64 Pn(x) = (x - x1)(x - x2)...(x - xn) where xi = complex roots
65
66 We can get the following:
67 ln|Pn(x)| = ln|x - x1| + ln|x - x2| + ... + ln|x - xn|
68
69 G := d ln|Pn(x)| / dx =
70 +1/(x-x1) + 1/(x-x2) + ... + 1/(x-xn)
71
72 and
73
74 H := - d2 ln|Pn(x)| / d2x =
75 +1/(x-x1)^2 + 1/(x-x2)^2 + ... + 1/(x-xn)^2
76
77 which gives
78 H = [Pn'/Pn]^2 - Pn''/Pn
79
80 Laguerre's formula guesses that the root we are seeking x1 is located
81 some distance a from our current guess x, and all the other roots are
82 located at distance b.
83
84 Using this:
85
86 1/a + (n-1)/b = G
87
88 and
89
90 1/a^2 + (n-1)/b^2 = H
91
92 which yields this solution for a:
93
94 a = n / G +- sqrt( (n-1)(nH - G^2) )
95
96 where +- is determined by which ever yields the largest magnitude for the denominator.
97 a can easily be complex since the factor inside the square-root can be negative.
98
99 This method iterates (x=x-a) until a is sufficiently small.
100 */
101
102 /* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial sum(i=0,m){a[i]x^i},
103 and given a complex value x, this routine improves x by laguerre's method until it converges,
104 within the achievable roundoff limit, to a root of the given polynomial. The number of iterations taken
105 is returned as `its'.
106 */
laguer(Complex a[],int m,Complex * x,int * its)107 void laguer(Complex a[], int m, Complex *x, int *its)
108 {
109 int iter,j;
110 float abx, abp, abm, err;
111 Complex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
112
113 //Fractions used to break a limit cycle
114 static float frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
115
116 for(iter = 1; iter <= MAXIT; ++iter)
117 {
118 *its = iter; //number of iterations so far
119
120 b = a[m]; //the highest coefficient
121 err = abs(b); //its magnitude
122
123 d = f = Complex(0,0); //clear variables for use
124 abx = abs(*x); //the magnitude of the current root
125
126 //Efficient computation of the polynomial and its first 2 derivatives
127 for(j = m-1; j >= 0; --j)
128 {
129 f = (*x)*f + d;
130 d = (*x)*d + b;
131 b = (*x)*b + a[j];
132
133 err = abs(b) + abx*err;
134 }
135
136 //Estimate the roundoff error in evaluation polynomial
137 err *= EPSS;
138
139 //Are we on the root?
140 if(abs(b) < err)
141 {
142 return;
143 }
144
145 //General case: use Laguerre's formula
146 //a = n / G +- sqrt( (n-1)(nH - G^2) )
147 //x = x - a
148
149 g = d / b; //get G
150 g2 = g * g; //for the sqrt calc
151
152 h = g2 - 2.0f * (f / b); //get H
153
154 sq = pow( (float)(m-1) * ((float)m*h - g2), 0.5f ); //get the sqrt
155
156 //get the denominator
157 gp = g + sq;
158 gm = g - sq;
159
160 abp = abs(gp);
161 abm = abs(gm);
162
163 //get the denominator with the highest magnitude
164 if(abp < abm)
165 {
166 abp = abm;
167 gp = gm;
168 }
169
170 //if the denominator is positive do one thing, otherwise do the other
171 dx = (abp > 0.0) ? (float)m / gp : polar((1+abx),(float)iter);
172 x1 = *x - dx;
173
174 //Have we converged?
175 if( *x == x1 )
176 {
177 return;
178 }
179
180 //Every so often take a fractional step, to break any limit cycle (itself a rare occurrence).
181 if( iter % MT )
182 {
183 *x = x1;
184 }else
185 {
186 *x = *x - (frac[iter/MT]*dx);
187 }
188 }
189
190 //very unusual - can occur only for complex roots. Try a different starting guess for the root.
191 //nrerror("too many iterations in laguer");
192 return;
193 }
194
195 #define EPS 2.0e-6
196 #define MAXM 100 //a small number, and maximum anticipated value of m..
197
198 /* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial a0 + a1*x +...+ an*x^n
199 the routine successively calls laguer and finds all m complex roots in roots[1..m].
200 The boolean variable polish should be input as true (1) if polishing (also by Laguerre's Method)
201 is desired, false (0) if the roots will be subsequently polished by other means.
202 */
find_all_roots(bool polish)203 void RootFinder::find_all_roots(bool polish)
204 {
205 int i,its,j,jj;
206 Complex x,b,c;
207 int m = coefs.size()-1;
208
209 //make sure roots is big enough
210 roots.resize(m);
211
212 if(workcoefs.size() < MAXM) workcoefs.resize(MAXM);
213
214 //Copy the coefficients for successive deflation
215 for(j = 0; j <= m; ++j)
216 {
217 workcoefs[j] = coefs[j];
218 }
219
220 //Loop over each root to be found
221 for(j = m-1; j >= 0; --j)
222 {
223 //Start at 0 to favor convergence to smallest remaining root, and find the root
224 x = Complex(0,0);
225 laguer(&workcoefs[0],j+1,&x,&its); //must add 1 to get the degree
226
227 //if it is close enough to a real root, then make it so
228 if(abs(x.imag()) <= 2.0*EPS*abs(x.real()))
229 {
230 x = Complex(x.real());
231 }
232
233 roots[j] = x;
234
235 //forward deflation
236
237 //the degree is j+1 since j(0,m-1)
238 b = workcoefs[j+1];
239 for(jj = j; jj >= 0; --jj)
240 {
241 c = workcoefs[jj];
242 workcoefs[jj] = b;
243 b = x*b + c;
244 }
245 }
246
247 //Polish the roots using the undeflated coefficients
248 if(polish)
249 {
250 for(j = 0; j < m; ++j)
251 {
252 laguer(&coefs[0],m,&roots[j],&its);
253 }
254 }
255
256 //Sort roots by their real parts by straight insertion
257 for(j = 1; j < m; ++j)
258 {
259 x = roots[j];
260 for( i = j-1; i >= 1; --i)
261 {
262 if(roots[i].real() <= x.real()) break;
263 roots[i+1] = roots[i];
264 }
265 roots[i+1] = x;
266 }
267 }
268