X-Git-Url: https://git.pterodactylus.net/?a=blobdiff_plain;f=synfig-core%2Ftags%2Fsynfig_0_61_05%2Fsynfig-core%2Fsrc%2Fsynfig%2Fpolynomial_root.cpp;fp=synfig-core%2Ftags%2Fsynfig_0_61_05%2Fsynfig-core%2Fsrc%2Fsynfig%2Fpolynomial_root.cpp;h=0000000000000000000000000000000000000000;hb=6fa8f2f38d4b0b35f8539bf94e27ae27015c7689;hp=ea2962b68ed55457fc1000720b422121ce24c2ec;hpb=47fce282611fbba1044921d22ca887f9b53ad91a;p=synfig.git diff --git a/synfig-core/tags/synfig_0_61_05/synfig-core/src/synfig/polynomial_root.cpp b/synfig-core/tags/synfig_0_61_05/synfig-core/src/synfig/polynomial_root.cpp deleted file mode 100644 index ea2962b..0000000 --- a/synfig-core/tags/synfig_0_61_05/synfig-core/src/synfig/polynomial_root.cpp +++ /dev/null @@ -1,267 +0,0 @@ -/* === S Y N F I G ========================================================= */ -/*! \file template.cpp -** \brief Template File -** -** $Id: polynomial_root.cpp,v 1.1.1.1 2005/01/04 01:23:14 darco Exp $ -** -** \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 acheivable roundoff limit, to a root of teh 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 it's 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 ad 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 teh 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; - } -}