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

   1 /* === S Y N F I G ========================================================= */
   2 /*!     \file polynomial_root.h
   3 **      \brief Polynomial Root Finder Header
   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 /* === S T A R T =========================================================== */
  24
  25 #ifndef __SYNFIG_POLYNOMIAL_ROOT_H
  26 #define __SYNFIG_POLYNOMIAL_ROOT_H
  27
  28 /* === H E A D E R S ======================================================= */
  29
  30 #include <complex>
  31 #include <vector>
  32
  33 /* === M A C R O S ========================================================= */
  34
  35 /* === T Y P E D E F S ===================================================== */
  36
  37 /* === C L A S S E S & S T R U C T S ======================================= */
  38 template < typename T = float, typename F = float >
  39 class Polynomial : public std::vector<T> //a0 + a1x + a2x^2 + ... + anx^n
  40 {
  41 public:
  42
  43         //Will maintain all lower constants
  44         void degree(unsigned int d, const T & def = (T)0) { resize(d+1,def); }
  45         unsigned int degree()const { return size() - 1; }
  46
  47         const Polynomial & operator+=(const Polynomial &p)
  48         {
  49                 if(p.size() > size())
  50                         resize(p.size(), (T)0);
  51
  52                 for(int i = 0; i < p.size(); ++i)
  53                 {
  54                         (*this)[i] += p[i];
  55                 }
  56                 return *this;
  57         }
  58
  59         const Polynomial & operator-=(const Polynomial &p)
  60         {
  61                 if(p.size() > size())
  62                         resize(p.size(), (T)0);
  63
  64                 for(int i = 0; i < p.size(); ++i)
  65                 {
  66                         (*this)[i] -= p[i];
  67                 }
  68                 return *this;
  69         }
  70
  71         const Polynomial & operator*=(const Polynomial &p)
  72         {
  73                 if(p.size() < 1)
  74                 {
  75                         resize(0);
  76                         return *this;
  77                 }
  78
  79                 unsigned int i,j;
  80                 std::vector<T> nc(*this);
  81
  82                 //in place for constant stuff
  83                 for(i = 0; i < nc.size(); ++i)
  84                 {
  85                         (*this)[i] *= p[0];
  86                 }
  87
  88                 if(p.size() < 2) return *this;
  89
  90                 resize(size() + p.degree());
  91                 for(int i = 0; i < nc.size(); ++i)
  92                 {
  93                         for(int j = 1; j < p.size(); ++j)
  94                         {
  95                                 nc[i+j] += nc[i]*p[j];
  96                         }
  97                 }
  98
  99                 return *this;
 100         }
 101 };
 102
 103 class RootFinder
 104 {
 105         std::vector< std::complex<float> >      workcoefs;
 106         int     its;
 107
 108 public:
 109         std::vector< std::complex<float> >      coefs; //the number of coefficients determines the degree of polynomial
 110
 111         std::vector< std::complex<float> >      roots;
 112
 113         void find_all_roots(bool polish);
 114 };
 115
 116
 117
 118 /* === E N D =============================================================== */
 119
 120 #endif