1 /* === S Y N F I G ========================================================= */
2 /*! \file polynomial_root.cpp
3 ** \brief Template File
8 ** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
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.
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.
21 /* ========================================================================= */
23 /* === H E A D E R S ======================================================= */
32 #include "polynomial_root.h"
37 /* === U S I N G =========================================================== */
40 //using namespace etl;
41 //using namespace synfig;
43 /* === M A C R O S ========================================================= */
45 /* === G L O B A L S ======================================================= */
46 typedef complex<float> Complex;
48 /* === P R O C E D U R E S ================================================= */
50 /* === M E T H O D S ======================================================= */
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.
63 A polynomial can be represented like so:
64 Pn(x) = (x - x1)(x - x2)...(x - xn) where xi = complex roots
66 We can get the following:
67 ln|Pn(x)| = ln|x - x1| + ln|x - x2| + ... + ln|x - xn|
69 G := d ln|Pn(x)| / dx =
70 +1/(x-x1) + 1/(x-x2) + ... + 1/(x-xn)
74 H := - d2 ln|Pn(x)| / d2x =
75 +1/(x-x1)^2 + 1/(x-x2)^2 + ... + 1/(x-xn)^2
78 H = [Pn'/Pn]^2 - Pn''/Pn
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.
92 which yields this solution for a:
94 a = n / G +- sqrt( (n-1)(nH - G^2) )
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.
99 This method iterates (x=x-a) until a is sufficiently small.
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'.
107 void laguer(Complex a[], int m, Complex *x, int *its)
110 float abx, abp, abm, err;
111 Complex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
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};
116 for(iter = 1; iter <= MAXIT; ++iter)
118 *its = iter; //number of iterations so far
120 b = a[m]; //the highest coefficient
121 err = abs(b); //its magnitude
123 d = f = Complex(0,0); //clear variables for use
124 abx = abs(*x); //the magnitude of the current root
126 //Efficent computation of the polynomial and its first 2 derivatives
127 for(j = m-1; j >= 0; --j)
133 err = abs(b) + abx*err;
136 //Estimate the roundoff error in evaluation polynomial
139 //Are we on the root?
145 //General case: use Laguerre's formula
146 //a = n / G +- sqrt( (n-1)(nH - G^2) )
150 g2 = g * g; //for the sqrt calc
152 h = g2 - 2.0f * (f / b); //get H
154 sq = pow( (float)(m-1) * ((float)m*h - g2), 0.5f ); //get the sqrt
156 //get the denominator
163 //get the denominator with the highest magnitude
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);
180 //Every so often take a fractional step, to break any limit cycle (itself a rare occurrence).
186 *x = *x - (frac[iter/MT]*dx);
190 //very unusual - can occur only for complex roots. Try a different starting guess for the root.
191 //nrerror("too many iterations in laguer");
196 #define MAXM 100 //a small number, and maximum anticipated value of m..
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.
203 void RootFinder::find_all_roots(bool polish)
207 int m = coefs.size()-1;
209 //make sure roots is big enough
212 if(workcoefs.size() < MAXM) workcoefs.resize(MAXM);
214 //Copy the coefficients for successive deflation
215 for(j = 0; j <= m; ++j)
217 workcoefs[j] = coefs[j];
220 //Loop over each root to be found
221 for(j = m-1; j >= 0; --j)
223 //Start at 0 to favor convergence to smallest remaining root, and find the root
225 laguer(&workcoefs[0],j+1,&x,&its); //must add 1 to get the degree
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()))
230 x = Complex(x.real());
237 //the degree is j+1 since j(0,m-1)
239 for(jj = j; jj >= 0; --jj)
247 //Polish the roots using the undeflated coefficients
250 for(j = 0; j < m; ++j)
252 laguer(&coefs[0],m,&roots[j],&its);
256 //Sort roots by their real parts by straight insertion
257 for(j = 1; j < m; ++j)
260 for( i = j-1; i >= 1; --i)
262 if(roots[i].real() <= x.real()) break;
263 roots[i+1] = roots[i];