1 #ifndef MSP_INTERPOLATE_POLYNOMIAL_H_
2 #define MSP_INTERPOLATE_POLYNOMIAL_H_
5 #include <msp/linal/vector.h>
8 namespace Interpolate {
11 Represents a polynomial of degree D with a single variable. The coefficients
12 are stored as a vector with the highest degree term first.
14 template<typename T, unsigned D>
18 LinAl::Vector<T, D+1> coeff;
22 Polynomial(const LinAl::Vector<T, D+1> &c): coeff(c) { }
23 template<typename U, unsigned G>
24 Polynomial(const Polynomial<U, G> &);
26 unsigned degree() const { return D; }
28 LinAl::Vector<T, D+1> coefficients() { return coeff; }
29 const LinAl::Vector<T, D+1> coefficients() const { return coeff; }
30 T &coefficient(unsigned i) { return coeff[i]; }
31 const T &coefficient(unsigned i) const { return coeff[i]; }
33 Polynomial &operator+=(const Polynomial &);
34 Polynomial &operator-=(const Polynomial &);
35 Polynomial &operator*=(T);
36 Polynomial &operator/=(T);
38 /// Calculates the value of the polynomial at a particular point.
39 T operator()(T) const;
41 /** Substitutes the variable of the polynomial with another polynomial. The
42 result is a new polynomial. */
44 Polynomial<T, D*G> operator()(const Polynomial<T, G> &) const;
47 template<typename T, unsigned D>
48 template<typename U, unsigned G>
49 inline Polynomial<T, D>::Polynomial(const Polynomial<U, G> &other)
51 char source_cannot_be_larger[D>=G ? 1 : -1];
52 for(unsigned i=0; i<=G; ++i)
53 coeff[D-i] = other.coefficient(G-i);
56 template<typename T, unsigned D>
57 inline Polynomial<T, D> &Polynomial<T, D>::operator+=(const Polynomial<T, D> &p)
63 template<typename T, unsigned D>
64 inline Polynomial<T, D> operator+(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
66 Polynomial<T, D> r(p);
70 template<typename T, unsigned D>
71 inline Polynomial<T, D> &Polynomial<T, D>::operator-=(const Polynomial<T, D> &p)
77 template<typename T, unsigned D>
78 inline Polynomial<T, D> operator-(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
80 Polynomial<T, D> r(p);
84 template<typename T, unsigned D>
85 inline Polynomial<T, D> &Polynomial<T, D>::operator*=(T s)
91 template<typename T, unsigned D>
92 inline Polynomial<T, D> operator*(const Polynomial<T, D> &p, T s)
94 Polynomial<T, D> r(p);
98 template<typename T, unsigned D>
99 inline Polynomial<T, D> operator*(T s, const Polynomial<T, D> &p)
104 template<typename T, unsigned D, unsigned G>
105 inline Polynomial<T, D+G> operator*(const Polynomial<T, D> &p, const Polynomial<T, G> &q)
107 LinAl::Vector<T, D+G+1> coeff;
108 for(unsigned i=0; i<=D; ++i)
109 for(unsigned j=0; j<=G; ++j)
110 coeff[i+j] += p.coefficient(i)*q.coefficient(j);
111 return Polynomial<T, D+G>(coeff);
114 template<typename T, unsigned D>
115 inline Polynomial<T, D> &Polynomial<T, D>::operator/=(T s)
121 template<typename T, unsigned D>
122 inline Polynomial<T, D> operator/(const Polynomial<T, D> &p, T s)
124 Polynomial<T, D> r(p);
128 template<typename T, unsigned D>
129 inline T Polynomial<T, D>::operator()(T x) const
133 for(int i=D; i>=0; --i, x_i*=x)
138 template<typename T, unsigned D>
140 inline Polynomial<T, D*G> Polynomial<T, D>::operator()(const Polynomial<T, G> &p) const
142 Polynomial<T, D*G> r;
143 r.coefficient(D*G) = coeff[D];
145 LinAl::Vector<T, D*G+1> cf;
148 for(unsigned i=1; i<=D; ++i)
150 for(int j=i*G; j>=0; --j)
153 for(unsigned k=0; (k<=static_cast<unsigned>(j) && k<=G); ++k)
154 c += cf[j-k]*p.coefficient(G-k);
158 for(unsigned j=0; j<=i*G; ++j)
159 r.coefficient(D*G-j) += coeff[D-i]*cf[j];
165 } // namespace Interpolate