/* === S Y N F I G ========================================================= */
-/*! \file template.cpp
+/*! \file polynomial_root.cpp
** \brief Template File
**
-** $Id: polynomial_root.cpp,v 1.1.1.1 2005/01/04 01:23:14 darco Exp $
+** $Id$
**
** \legal
-** Copyright (c) 2002 Robert B. Quattlebaum Jr.
+** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
**
-** This software and associated documentation
-** are CONFIDENTIAL and PROPRIETARY property of
-** the above-mentioned copyright holder.
+** 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.
**
-** You may not copy, print, publish, or in any
-** other way distribute this software without
-** a prior written agreement with
-** the copyright holder.
+** 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
*/
/* ========================================================================= */
*/
/* 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|
-
+ 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
+
+ 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
-
+
+ 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.
+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 it's first 2 derivatives
+
+ //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);
-
+ 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 = *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
+/* 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 teh roots will be subsequently polished by other means.
+ 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
{
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);
+ 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)
b = x*b + c;
}
}
-
+
//Polish the roots using the undeflated coefficients
if(polish)
{
laguer(&coefs[0],m,&roots[j],&its);
}
}
-
+
//Sort roots by their real parts by straight insertion
for(j = 1; j < m; ++j)
{
projects / synfig.git / blobdiff