--- /dev/null
+/* === S Y N F I G ========================================================= */
+/*! \file polynomial_root.cpp
+** \brief Template File
+**
+** $Id$
+**
+** \legal
+** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
+**
+** This package is free software; you can redistribute it and/or
+** modify it under the terms of the GNU General Public License as
+** published by the Free Software Foundation; either version 2 of
+** the License, or (at your option) any later version.
+**
+** This package is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** General Public License for more details.
+** \endlegal
+*/
+/* ========================================================================= */
+
+/* === H E A D E R S ======================================================= */
+
+#ifdef USING_PCH
+# include "pch.h"
+#else
+#ifdef HAVE_CONFIG_H
+# include <config.h>
+#endif
+
+#include "polynomial_root.h"
+#include <complex>
+
+#endif
+
+/* === U S I N G =========================================================== */
+
+using namespace std;
+//using namespace etl;
+//using namespace synfig;
+
+/* === M A C R O S ========================================================= */
+
+/* === G L O B A L S ======================================================= */
+typedef complex<float> Complex;
+
+/* === P R O C E D U R E S ================================================= */
+
+/* === M E T H O D S ======================================================= */
+
+#define EPSS 1.0e-7
+#define MR 8
+#define MT 10
+#define MAXIT (MT*MR)
+
+/*EPSS is the estimated fractional roundoff error. We try to break (rare) limit
+cycles with MR different fractional values, once every MT steps, for MAXIT total allowed iterations.
+*/
+
+/* Explanation:
+
+ A polynomial can be represented like so:
+ Pn(x) = (x - x1)(x - x2)...(x - xn) where xi = complex roots
+
+ We can get the following:
+ ln|Pn(x)| = ln|x - x1| + ln|x - x2| + ... + ln|x - xn|
+
+ G := d ln|Pn(x)| / dx =
+ +1/(x-x1) + 1/(x-x2) + ... + 1/(x-xn)
+
+ and
+
+ H := - d2 ln|Pn(x)| / d2x =
+ +1/(x-x1)^2 + 1/(x-x2)^2 + ... + 1/(x-xn)^2
+
+ which gives
+ H = [Pn'/Pn]^2 - Pn''/Pn
+
+ Laguerre's formula guesses that the root we are seeking x1 is located
+ some distance a from our current guess x, and all the other roots are
+ located at distance b.
+
+ Using this:
+
+ 1/a + (n-1)/b = G
+
+ and
+
+ 1/a^2 + (n-1)/b^2 = H
+
+ which yields this solution for a:
+
+ a = n / G +- sqrt( (n-1)(nH - G^2) )
+
+ where +- is determined by which ever yields the largest magnitude for the denominator.
+ a can easily be complex since the factor inside the square-root can be negative.
+
+ This method iterates (x=x-a) until a is sufficiently small.
+*/
+
+/* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial sum(i=0,m){a[i]x^i},
+and given a complex value x, this routine improves x by laguerre's method until it converges,
+within the achievable roundoff limit, to a root of the given polynomial. The number of iterations taken
+is returned as `its'.
+*/
+void laguer(Complex a[], int m, Complex *x, int *its)
+{
+ int iter,j;
+ float abx, abp, abm, err;
+ Complex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
+
+ //Fractions used to break a limit cycle
+ static float frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
+
+ for(iter = 1; iter <= MAXIT; ++iter)
+ {
+ *its = iter; //number of iterations so far
+
+ b = a[m]; //the highest coefficient
+ err = abs(b); //its magnitude
+
+ d = f = Complex(0,0); //clear variables for use
+ abx = abs(*x); //the magnitude of the current root
+
+ //Efficient computation of the polynomial and its first 2 derivatives
+ for(j = m-1; j >= 0; --j)
+ {
+ f = (*x)*f + d;
+ d = (*x)*d + b;
+ b = (*x)*b + a[j];
+
+ err = abs(b) + abx*err;
+ }
+
+ //Estimate the roundoff error in evaluation polynomial
+ err *= EPSS;
+
+ //Are we on the root?
+ if(abs(b) < err)
+ {
+ return;
+ }
+
+ //General case: use Laguerre's formula
+ //a = n / G +- sqrt( (n-1)(nH - G^2) )
+ //x = x - a
+
+ g = d / b; //get G
+ g2 = g * g; //for the sqrt calc
+
+ h = g2 - 2.0f * (f / b); //get H
+
+ sq = pow( (float)(m-1) * ((float)m*h - g2), 0.5f ); //get the sqrt
+
+ //get the denominator
+ gp = g + sq;
+ gm = g - sq;
+
+ abp = abs(gp);
+ abm = abs(gm);
+
+ //get the denominator with the highest magnitude
+ if(abp < abm)
+ {
+ abp = abm;
+ gp = gm;
+ }
+
+ //if the denominator is positive do one thing, otherwise do the other
+ dx = (abp > 0.0) ? (float)m / gp : polar((1+abx),(float)iter);
+ x1 = *x - dx;
+
+ //Have we converged?
+ if( *x == x1 )
+ {
+ return;
+ }
+
+ //Every so often take a fractional step, to break any limit cycle (itself a rare occurrence).
+ if( iter % MT )
+ {
+ *x = x1;
+ }else
+ {
+ *x = *x - (frac[iter/MT]*dx);
+ }
+ }
+
+ //very unusual - can occur only for complex roots. Try a different starting guess for the root.
+ //nrerror("too many iterations in laguer");
+ return;
+}
+
+#define EPS 2.0e-6
+#define MAXM 100 //a small number, and maximum anticipated value of m..
+
+/* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial a0 + a1*x +...+ an*x^n
+ the routine successively calls laguer and finds all m complex roots in roots[1..m].
+ The boolean variable polish should be input as true (1) if polishing (also by Laguerre's Method)
+ is desired, false (0) if the roots will be subsequently polished by other means.
+*/
+void RootFinder::find_all_roots(bool polish)
+{
+ int i,its,j,jj;
+ Complex x,b,c;
+ int m = coefs.size()-1;
+
+ //make sure roots is big enough
+ roots.resize(m);
+
+ if(workcoefs.size() < MAXM) workcoefs.resize(MAXM);
+
+ //Copy the coefficients for successive deflation
+ for(j = 0; j <= m; ++j)
+ {
+ workcoefs[j] = coefs[j];
+ }
+
+ //Loop over each root to be found
+ for(j = m-1; j >= 0; --j)
+ {
+ //Start at 0 to favor convergence to smallest remaining root, and find the root
+ x = Complex(0,0);
+ laguer(&workcoefs[0],j+1,&x,&its); //must add 1 to get the degree
+
+ //if it is close enough to a real root, then make it so
+ if(abs(x.imag()) <= 2.0*EPS*abs(x.real()))
+ {
+ x = Complex(x.real());
+ }
+
+ roots[j] = x;
+
+ //forward deflation
+
+ //the degree is j+1 since j(0,m-1)
+ b = workcoefs[j+1];
+ for(jj = j; jj >= 0; --jj)
+ {
+ c = workcoefs[jj];
+ workcoefs[jj] = b;
+ b = x*b + c;
+ }
+ }
+
+ //Polish the roots using the undeflated coefficients
+ if(polish)
+ {
+ for(j = 0; j < m; ++j)
+ {
+ laguer(&coefs[0],m,&roots[j],&its);
+ }
+ }
+
+ //Sort roots by their real parts by straight insertion
+ for(j = 1; j < m; ++j)
+ {
+ x = roots[j];
+ for( i = j-1; i >= 1; --i)
+ {
+ if(roots[i].real() <= x.real()) break;
+ roots[i+1] = roots[i];
+ }
+ roots[i+1] = x;
+ }
+}
projects / synfig.git / blobdiff