synfig-core/tags/synfig_0_61_07_rc1/src/synfig/polynomial_root.cpp

   1 /* === S Y N F I G ========================================================= */
   2 /*!     \file polynomial_root.cpp
   3 **      \brief Template File
   4 **
   5 **      $Id$
   6 **
   7 **      \legal
   8 **      Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
   9 **
  10 **      This package is free software; you can redistribute it and/or
  11 **      modify it under the terms of the GNU General Public License as
  12 **      published by the Free Software Foundation; either version 2 of
  13 **      the License, or (at your option) any later version.
  14 **
  15 **      This package is distributed in the hope that it will be useful,
  16 **      but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 **      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  18 **      General Public License for more details.
  19 **      \endlegal
  20 */
  21 /* ========================================================================= */
  22
  23 /* === H E A D E R S ======================================================= */
  24
  25 #ifdef USING_PCH
  26 #       include "pch.h"
  27 #else
  28 #ifdef HAVE_CONFIG_H
  29 #       include <config.h>
  30 #endif
  31
  32 #include "polynomial_root.h"
  33 #include <complex>
  34
  35 #endif
  36
  37 /* === U S I N G =========================================================== */
  38
  39 using namespace std;
  40 //using namespace etl;
  41 //using namespace synfig;
  42
  43 /* === M A C R O S ========================================================= */
  44
  45 /* === G L O B A L S ======================================================= */
  46 typedef complex<float>  Complex;
  47
  48 /* === P R O C E D U R E S ================================================= */
  49
  50 /* === M E T H O D S ======================================================= */
  51
  52 #define EPSS 1.0e-7
  53 #define MR 8
  54 #define MT 10
  55 #define MAXIT (MT*MR)
  56
  57 /*EPSS is the estimated fractional roundoff error.  We try to break (rare) limit
  58 cycles with MR different fractional values, once every MT steps, for MAXIT total allowed iterations.
  59 */
  60
  61 /*      Explanation:
  62
  63         A polynomial can be represented like so:
  64         Pn(x) = (x - x1)(x - x2)...(x - xn)     where xi = complex roots
  65
  66         We can get the following:
  67         ln|Pn(x)| = ln|x - x1| + ln|x - x2| + ... + ln|x - xn|
  68
  69         G := d ln|Pn(x)| / dx =
  70         +1/(x-x1) + 1/(x-x2) + ... + 1/(x-xn)
  71
  72         and
  73
  74         H := - d2 ln|Pn(x)| / d2x =
  75         +1/(x-x1)^2 + 1/(x-x2)^2 + ... + 1/(x-xn)^2
  76
  77         which gives
  78         H = [Pn'/Pn]^2 - Pn''/Pn
  79
  80         Laguerre's formula guesses that the root we are seeking x1 is located
  81         some distance a from our current guess x, and all the other roots are
  82         located at distance b.
  83
  84         Using this:
  85
  86         1/a + (n-1)/b = G
  87
  88         and
  89
  90         1/a^2 + (n-1)/b^2 = H
  91
  92         which yields this solution for a:
  93
  94         a = n / G +- sqrt( (n-1)(nH - G^2) )
  95
  96         where +- is determined by which ever yields the largest magnitude for the denominator.
  97         a can easily be complex since the factor inside the square-root can be negative.
  98
  99         This method iterates (x=x-a) until a is sufficiently small.
 100 */
 101
 102 /* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial sum(i=0,m){a[i]x^i},
 103 and given a complex value x, this routine improves x by laguerre's method until it converges,
 104 within the achievable roundoff limit, to a root of the given polynomial.  The number of iterations taken
 105 is returned as `its'.
 106 */
 107 void laguer(Complex a[], int m, Complex *x, int *its)
 108 {
 109         int iter,j;
 110         float abx, abp, abm, err;
 111         Complex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
 112
 113         //Fractions used to break a limit cycle
 114         static float frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
 115
 116         for(iter = 1; iter <= MAXIT; ++iter)
 117         {
 118                 *its = iter; //number of iterations so far
 119
 120                 b       = a[m];         //the highest coefficient
 121                 err     = abs(b);       //its magnitude
 122
 123                 d = f = Complex(0,0); //clear variables for use
 124                 abx = abs(*x);  //the magnitude of the current root
 125
 126                 //Efficent computation of the polynomial and its first 2 derivatives
 127                 for(j = m-1; j >= 0; --j)
 128                 {
 129                         f = (*x)*f + d;
 130                         d = (*x)*d + b;
 131                         b = (*x)*b + a[j];
 132
 133                         err = abs(b) + abx*err;
 134                 }
 135
 136                 //Estimate the roundoff error in evaluation polynomial
 137                 err *= EPSS;
 138
 139                 //Are we on the root?
 140                 if(abs(b) < err)
 141                 {
 142                         return;
 143                 }
 144
 145                 //General case: use Laguerre's formula
 146                 //a = n / G +- sqrt( (n-1)(nH - G^2) )
 147                 //x = x - a
 148
 149                 g = d / b;      //get G
 150                 g2 = g * g; //for the sqrt calc
 151
 152                 h = g2 - 2.0f * (f / b);        //get H
 153
 154                 sq = pow( (float)(m-1) * ((float)m*h - g2), 0.5f ); //get the sqrt
 155
 156                 //get the denominator
 157                 gp = g + sq;
 158                 gm = g - sq;
 159
 160                 abp = abs(gp);
 161                 abm = abs(gm);
 162
 163                 //get the denominator with the highest magnitude
 164                 if(abp < abm)
 165                 {
 166                         abp = abm;
 167                         gp = gm;
 168                 }
 169
 170                 //if the denominator is positive do one thing, otherwise do the other
 171                 dx = (abp > 0.0) ? (float)m / gp : polar((1+abx),(float)iter);
 172                 x1 = *x - dx;
 173
 174                 //Have we converged?
 175                 if( *x == x1 )
 176                 {
 177                         return;
 178                 }
 179
 180                 //Every so often take a fractional step, to break any limit cycle (itself a rare occurrence).
 181                 if( iter % MT )
 182                 {
 183                         *x = x1;
 184                 }else
 185                 {
 186                         *x = *x - (frac[iter/MT]*dx);
 187                 }
 188         }
 189
 190         //very unusual - can occur only for complex roots.  Try a different starting guess for the root.
 191         //nrerror("too many iterations in laguer");
 192         return;
 193 }
 194
 195 #define EPS 2.0e-6
 196 #define MAXM 100        //a small number, and maximum anticipated value of m..
 197
 198 /*      Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial a0 + a1*x +...+ an*x^n
 199         the routine successively calls laguer and finds all m complex roots in roots[1..m].
 200         The boolean variable polish should be input as true (1) if polishing (also by Laguerre's Method)
 201         is desired, false (0) if the roots will be subsequently polished by other means.
 202 */
 203 void RootFinder::find_all_roots(bool polish)
 204 {
 205         int i,its,j,jj;
 206         Complex x,b,c;
 207         int m = coefs.size()-1;
 208
 209         //make sure roots is big enough
 210         roots.resize(m);
 211
 212         if(workcoefs.size() < MAXM) workcoefs.resize(MAXM);
 213
 214         //Copy the coefficients for successive deflation
 215         for(j = 0; j <= m; ++j)
 216         {
 217                 workcoefs[j] = coefs[j];
 218         }
 219
 220         //Loop over each root to be found
 221         for(j = m-1; j >= 0; --j)
 222         {
 223                 //Start at 0 to favor convergence to smallest remaining root, and find the root
 224                 x = Complex(0,0);
 225                 laguer(&workcoefs[0],j+1,&x,&its); //must add 1 to get the degree
 226
 227                 //if it is close enough to a real root, then make it so
 228                 if(abs(x.imag()) <= 2.0*EPS*abs(x.real()))
 229                 {
 230                         x = Complex(x.real());
 231                 }
 232
 233                 roots[j] = x;
 234
 235                 //forward deflation
 236
 237                 //the degree is j+1 since j(0,m-1)
 238                 b = workcoefs[j+1];
 239                 for(jj = j; jj >= 0; --jj)
 240                 {
 241                         c = workcoefs[jj];
 242                         workcoefs[jj] = b;
 243                         b = x*b + c;
 244                 }
 245         }
 246
 247         //Polish the roots using the undeflated coefficients
 248         if(polish)
 249         {
 250                 for(j = 0; j < m; ++j)
 251                 {
 252                         laguer(&coefs[0],m,&roots[j],&its);
 253                 }
 254         }
 255
 256         //Sort roots by their real parts by straight insertion
 257         for(j = 1; j < m; ++j)
 258         {
 259                 x = roots[j];
 260                 for( i = j-1; i >= 1; --i)
 261                 {
 262                         if(roots[i].real() <= x.real()) break;
 263                         roots[i+1] = roots[i];
 264                 }
 265                 roots[i+1] = x;
 266         }
 267 }