source/interpolate/polynomial.h

   1 #ifndef MSP_INTERPOLATE_POLYNOMIAL_H_
   2 #define MSP_INTERPOLATE_POLYNOMIAL_H_
   3
   4 #include <algorithm>
   5 #include <msp/linal/vector.h>
   6
   7 namespace Msp {
   8 namespace Interpolate {
   9
  10 /**
  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.
  13 */
  14 template<typename T, unsigned D>
  15 class Polynomial
  16 {
  17 private:
  18         LinAl::Vector<T, D+1> coeff;
  19
  20 public:
  21         Polynomial() { }
  22         Polynomial(const LinAl::Vector<T, D+1> &c): coeff(c) { }
  23         template<typename U, unsigned G>
  24         Polynomial(const Polynomial<U, G> &);
  25
  26         unsigned degree() const { return D; }
  27
  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]; }
  32
  33         Polynomial &operator+=(const Polynomial &);
  34         Polynomial &operator-=(const Polynomial &);
  35         Polynomial &operator*=(T);
  36         Polynomial &operator/=(T);
  37
  38         /// Calculates the value of the polynomial at a particular point.
  39         T operator()(T) const;
  40
  41         /** Substitutes the variable of the polynomial with another polynomial.  The
  42         result is a new polynomial. */
  43         template<unsigned G>
  44         Polynomial<T, D*G> operator()(const Polynomial<T, G> &) const;
  45 };
  46
  47 template<typename T, unsigned D>
  48 template<typename U, unsigned G>
  49 inline Polynomial<T, D>::Polynomial(const Polynomial<U, G> &other)
  50 {
  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);
  54 }
  55
  56 template<typename T, unsigned D>
  57 inline Polynomial<T, D> &Polynomial<T, D>::operator+=(const Polynomial<T, D> &p)
  58 {
  59         coeff += p.coeff;
  60         return *this;
  61 }
  62
  63 template<typename T, unsigned D>
  64 inline Polynomial<T, D> operator+(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
  65 {
  66         Polynomial<T, D> r(p);
  67         return r += q;
  68 }
  69
  70 template<typename T, unsigned D>
  71 inline Polynomial<T, D> &Polynomial<T, D>::operator-=(const Polynomial<T, D> &p)
  72 {
  73         coeff -= p.coeff;
  74         return *this;
  75 }
  76
  77 template<typename T, unsigned D>
  78 inline Polynomial<T, D> operator-(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
  79 {
  80         Polynomial<T, D> r(p);
  81         return r -= q;
  82 }
  83
  84 template<typename T, unsigned D>
  85 inline Polynomial<T, D> &Polynomial<T, D>::operator*=(T s)
  86 {
  87         coeff *= s;
  88         return *this;
  89 }
  90
  91 template<typename T, unsigned D>
  92 inline Polynomial<T, D> operator*(const Polynomial<T, D> &p, T s)
  93 {
  94         Polynomial<T, D> r(p);
  95         return r *= s;
  96 }
  97
  98 template<typename T, unsigned D>
  99 inline Polynomial<T, D> operator*(T s, const Polynomial<T, D> &p)
 100 {
 101         return p*s;
 102 }
 103
 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)
 106 {
 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);
 112 }
 113
 114 template<typename T, unsigned D>
 115 inline Polynomial<T, D> &Polynomial<T, D>::operator/=(T s)
 116 {
 117         coeff /= s;
 118         return *this;
 119 }
 120
 121 template<typename T, unsigned D>
 122 inline Polynomial<T, D> operator/(const Polynomial<T, D> &p, T s)
 123 {
 124         Polynomial<T, D> r(p);
 125         return r /= s;
 126 }
 127
 128 template<typename T, unsigned D>
 129 inline T Polynomial<T, D>::operator()(T x) const
 130 {
 131         T r = T();
 132         T x_i = T(1);
 133         for(int i=D; i>=0; --i, x_i*=x)
 134                 r += coeff[i]*x_i;
 135         return r;
 136 }
 137
 138 template<typename T, unsigned D>
 139 template<unsigned G>
 140 inline Polynomial<T, D*G> Polynomial<T, D>::operator()(const Polynomial<T, G> &p) const
 141 {
 142         Polynomial<T, D*G> r;
 143         r.coefficient(D*G) = coeff[D];
 144
 145         LinAl::Vector<T, D*G+1> cf;
 146         cf[0] = T(1);
 147
 148         for(unsigned i=1; i<=D; ++i)
 149         {
 150                 for(int j=i*G; j>=0; --j)
 151                 {
 152                         T c = T();
 153                         for(unsigned k=0; (k<=static_cast<unsigned>(j) && k<=G); ++k)
 154                                 c += cf[j-k]*p.coefficient(G-k);
 155                         cf[j] = c;
 156                 }
 157
 158                 for(unsigned j=0; j<=i*G; ++j)
 159                         r.coefficient(D*G-j) += coeff[D-i]*cf[j];
 160         }
 161
 162         return r;
 163 }
 164
 165 } // namespace Interpolate
 166 } // namespace Msp
 167
 168 #endif