[Boost-commit] svn:boost r48261 - sandbox/SOC/2008/polynomial/boost

From: lord.zerom1nd_at_[hidden]
Date: 2008-08-20 09:26:49

Author: kieliszczyk
Date: 2008-08-20 09:26:49 EDT (Wed, 20 Aug 2008)
New Revision: 48261
URL: http://svn.boost.org/trac/boost/changeset/48261

Log:
poly
Added:
   sandbox/SOC/2008/polynomial/boost/polynomial.hpp   (contents, props changed)

Added: sandbox/SOC/2008/polynomial/boost/polynomial.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2008/polynomial/boost/polynomial.hpp	2008-08-20 09:26:49 EDT (Wed, 20 Aug 2008)
@@ -0,0 +1,678 @@
+#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP
+#define BOOST_MATH_TOOLS_POLYNOMIAL_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <algorithm>
+#include <complex>
+#include <functional>
+#include <iosfwd>
+#include <vector>
+#include <boost/assert.hpp>
+#include <boost/math/common_factor.hpp>
+#include <boost/numeric/conversion/cast.hpp>
+#include <boost/math/special_functions/round.hpp>
+#include "polynomial/fast_fourier_transform.hpp"
+#include "polynomial/polynomial_special_forms.hpp"
+#include "polynomial/polynomial_special_functions.hpp"
+
+namespace boost {
+namespace math {
+namespace tools {
+
+  template<typename FieldType>
+  class polynomial;
+  template<typename IntegerType>
+  polynomial<IntegerType> gcd(polynomial<IntegerType> u, polynomial<IntegerType> v);
+
+  namespace detail {
+
+    template<typename InputIterator>
+    inline typename std::iterator_traits<InputIterator>::value_type cont_evaluate(InputIterator first, InputIterator last) {
+      BOOST_ASSERT(first != last);
+      typename std::iterator_traits<InputIterator>::value_type cont = *first++;
+      while(first != last)
+        cont = boost::math::gcd(cont, *first++);
+      return cont;
+    }
+
+    template<typename InputIterator1, typename InputIterator2, typename OutputIterator>
+    void polynomial_add_in_place(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2, OutputIterator result) {
+      while(first1 != last1 && first2 != last2)
+        *result++ = (*first1++) + (*first2++);
+      while(first1 != last1)
+        *result++ = *first1++;
+      while(first2 != last2)
+        *result++ = *first2++;
+    }
+
+    template<typename InputIterator1, typename InputIterator2, typename OutputIterator>
+    void polynomial_sub_in_place(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2, OutputIterator result) {
+      while(first1 != last1 && first2 != last2)
+        *result++ = (*first1++) - (*first2++);
+      while(first1 != last1)
+        *result++ = *first1++;
+      while(first2 != last2)
+        *result++ = -(*first2++);
+    }
+
+    template<typename IntegerType>
+    polynomial<IntegerType> modulo_for_gcd(std::vector<IntegerType> u, const std::vector<IntegerType>& v) {
+      typedef typename polynomial<IntegerType>::size_type size_type;
+      if(u.size() < v.size())
+        return u;
+      if(v.size() == 1)
+        return polynomial<IntegerType>();
+      size_type k = u.size() - v.size();
+      size_type n = v.size() - 1;
+      do {
+        size_type j = n + k - 1;
+        do {
+          if(j < k)
+            u[j] = v[n] * u[j];
+          else
+            u[j] = v[n] * u[j] - u[n + k] * v[j - k];
+        }
+        while(j--);
+      }
+      while(k--);
+      while(u.size() >= v.size())
+        u.pop_back();
+      while(u.size() > 1 && u.back() == boost::numeric_cast<IntegerType>(0))
+        u.pop_back();
+      return polynomial<IntegerType>(u);
+    }
+
+  } // namespace detail
+
+  template<typename FieldType>
+  class polynomial {
+    std::vector<FieldType> coefficients;
+    std::vector<double> preconditions;
+    template<typename U>
+    U round(const U& v);
+    void normalize();
+    template<typename U>
+    static FieldType cvt(U const& n) {
+      return boost::numeric_cast<FieldType>(n);
+    }
+    static FieldType zero() {
+      return cvt(0.0);
+    }
+    friend polynomial<FieldType> boost::math::tools::gcd<>(polynomial<FieldType> u, polynomial<FieldType> v);
+  public:
+    // typedefs:
+    typedef typename std::vector<FieldType>::value_type value_type;
+    typedef typename std::vector<FieldType>::size_type size_type;
+
+    // constructors:
+    polynomial() {
+      coefficients.push_back(zero());
+    }
+    template<typename InputIterator>
+    polynomial(InputIterator first, InputIterator last) : coefficients(first, last) {
+      normalize();
+    }
+    template<typename U>
+    polynomial(const U* coef, const size_type n) : coefficients(coef, coef + n) {
+      normalize();
+    }
+    template<typename U>
+    polynomial(const U& v) {
+      coefficients.push_back(v);
+      normalize();
+    }
+    template<typename InputIterator>
+    polynomial(InputIterator first1, InputIterator last1, InputIterator first2) : coefficients(std::distance(first1, last1)) {
+      boost::math::tools::interpolate_polynomial(first1, last1, first2, coefficients.begin());
+      normalize();
+    }
+    polynomial(std::vector<FieldType>& c) {
+      coefficients.swap(c);
+      normalize();
+    }
+
+    // copy constructor:
+    template<typename U>
+    polynomial(const polynomial<U>& poly) {
+      std::transform(poly.begin(), poly.end(), std::back_inserter(coefficients), std::ptr_fun<U, FieldType>(boost::numeric_cast));
+      normalize();
+    }
+
+    // modification:
+    template<typename U>
+    polynomial<FieldType>& operator=(const polynomial<U>& poly);
+    polynomial<FieldType>& operator=(std::vector<FieldType>& c) {
+      coefficients.swap(c);
+      normalize();
+      return *this;
+    }
+    template<typename InputIterator>
+    void interpolate(InputIterator first1, InputIterator last1, InputIterator first2) {
+      coefficients.resize(std::distance(first1, last1));
+      boost::math::tools::interpolate_polynomial(first1, last1, first2, coefficients.begin());
+      normalize();
+    }
+    void swap(polynomial<FieldType>& poly) {
+      coefficients.swap(poly.coefficients);
+    }
+
+    // access:
+    size_type size() const {
+      return coefficients.size();
+    }
+    size_type degree() const {
+      return coefficients.size() - 1;
+    }
+    const value_type& operator[](const size_type i) const {
+      return coefficients[i];
+    }
+    value_type& operator[](const size_type i) {
+      return coefficients[i];
+    }
+    typename std::vector<FieldType>::const_iterator begin() const {
+      return coefficients.begin();
+    }
+    typename std::vector<FieldType>::iterator begin() {
+      return coefficients.begin();
+    }
+    typename std::vector<FieldType>::const_iterator end() const {
+      return coefficients.end();
+    }
+    typename std::vector<FieldType>::iterator end() {
+      return coefficients.end();
+    }
+
+    // evaluation:
+    FieldType operator()(const FieldType& z) const {
+      return boost::math::tools::evaluate_polynomial(coefficients.rbegin(), coefficients.rend(), z);
+    }
+    FieldType evaluate(const FieldType& z) const {
+      return boost::math::tools::evaluate_polynomial(coefficients.rbegin(), coefficients.rend(), z);
+    }
+    FieldType evaluate_faithfully(const FieldType& z) const {
+      return boost::math::tools::evaluate_polynomial_faithfully(coefficients.rbegin(), coefficients.rend(), z);
+    }
+    void precondition_coefficients() {
+      BOOST_ASSERT(size() > 4 && size() < 8);
+      preconditions.resize(size());
+      boost::math::tools::precondition_polynomial_coefficients(begin(), end(), preconditions.begin());
+    }
+    FieldType evaluate_by_preconditioning(const FieldType& z) const {
+      return boost::math::tools::evaluate_polynomial_by_preconditioning(preconditions.begin(), preconditions.end(), z);
+    }
+    // derivatives and integrals:
+    polynomial<FieldType> derivative() const;
+    polynomial<FieldType> integral(const FieldType& c) const;
+    FieldType definite_integral(const FieldType& a, const FieldType& b) const;
+
+    // operators:
+
+    const polynomial<FieldType>& operator+() const;
+    polynomial<FieldType> operator-() const;
+
+    // addition:
+    template<typename U>
+    polynomial<FieldType>& operator+=(const polynomial<U>& poly);
+    template<typename U>
+    polynomial<FieldType>& operator+=(const U& v);
+
+    // substraction:
+    template<typename U>
+    polynomial<FieldType>& operator-=(const polynomial<U>& poly);
+    template<typename U>
+    polynomial<FieldType>& operator-=(const U& v);
+
+    // multiplication:
+    template<typename U>
+    polynomial<FieldType>& operator*=(const polynomial<U>& poly);
+    template<typename U>
+    polynomial<FieldType>& operator*=(const U& v);
+
+    // division:
+    template<typename U>
+    polynomial<FieldType>& operator/=(const polynomial<U>& poly);
+    template<typename U>
+    polynomial<FieldType>& operator/=(const U& v);
+
+    // modulo:
+    template<typename U>
+    polynomial<FieldType>& operator%=(const polynomial<U>& poly);
+  };
+
+  template<typename FieldType>
+  template<typename U>
+  inline U polynomial<FieldType>::round(const U& v) {
+    return boost::math::round(v);
+  }
+
+  // specializations for floating-point types (rounding shouldn't be done):
+  template<>
+  template<typename U>
+  inline U polynomial<float>::round(const U& v) {
+    return v;
+  }
+
+  template<>
+  template<typename U>
+  inline U polynomial<double>::round(const U& v) {
+    return v;
+  }
+
+  template<>
+  template<typename U>
+  inline U polynomial<long double>::round(const U& v) {
+    return v;
+  }
+
+
+  template<typename FieldType>
+  inline void polynomial<FieldType>::normalize() {
+    while(coefficients.size() > 1 && coefficients.back() == zero())
+      coefficients.pop_back();
+    if(coefficients.empty())
+      coefficients.push_back(zero());
+  }
+
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator=(const polynomial<U>& poly) {
+    if(reinterpret_cast<void const*>(this) != reinterpret_cast<void const*>(&poly)) {
+      std::vector<FieldType> tmp(poly.size());
+      typename std::vector<FieldType>::iterator tmp_it = tmp.begin();
+      typename std::vector<U>::const_iterator poly_it = poly.begin();
+      while(tmp_it != tmp.end())
+        *tmp_it++ = cvt(*poly_it++);
+      coefficients.swap(tmp);
+      normalize();
+    }
+    return *this;
+  }
+
+
+  // derivatives and integrals:
+  template<typename FieldType>
+  inline polynomial<FieldType> polynomial<FieldType>::derivative() const {
+    std::vector<FieldType> tmp(degree());
+    boost::math::tools::copy_polynomial_derivative(coefficients.begin(), coefficients.end(), tmp.begin());
+    return polynomial<FieldType>(tmp.begin(), tmp.end());
+  }
+
+  template<typename FieldType>
+  inline polynomial<FieldType> polynomial<FieldType>::integral(const FieldType& c = FieldType(0.0)) const {
+    std::vector<FieldType> tmp(size() + 1);
+    boost::math::tools::copy_polynomial_integral(coefficients.begin(), coefficients.end(), tmp.begin(), c);
+    return polynomial<FieldType>(tmp.begin(), tmp.end());
+  }
+
+  template<typename FieldType>
+  inline FieldType polynomial<FieldType>::definite_integral(const FieldType& a, const FieldType& b) const {
+    return boost::math::tools::evaluate_polynomial_definite_integral(coefficients.begin(), coefficients.end(), a, b);
+  }
+
+
+  // addition:
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator+=(const polynomial<U>& poly) {
+    std::vector<FieldType> tmp(std::max(size(), poly.size()));
+    boost::math::tools::detail::polynomial_add_in_place(begin(), end(), poly.begin(), poly.end(), tmp.begin());
+    std::swap(coefficients, tmp);
+    normalize();
+    return *this;
+  }
+
+  template<typename FieldType>
+  template<typename U>
+  inline polynomial<FieldType>& polynomial<FieldType>::operator+=(const U& v) {
+    coefficients.front() += cvt(v);
+    return *this;
+  }
+
+  // substraction:
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator-=(const polynomial<U>& poly) {
+    if(reinterpret_cast<void const*>(this) != reinterpret_cast<void const*>(&poly)) {
+      std::vector<FieldType> tmp(std::max(size(), poly.size()));
+      boost::math::tools::detail::polynomial_sub_in_place(begin(), end(), poly.begin(), poly.end(), tmp.begin());
+      std::swap(coefficients, tmp);
+      normalize();
+    }
+    else {
+      coefficients.clear();
+      coefficients.push_back(zero());
+    }
+    return *this;
+  }
+
+  template<typename FieldType>
+  template<typename U>
+  inline polynomial<FieldType>& polynomial<FieldType>::operator-=(const U& v) {
+    coefficients.front() -= cvt(v);
+    return *this;
+  }
+
+  // multiplication:
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator*=(const polynomial<U>& poly) {
+    double l = log2(std::max(coefficients.size(), poly.size()) * 2);
+    size_type new_sz = static_cast<size_type>(floor(l));
+    if(l > new_sz)
+      new_sz = static_cast<size_type>(pow(2.0, static_cast<int>(new_sz + 1)));
+    else
+      new_sz = static_cast<size_type>(pow(2.0, static_cast<int>(new_sz)));
+    std::vector<std::complex<double> > dft1, dft2;
+    dft1.reserve(new_sz);
+    dft2.reserve(new_sz);
+    std::transform(begin(), end(), std::back_inserter(dft1), std::ptr_fun<double, FieldType>(boost::numeric_cast));
+    std::fill_n(std::back_inserter(dft1), new_sz - dft1.size(), zero());
+    std::transform(poly.begin(), poly.end(), std::back_inserter(dft2), std::ptr_fun<double, FieldType>(boost::numeric_cast));
+    std::fill_n(std::back_inserter(dft2), new_sz - dft2.size(), zero());
+    boost::math::tools::fast_fourier_transform(dft1.begin(), dft1.end(), dft1.begin());
+    boost::math::tools::fast_fourier_transform(dft2.begin(), dft2.end(), dft2.begin());
+    std::transform(dft1.begin(), dft1.end(), dft2.begin(), dft1.begin(), std::multiplies<std::complex<double> >());
+    boost::math::tools::inverse_fast_fourier_transform(dft1.begin(), dft1.end(), dft1.begin());
+    coefficients.resize(coefficients.size() + poly.size() - 1);
+    typename std::vector<std::complex<double> >::iterator dft1_it = dft1.begin();
+    typename std::vector<FieldType>::iterator c_it = begin();
+    // probably std::transform might be used here but I can't make it work properly:
+    while(c_it != end())
+      *c_it++ = cvt(this->round((dft1_it++)->real()));
+    normalize();
+    return *this;
+  }
+
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator*=(const U& v) {
+    std::vector<FieldType> tmp;
+    for(typename std::vector<FieldType>::iterator c_it = begin(); c_it != end(); ++c_it)
+      tmp.push_back(cvt(*c_it * v));
+    std::swap(coefficients, tmp);
+    normalize();
+    return *this;
+  }
+
+  // division:
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator/=(const polynomial<U>& poly) {
+    if(reinterpret_cast<void const*>(this) != reinterpret_cast<void const*>(&poly)) {
+      std::vector<FieldType> tmp(coefficients);
+      if(size() < poly.size()) {
+        tmp.clear();
+        tmp.push_back(zero());
+        std::swap(coefficients, tmp);
+        return *this;
+      }
+      size_type k = tmp.size() - poly.size();
+      std::vector<double> new_coefficients;
+      new_coefficients.resize(k + 1);
+      while(true) {
+        new_coefficients[k] = boost::numeric_cast<double>(tmp[poly.size() + k - 1]) / poly.coefficients.back();
+        for(size_type j = poly.size() + k - 1; j > k; --j)
+          tmp[j] = cvt(tmp[j] - new_coefficients[k] * poly[j - k]);
+        tmp[k] = cvt(tmp[k] - new_coefficients[k] * poly[0]);
+        if(k-- == 0)
+          break;
+      }
+      tmp.resize(new_coefficients.size());
+      std::transform(new_coefficients.begin(), new_coefficients.end(), tmp.begin(), std::ptr_fun<double, FieldType>(boost::numeric_cast));
+      std::swap(coefficients, tmp);
+      normalize();
+    }
+    else {
+      coefficients.clear();
+      coefficients.push_back(cvt(1.0));
+    }
+    return *this;
+  }
+
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator/=(const U& v) {
+    std::vector<FieldType> tmp;
+    for(typename std::vector<FieldType>::iterator c_it = begin(); c_it != end(); ++c_it)
+      tmp.push_back(cvt(*c_it / v));
+    std::swap(coefficients, tmp);
+    normalize();
+    return *this;
+  }
+
+  // modulo:
+  template<typename FieldType>
+  template<typename U>
+  polynomial<FieldType>& polynomial<FieldType>::operator%=(const polynomial<U>& poly) {
+    if(reinterpret_cast<void const*>(this) != reinterpret_cast<void const*>(&poly)) {
+      if(size() < poly.size()) {
+        return *this;
+      }
+      std::vector<FieldType> tmp(coefficients);
+      size_type k = tmp.size() - poly.size();
+      double q;
+      while(true) {
+        q = boost::numeric_cast<double>(tmp[poly.size() + k - 1]) / poly.coefficients.back();
+        for(size_type j = poly.size() + k - 1; j > k; --j)
+          tmp[j] = cvt(tmp[j] - q * poly[j - k]);
+        tmp[k] = cvt(tmp[k] - q * poly[0]);
+        if(k-- == 0)
+          break;
+      }
+      while(tmp.size() >= poly.size())
+        tmp.pop_back();
+      std::swap(coefficients, tmp);
+      normalize();
+    }
+    else {
+      coefficients.clear();
+      coefficients.push_back(zero());
+    }
+    return *this;
+  }
+
+  // greatest common divisor:
+  template<typename IntegerType>
+  polynomial<IntegerType> gcd(polynomial<IntegerType> u, polynomial<IntegerType> v) {
+    IntegerType cont_u = boost::math::tools::detail::cont_evaluate(u.begin(), u.end());
+    u /= cont_u;
+    IntegerType cont_v = boost::math::tools::detail::cont_evaluate(v.begin(), v.end());
+    v /= cont_v;
+    while(true) {
+      polynomial<IntegerType> r = boost::math::tools::detail::modulo_for_gcd(u.coefficients, v.coefficients);
+      if(r.size() == 1) {
+        if(*(r.begin()) == r.zero())
+          break;
+        v = polynomial<IntegerType>(boost::numeric_cast<IntegerType>(1.0));
+        break;
+      }
+      u = v;
+      v = r / boost::math::tools::detail::cont_evaluate(r.begin(), r.end());
+    }
+    if(v.coefficients.back() < v.zero())
+      return -boost::math::gcd(cont_u, cont_v) * v;
+    return boost::math::gcd(cont_u, cont_v) * v;
+  }
+
+
+  // operators:
+
+  template<typename FieldType>
+  inline const polynomial<FieldType>& polynomial<FieldType>::operator+() const {
+    return *this;
+  }
+
+  template<typename FieldType>
+  inline polynomial<FieldType> polynomial<FieldType>::operator-() const {
+    std::vector<FieldType> tmp;
+    tmp.reserve(size());
+    std::transform(begin(), end(), std::back_inserter(tmp), std::negate<FieldType>());
+    return polynomial<FieldType>(tmp);
+  }
+
+  // addition:
+  template<typename FieldType>
+  inline polynomial<FieldType> operator+(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    polynomial<FieldType> result(left);
+    result += right;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator+(const polynomial<FieldType>& poly, const U& val) {
+    polynomial<FieldType> result(poly);
+    result += val;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator+(const U& val, const polynomial<FieldType>& poly) {
+    polynomial<FieldType> result(poly);
+    result += val;
+    return result;
+  }
+
+  // substraction:
+  template<typename FieldType>
+  inline polynomial<FieldType> operator-(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    polynomial<FieldType> result(left);
+    result -= right;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator-(const polynomial<FieldType>& poly, const U& val) {
+    polynomial<FieldType> result(poly);
+    result -= val;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator-(const U& val, const polynomial<FieldType>& poly) {
+    polynomial<FieldType> result(-poly);
+    result += val;
+    return result;
+  }
+
+  // multiplication:
+  template<typename FieldType>
+  inline polynomial<FieldType> operator*(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    polynomial<FieldType> result(left);
+    result *= right;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator*(const polynomial<FieldType>& poly, const U& val) {
+    polynomial<FieldType> result(poly);
+    result *= val;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator*(const U& val, const polynomial<FieldType>& poly) {
+    polynomial<FieldType> result(poly);
+    result *= val;
+    return result;
+  }
+
+  // division:
+  template<typename FieldType>
+  inline polynomial<FieldType> operator/(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    polynomial<FieldType> result(left);
+    result /= right;
+    return result;
+  }
+
+  template<typename FieldType, typename U>
+  inline polynomial<FieldType> operator/(const polynomial<FieldType>& poly, const U& val) {
+    polynomial<FieldType> result(poly);
+    result /= val;
+    return result;
+  }
+
+  // modulo:
+  template<typename FieldType>
+  inline polynomial<FieldType> operator%(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    polynomial<FieldType> result(left);
+    result %= right;
+    return result;
+  }
+
+  // comparisons:
+  template<typename FieldType>
+  bool operator==(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    if(left.size() != right.size())
+      return false;
+    typename std::vector<FieldType>::const_iterator left_it = left.begin();
+    typename std::vector<FieldType>::const_iterator right_it = right.begin();
+    while(left_it != left.end())
+      if(*left_it++ != *right_it++)
+        return false;
+    return true;
+  }
+
+  template<typename FieldType>
+  bool operator!=(const polynomial<FieldType>& left, const polynomial<FieldType>& right) {
+    if(left.size() != right.size())
+      return true;
+    typename std::vector<FieldType>::const_iterator left_it = left.begin();
+    typename std::vector<FieldType>::const_iterator right_it = right.begin();
+    while(left_it != left.end())
+      if(*left_it++ != *right_it++)
+        return true;
+    return false;
+  }
+
+  // special forms:
+  template<typename FieldType, typename Size>
+  polynomial<FieldType> chebyshev_form(const Size sz) {
+    std::vector<FieldType> tmp(sz + 1);
+    boost::math::tools::polynomial_chebyshev_form(tmp.begin(), tmp.end());
+    return boost::math::tools::polynomial<FieldType>(tmp);
+  }
+
+  template<typename FieldType, typename Size>
+  polynomial<FieldType> hermite_form(const Size sz) {
+    std::vector<FieldType> tmp(sz + 1);
+    boost::math::tools::polynomial_hermite_form(tmp.begin(), tmp.end());
+    return boost::math::tools::polynomial<FieldType>(tmp);
+  }
+
+  template<typename FieldType, typename Size>
+  polynomial<FieldType> laguerre_form(const Size sz) {
+    std::vector<FieldType> tmp(sz + 1);
+    boost::math::tools::polynomial_laguerre_form(tmp.begin(), tmp.end());
+    return boost::math::tools::polynomial<FieldType>(tmp);
+  }
+
+  template<typename FieldType, typename Size>
+  polynomial<FieldType> legendre_form(const Size sz) {
+    std::vector<FieldType> tmp(sz + 1);
+    boost::math::tools::polynomial_legendre_form(tmp.begin(), tmp.end());
+    return boost::math::tools::polynomial<FieldType>(tmp);
+  }
+
+  // output stream:
+  template<typename charT, typename traits, typename FieldType>
+  inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const polynomial<FieldType>& poly) {
+    os << "{";
+    typename std::vector<FieldType>::const_iterator it = poly.begin();
+    while(it != poly.end()) {
+      if(it != poly.begin())
+        os << ", ";
+      os << *it++;
+    }
+    os << "}";
+    return os;
+  }
+
+} // namespace tools
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP
+