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