X-Git-Url: https://git.pterodactylus.net/?a=blobdiff_plain;f=synfig-core%2Ftags%2Fstable%2Fsrc%2Fsynfig%2Fpolynomial_root.cpp;fp=synfig-core%2Ftags%2Fstable%2Fsrc%2Fsynfig%2Fpolynomial_root.cpp;h=ab792f8899d548bff9300a8d0bbe86e5b9369e70;hb=684cf1db661a9a5cbc142238cf05d6d2f7aa3f89;hp=0000000000000000000000000000000000000000;hpb=829dade3901a0f6e606075811b8b5c0d1bb5393e;p=synfig.git diff --git a/synfig-core/tags/stable/src/synfig/polynomial_root.cpp b/synfig-core/tags/stable/src/synfig/polynomial_root.cpp new file mode 100644 index 0000000..ab792f8 --- /dev/null +++ b/synfig-core/tags/stable/src/synfig/polynomial_root.cpp @@ -0,0 +1,267 @@ +/* === 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 +#endif + +#include "polynomial_root.h" +#include + +#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 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 + + //Efficent 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; + } +}