Add an interpolation sub-library

[libs/math.git] / source / interpolate / polynomial.h

diff --git a/source/interpolate/polynomial.h b/source/interpolate/polynomial.h
new file mode 100644 (file)
index 0000000..859eea8
--- /dev/null
+++ b/source/interpolate/polynomial.h
@@ -0,0 +1,168 @@
+#ifndef MSP_INTERPOLATE_POLYNOMIAL_H_
+#define MSP_INTERPOLATE_POLYNOMIAL_H_
+
+#include <algorithm>
+#include <msp/linal/vector.h>
+
+namespace Msp {
+namespace Interpolate {
+
+/**
+Represents a polynomial of degree D with a single variable.  The coefficients
+are stored as a vector with the highest degree term first.
+*/
+template<typename T, unsigned D>
+class Polynomial
+{
+private:
+       LinAl::Vector<T, D+1> coeff;
+
+public:
+       Polynomial() { }
+       Polynomial(const LinAl::Vector<T, D+1> &c): coeff(c) { }
+       template<typename U, unsigned G>
+       Polynomial(const Polynomial<U, G> &);
+
+       unsigned degree() const { return D; }
+
+       LinAl::Vector<T, D+1> coefficients() { return coeff; }
+       const LinAl::Vector<T, D+1> coefficients() const { return coeff; }
+       T &coefficient(unsigned i) { return coeff[i]; }
+       const T &coefficient(unsigned i) const { return coeff[i]; }
+
+       Polynomial &operator+=(const Polynomial &);
+       Polynomial &operator-=(const Polynomial &);
+       Polynomial &operator*=(T);
+       Polynomial &operator/=(T);
+
+       /// Calculates the value of the polynomial at a particular point.
+       T operator()(T) const;
+
+       /** Substitutes the variable of the polynomial with another polynomial.  The
+       result is a new polynomial. */
+       template<unsigned G>
+       Polynomial<T, D*G> operator()(const Polynomial<T, G> &) const;
+};
+
+template<typename T, unsigned D>
+template<typename U, unsigned G>
+inline Polynomial<T, D>::Polynomial(const Polynomial<U, G> &other)
+{
+       char source_cannot_be_larger[D>=G ? 1 : -1];
+       for(unsigned i=0; i<=G; ++i)
+               coeff[D-i] = other.coefficient(G-i);
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator+=(const Polynomial<T, D> &p)
+{
+       coeff += p.coeff;
+       return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator+(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
+{
+       Polynomial<T, D> r(p);
+       return r += q;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator-=(const Polynomial<T, D> &p)
+{
+       coeff -= p.coeff;
+       return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator-(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
+{
+       Polynomial<T, D> r(p);
+       return r -= q;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator*=(T s)
+{
+       coeff *= s;
+       return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator*(const Polynomial<T, D> &p, T s)
+{
+       Polynomial<T, D> r(p);
+       return r *= s;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator*(T s, const Polynomial<T, D> &p)
+{
+       return p*s;
+}
+
+template<typename T, unsigned D, unsigned G>
+inline Polynomial<T, D+G> operator*(const Polynomial<T, D> &p, const Polynomial<T, G> &q)
+{
+       LinAl::Vector<T, D+G+1> coeff;
+       for(unsigned i=0; i<=D; ++i)
+               for(unsigned j=0; j<=G; ++j)
+                       coeff[i+j] += p.coefficient(i)*q.coefficient(j);
+       return Polynomial<T, D+G>(coeff);
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator/=(T s)
+{
+       coeff /= s;
+       return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator/(const Polynomial<T, D> &p, T s)
+{
+       Polynomial<T, D> r(p);
+       return r /= s;
+}
+
+template<typename T, unsigned D>
+inline T Polynomial<T, D>::operator()(T x) const
+{
+       T r = T();
+       T x_i = T(1);
+       for(int i=D; i>=0; --i, x_i*=x)
+               r += coeff[i]*x_i;
+       return r;
+}
+
+template<typename T, unsigned D>
+template<unsigned G>
+inline Polynomial<T, D*G> Polynomial<T, D>::operator()(const Polynomial<T, G> &p) const
+{
+       Polynomial<T, D*G> r;
+       r.coefficient(D*G) = coeff[D];
+
+       LinAl::Vector<T, D*G+1> cf;
+       cf[0] = T(1);
+
+       for(unsigned i=1; i<=D; ++i)
+       {
+               for(int j=i*G; j>=0; --j)
+               {
+                       T c = T();
+                       for(unsigned k=0; (k<=static_cast<unsigned>(j) && k<=G); ++k)
+                               c += cf[j-k]*p.coefficient(G-k);
+                       cf[j] = c;
+               }
+
+               for(unsigned j=0; j<=i*G; ++j)
+                       r.coefficient(D*G-j) += coeff[D-i]*cf[j];
+       }
+
+       return r;
+}
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif