synfig-core/tags/0.61.09/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 **      Copyright (c) 2007 Chris Moore
  10 **
  11 **      This package is free software; you can redistribute it and/or
  12 **      modify it under the terms of the GNU General Public License as
  13 **      published by the Free Software Foundation; either version 2 of
  14 **      the License, or (at your option) any later version.
  15 **
  16 **      This package is distributed in the hope that it will be useful,
  17 **      but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 **      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  19 **      General Public License for more details.
  20 **      \endlegal
  21 */
  22 /* ========================================================================= */
  23
  24 /* === S T A R T =========================================================== */
  25
  26 #ifndef __SYNFIG_POLYNOMIAL_ROOT_H
  27 #define __SYNFIG_POLYNOMIAL_ROOT_H
  28
  29 /* === H E A D E R S ======================================================= */
  30
  31 #include <complex>
  32 #include <vector>
  33
  34 /* === M A C R O S ========================================================= */
  35
  36 /* === T Y P E D E F S ===================================================== */
  37
  38 /* === C L A S S E S & S T R U C T S ======================================= */
  39 template < typename T = float, typename F = float >
  40 class Polynomial : public std::vector<T> //a0 + a1x + a2x^2 + ... + anx^n
  41 {
  42 public:
  43
  44         //Will maintain all lower constants
  45         void degree(unsigned int d, const T & def = (T)0) { resize(d+1,def); }
  46         unsigned int degree()const { return this->size() - 1; }
  47
  48         const Polynomial & operator+=(const Polynomial &p)
  49         {
  50                 if(p.size() > this->size())
  51                         resize(p.size(), (T)0);
  52
  53                 for(int i = 0; i < p.size(); ++i)
  54                 {
  55                         (*this)[i] += p[i];
  56                 }
  57                 return *this;
  58         }
  59
  60         const Polynomial & operator-=(const Polynomial &p)
  61         {
  62                 if(p.size() > this->size())
  63                         resize(p.size(), (T)0);
  64
  65                 for(int i = 0; i < p.size(); ++i)
  66                 {
  67                         (*this)[i] -= p[i];
  68                 }
  69                 return *this;
  70         }
  71
  72         const Polynomial & operator*=(const Polynomial &p)
  73         {
  74                 if(p.size() < 1)
  75                 {
  76                         this->resize(0);
  77                         return *this;
  78                 }
  79
  80                 unsigned int i,j;
  81                 std::vector<T> nc(*this);
  82
  83                 //in place for constant stuff
  84                 for(i = 0; i < nc.size(); ++i)
  85                 {
  86                         (*this)[i] *= p[0];
  87                 }
  88
  89                 if(p.size() < 2) return *this;
  90
  91                 this->resize(this->size() + p.degree());
  92                 for(int i = 0; i < nc.size(); ++i)
  93                 {
  94                         for(int j = 1; j < p.size(); ++j)
  95                         {
  96                                 nc[i+j] += nc[i]*p[j];
  97                         }
  98                 }
  99
 100                 return *this;
 101         }
 102 };
 103
 104 class RootFinder
 105 {
 106         std::vector< std::complex<float> >      workcoefs;
 107         int     its;
 108
 109 public:
 110         std::vector< std::complex<float> >      coefs; //the number of coefficients determines the degree of polynomial
 111
 112         std::vector< std::complex<float> >      roots;
 113
 114         void find_all_roots(bool polish);
 115 };
 116
 117
 118
 119 /* === E N D =============================================================== */
 120
 121 #endif