00001
00021
00022
00023
00024
00025 #ifndef __SYNFIG_POLYNOMIAL_ROOT_H
00026 #define __SYNFIG_POLYNOMIAL_ROOT_H
00027
00028
00029
00030 #include <complex>
00031 #include <vector>
00032
00033
00034
00035
00036
00037
00038 template < typename T = float, typename F = float >
00039 class Polynomial : public std::vector<T>
00040 {
00041 public:
00042
00043
00044 void degree(unsigned int d, const T & def = (T)0) { resize(d+1,def); }
00045 unsigned int degree()const { return size() - 1; }
00046
00047 const Polynomial & operator+=(const Polynomial &p)
00048 {
00049 if(p.size() > size())
00050 resize(p.size(), (T)0);
00051
00052 for(int i = 0; i < p.size(); ++i)
00053 {
00054 (*this)[i] += p[i];
00055 }
00056 return *this;
00057 }
00058
00059 const Polynomial & operator-=(const Polynomial &p)
00060 {
00061 if(p.size() > size())
00062 resize(p.size(), (T)0);
00063
00064 for(int i = 0; i < p.size(); ++i)
00065 {
00066 (*this)[i] -= p[i];
00067 }
00068 return *this;
00069 }
00070
00071 const Polynomial & operator*=(const Polynomial &p)
00072 {
00073 if(p.size() < 1)
00074 {
00075 resize(0);
00076 return *this;
00077 }
00078
00079 unsigned int i,j;
00080 std::vector<T> nc(*this);
00081
00082
00083 for(i = 0; i < nc.size(); ++i)
00084 {
00085 (*this)[i] *= p[0];
00086 }
00087
00088 if(p.size() < 2) return *this;
00089
00090 resize(size() + p.degree());
00091 for(int i = 0; i < nc.size(); ++i)
00092 {
00093 for(int j = 1; j < p.size(); ++j)
00094 {
00095 nc[i+j] += nc[i]*p[j];
00096 }
00097 }
00098
00099 return *this;
00100 }
00101 };
00102
00103 class RootFinder
00104 {
00105 std::vector< std::complex<float> > workcoefs;
00106 int its;
00107
00108 public:
00109 std::vector< std::complex<float> > coefs;
00110
00111 std::vector< std::complex<float> > roots;
00112
00113 void find_all_roots(bool polish);
00114 };
00115
00116
00117
00118
00119
00120 #endif