synfig-core/trunk/src/synfig/polynomial_root.cpp

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