Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP
Chris@16: #define BOOST_MATH_TOOLS_POLYNOMIAL_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/math/tools/rational.hpp>
Chris@16: #include <boost/math/tools/real_cast.hpp>
Chris@16: #include <boost/math/special_functions/binomial.hpp>
Chris@16: 
Chris@16: #include <vector>
Chris@16: #include <ostream>
Chris@16: #include <algorithm>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace tools{
Chris@16: 
Chris@16: template <class T>
Chris@16: T chebyshev_coefficient(unsigned n, unsigned m)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    if(m > n)
Chris@16:       return 0;
Chris@16:    if((n & 1) != (m & 1))
Chris@16:       return 0;
Chris@16:    if(n == 0)
Chris@16:       return 1;
Chris@16:    T result = T(n) / 2;
Chris@16:    unsigned r = n - m;
Chris@16:    r /= 2;
Chris@16: 
Chris@16:    BOOST_ASSERT(n - 2 * r == m);
Chris@16: 
Chris@16:    if(r & 1)
Chris@16:       result = -result;
Chris@16:    result /= n - r;
Chris@16:    result *= boost::math::binomial_coefficient<T>(n - r, r);
Chris@16:    result *= ldexp(1.0f, m);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class Seq>
Chris@16: Seq polynomial_to_chebyshev(const Seq& s)
Chris@16: {
Chris@16:    // Converts a Polynomial into Chebyshev form:
Chris@16:    typedef typename Seq::value_type value_type;
Chris@16:    typedef typename Seq::difference_type difference_type;
Chris@16:    Seq result(s);
Chris@16:    difference_type order = s.size() - 1;
Chris@16:    difference_type even_order = order & 1 ? order - 1 : order;
Chris@16:    difference_type odd_order = order & 1 ? order : order - 1;
Chris@16: 
Chris@16:    for(difference_type i = even_order; i >= 0; i -= 2)
Chris@16:    {
Chris@16:       value_type val = s[i];
Chris@16:       for(difference_type k = even_order; k > i; k -= 2)
Chris@16:       {
Chris@16:          val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
Chris@16:       }
Chris@16:       val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
Chris@16:       result[i] = val;
Chris@16:    }
Chris@16:    result[0] *= 2;
Chris@16: 
Chris@16:    for(difference_type i = odd_order; i >= 0; i -= 2)
Chris@16:    {
Chris@16:       value_type val = s[i];
Chris@16:       for(difference_type k = odd_order; k > i; k -= 2)
Chris@16:       {
Chris@16:          val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
Chris@16:       }
Chris@16:       val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
Chris@16:       result[i] = val;
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class Seq, class T>
Chris@16: T evaluate_chebyshev(const Seq& a, const T& x)
Chris@16: {
Chris@16:    // Clenshaw's formula:
Chris@16:    typedef typename Seq::difference_type difference_type;
Chris@16:    T yk2 = 0;
Chris@16:    T yk1 = 0;
Chris@16:    T yk = 0;
Chris@16:    for(difference_type i = a.size() - 1; i >= 1; --i)
Chris@16:    {
Chris@16:       yk2 = yk1;
Chris@16:       yk1 = yk;
Chris@16:       yk = 2 * x * yk1 - yk2 + a[i];
Chris@16:    }
Chris@16:    return a[0] / 2 + yk * x - yk1;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: class polynomial
Chris@16: {
Chris@16: public:
Chris@16:    // typedefs:
Chris@16:    typedef typename std::vector<T>::value_type value_type;
Chris@16:    typedef typename std::vector<T>::size_type size_type;
Chris@16: 
Chris@16:    // construct:
Chris@16:    polynomial(){}
Chris@16:    template <class U>
Chris@16:    polynomial(const U* data, unsigned order)
Chris@16:       : m_data(data, data + order + 1)
Chris@16:    {
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial(const U& point)
Chris@16:    {
Chris@16:       m_data.push_back(point);
Chris@16:    }
Chris@16: 
Chris@16:    // copy:
Chris@16:    polynomial(const polynomial& p)
Chris@16:       : m_data(p.m_data) { }
Chris@16: 
Chris@16:    template <class U>
Chris@16:    polynomial(const polynomial<U>& p)
Chris@16:    {
Chris@16:       for(unsigned i = 0; i < p.size(); ++i)
Chris@16:       {
Chris@16:          m_data.push_back(boost::math::tools::real_cast<T>(p[i]));
Chris@16:       }
Chris@16:    }
Chris@16: 
Chris@16:    // access:
Chris@16:    size_type size()const { return m_data.size(); }
Chris@16:    size_type degree()const { return m_data.size() - 1; }
Chris@16:    value_type& operator[](size_type i)
Chris@16:    {
Chris@16:       return m_data[i];
Chris@16:    }
Chris@16:    const value_type& operator[](size_type i)const
Chris@16:    {
Chris@16:       return m_data[i];
Chris@16:    }
Chris@16:    T evaluate(T z)const
Chris@16:    {
Chris@16:       return boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size());;
Chris@16:    }
Chris@16:    std::vector<T> chebyshev()const
Chris@16:    {
Chris@16:       return polynomial_to_chebyshev(m_data);
Chris@16:    }
Chris@16: 
Chris@16:    // operators:
Chris@16:    template <class U>
Chris@16:    polynomial& operator +=(const U& value)
Chris@16:    {
Chris@16:       if(m_data.size() == 0)
Chris@16:          m_data.push_back(value);
Chris@16:       else
Chris@16:       {
Chris@16:          m_data[0] += value;
Chris@16:       }
Chris@16:       return *this;
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial& operator -=(const U& value)
Chris@16:    {
Chris@16:       if(m_data.size() == 0)
Chris@16:          m_data.push_back(-value);
Chris@16:       else
Chris@16:       {
Chris@16:          m_data[0] -= value;
Chris@16:       }
Chris@16:       return *this;
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial& operator *=(const U& value)
Chris@16:    {
Chris@16:       for(size_type i = 0; i < m_data.size(); ++i)
Chris@16:          m_data[i] *= value;
Chris@16:       return *this;
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial& operator +=(const polynomial<U>& value)
Chris@16:    {
Chris@16:       size_type s1 = (std::min)(m_data.size(), value.size());
Chris@16:       for(size_type i = 0; i < s1; ++i)
Chris@16:          m_data[i] += value[i];
Chris@16:       for(size_type i = s1; i < value.size(); ++i)
Chris@16:          m_data.push_back(value[i]);
Chris@16:       return *this;
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial& operator -=(const polynomial<U>& value)
Chris@16:    {
Chris@16:       size_type s1 = (std::min)(m_data.size(), value.size());
Chris@16:       for(size_type i = 0; i < s1; ++i)
Chris@16:          m_data[i] -= value[i];
Chris@16:       for(size_type i = s1; i < value.size(); ++i)
Chris@16:          m_data.push_back(-value[i]);
Chris@16:       return *this;
Chris@16:    }
Chris@16:    template <class U>
Chris@16:    polynomial& operator *=(const polynomial<U>& value)
Chris@16:    {
Chris@16:       // TODO: FIXME: use O(N log(N)) algorithm!!!
Chris@16:       BOOST_ASSERT(value.size());
Chris@16:       polynomial base(*this);
Chris@16:       *this *= value[0];
Chris@16:       for(size_type i = 1; i < value.size(); ++i)
Chris@16:       {
Chris@16:          polynomial t(base);
Chris@16:          t *= value[i];
Chris@16:          size_type s = size() - i;
Chris@16:          for(size_type j = 0; j < s; ++j)
Chris@16:          {
Chris@16:             m_data[i+j] += t[j];
Chris@16:          }
Chris@16:          for(size_type j = s; j < t.size(); ++j)
Chris@16:             m_data.push_back(t[j]);
Chris@16:       }
Chris@16:       return *this;
Chris@16:    }
Chris@16: 
Chris@16: private:
Chris@16:    std::vector<T> m_data;
Chris@16: };
Chris@16: 
Chris@16: template <class T>
Chris@16: inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result += b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result -= b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result *= b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class U>
Chris@16: inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result += b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class U>
Chris@16: inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result -= b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class U>
Chris@16: inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result *= b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class U, class T>
Chris@16: inline polynomial<T> operator + (const U& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(b);
Chris@16:    result += a;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class U, class T>
Chris@16: inline polynomial<T> operator - (const U& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(a);
Chris@16:    result -= b;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class U, class T>
Chris@16: inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
Chris@16: {
Chris@16:    polynomial<T> result(b);
Chris@16:    result *= a;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class charT, class traits, class T>
Chris@16: inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
Chris@16: {
Chris@16:    os << "{ ";
Chris@16:    for(unsigned i = 0; i < poly.size(); ++i)
Chris@16:    {
Chris@16:       if(i) os << ", ";
Chris@16:       os << poly[i];
Chris@16:    }
Chris@16:    os << " }";
Chris@16:    return os;
Chris@16: }
Chris@16: 
Chris@16: } // namespace tools
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP
Chris@16: 
Chris@16: 
Chris@16: