+++ /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