1 // (C) Copyright John Maddock 2006.
2 // (C) Copyright Jeremy William Murphy 2015.
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0. (See accompanying file
7 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
9 #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP
10 #define BOOST_MATH_TOOLS_POLYNOMIAL_HPP
16 #include <boost/assert.hpp>
17 #include <boost/config.hpp>
18 #ifdef BOOST_NO_CXX11_LAMBDAS
19 #include <boost/lambda/lambda.hpp>
21 #include <boost/math/tools/rational.hpp>
22 #include <boost/math/tools/real_cast.hpp>
23 #include <boost/math/policies/error_handling.hpp>
24 #include <boost/math/special_functions/binomial.hpp>
25 #include <boost/core/enable_if.hpp>
26 #include <boost/type_traits/is_convertible.hpp>
27 #include <boost/math/tools/detail/is_const_iterable.hpp>
32 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
33 #include <initializer_list>
36 namespace boost{ namespace math{ namespace tools{
39 T chebyshev_coefficient(unsigned n, unsigned m)
44 if((n & 1) != (m & 1))
52 BOOST_ASSERT(n - 2 * r == m);
57 result *= boost::math::binomial_coefficient<T>(n - r, r);
58 result *= ldexp(1.0f, m);
63 Seq polynomial_to_chebyshev(const Seq& s)
65 // Converts a Polynomial into Chebyshev form:
66 typedef typename Seq::value_type value_type;
67 typedef typename Seq::difference_type difference_type;
69 difference_type order = s.size() - 1;
70 difference_type even_order = order & 1 ? order - 1 : order;
71 difference_type odd_order = order & 1 ? order : order - 1;
73 for(difference_type i = even_order; i >= 0; i -= 2)
75 value_type val = s[i];
76 for(difference_type k = even_order; k > i; k -= 2)
78 val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
80 val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
85 for(difference_type i = odd_order; i >= 0; i -= 2)
87 value_type val = s[i];
88 for(difference_type k = odd_order; k > i; k -= 2)
90 val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
92 val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
98 template <class Seq, class T>
99 T evaluate_chebyshev(const Seq& a, const T& x)
101 // Clenshaw's formula:
102 typedef typename Seq::difference_type difference_type;
106 for(difference_type i = a.size() - 1; i >= 1; --i)
110 yk = 2 * x * yk1 - yk2 + a[i];
112 return a[0] / 2 + yk * x - yk1;
116 template <typename T>
122 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
123 * Chapter 4.6.1, Algorithm D: Division of polynomials over a field.
125 * @tparam T Coefficient type, must be not be an integer.
127 * Template-parameter T actually must be a field but we don't currently have that
128 * subtlety of distinction.
130 template <typename T, typename N>
131 BOOST_DEDUCED_TYPENAME disable_if_c<std::numeric_limits<T>::is_integer, void >::type
132 division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
134 q[k] = u[n + k] / v[n];
135 for (N j = n + k; j > k;)
138 u[j] -= q[k] * v[j - k];
142 template <class T, class N>
143 T integer_power(T t, N n)
148 return static_cast<T>(1u);
156 T result = integer_power(t, n / 2);
165 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
166 * Chapter 4.6.1, Algorithm R: Pseudo-division of polynomials.
168 * @tparam T Coefficient type, must be an integer.
170 * Template-parameter T actually must be a unique factorization domain but we
171 * don't currently have that subtlety of distinction.
173 template <typename T, typename N>
174 BOOST_DEDUCED_TYPENAME enable_if_c<std::numeric_limits<T>::is_integer, void >::type
175 division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
177 q[k] = u[n + k] * integer_power(v[n], k);
178 for (N j = n + k; j > 0;)
181 u[j] = v[n] * u[j] - (j < k ? T(0) : u[n + k] * v[j - k]);
187 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
188 * Chapter 4.6.1, Algorithm D and R: Main loop.
193 template <typename T>
194 std::pair< polynomial<T>, polynomial<T> >
195 division(polynomial<T> u, const polynomial<T>& v)
197 BOOST_ASSERT(v.size() <= u.size());
201 typedef typename polynomial<T>::size_type N;
203 N const m = u.size() - 1, n = v.size() - 1;
206 q.data().resize(m - n + 1);
210 division_impl(q, u, v, n, k);
214 u.normalize(); // Occasionally, the remainder is zeroes.
215 return std::make_pair(q, u);
219 // These structures are the same as the void specializations of the functors of the same name
220 // in the std lib from C++14 onwards:
225 T operator()(T const &x) const
233 template <class T, class U>
234 T operator()(T const &x, U const& y) const
242 template <class T, class U>
243 T operator()(T const &x, U const& y) const
249 } // namespace detail
252 * Returns the zero element for multiplication of polynomials.
255 polynomial<T> zero_element(std::multiplies< polynomial<T> >)
257 return polynomial<T>();
261 polynomial<T> identity_element(std::multiplies< polynomial<T> >)
263 return polynomial<T>(T(1));
266 /* Calculates a / b and a % b, returning the pair (quotient, remainder) together
267 * because the same amount of computation yields both.
268 * This function is not defined for division by zero: user beware.
270 template <typename T>
271 std::pair< polynomial<T>, polynomial<T> >
272 quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
274 BOOST_ASSERT(divisor);
275 if (dividend.size() < divisor.size())
276 return std::make_pair(polynomial<T>(), dividend);
277 return detail::division(dividend, divisor);
286 typedef typename std::vector<T>::value_type value_type;
287 typedef typename std::vector<T>::size_type size_type;
293 polynomial(const U* data, unsigned order)
294 : m_data(data, data + order + 1)
300 polynomial(I first, I last)
301 : m_data(first, last)
306 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
307 polynomial(std::vector<T>&& p) : m_data(std::move(p))
314 explicit polynomial(const U& point, typename boost::enable_if<boost::is_convertible<U, T> >::type* = 0)
317 m_data.push_back(point);
320 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
322 polynomial(polynomial&& p) BOOST_NOEXCEPT
323 : m_data(std::move(p.m_data)) { }
327 polynomial(const polynomial& p)
328 : m_data(p.m_data) { }
331 polynomial(const polynomial<U>& p)
333 m_data.resize(p.size());
334 for(unsigned i = 0; i < p.size(); ++i)
336 m_data[i] = boost::math::tools::real_cast<T>(p[i]);
339 #ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
340 template <class Range>
341 explicit polynomial(const Range& r, typename boost::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = 0)
342 : polynomial(r.begin(), r.end())
346 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
347 polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l))
352 operator=(std::initializer_list<T> l)
354 m_data.assign(std::begin(l), std::end(l));
362 size_type size() const { return m_data.size(); }
363 size_type degree() const
366 throw std::logic_error("degree() is undefined for the zero polynomial.");
367 return m_data.size() - 1;
369 value_type& operator[](size_type i)
373 const value_type& operator[](size_type i) const
378 T evaluate(T z) const
380 return this->operator()(z);
383 T operator()(T z) const
385 return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : T(0);
387 std::vector<T> chebyshev() const
389 return polynomial_to_chebyshev(m_data);
392 std::vector<T> const& data() const
397 std::vector<T> & data()
402 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
403 polynomial<T> prime() const
406 // Disable int->float conversion warning:
407 #pragma warning(push)
408 #pragma warning(disable:4244)
410 if (m_data.size() == 0)
412 return polynomial<T>({});
415 std::vector<T> p_data(m_data.size() - 1);
416 for (size_t i = 0; i < p_data.size(); ++i) {
417 p_data[i] = m_data[i+1]*static_cast<T>(i+1);
419 return polynomial<T>(std::move(p_data));
425 polynomial<T> integrate() const
427 std::vector<T> i_data(m_data.size() + 1);
428 // Choose integration constant such that P(0) = 0.
430 for (size_t i = 1; i < i_data.size(); ++i)
432 i_data[i] = m_data[i-1]/static_cast<T>(i);
434 return polynomial<T>(std::move(i_data));
438 polynomial& operator =(polynomial&& p) BOOST_NOEXCEPT
440 m_data = std::move(p.m_data);
444 polynomial& operator =(const polynomial& p)
451 typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator +=(const U& value)
459 typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator -=(const U& value)
467 typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator *=(const U& value)
469 multiplication(value);
475 typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator /=(const U& value)
483 typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator %=(const U& /*value*/)
485 // We can always divide by a scalar, so there is no remainder:
491 polynomial& operator +=(const polynomial<U>& value)
499 polynomial& operator -=(const polynomial<U>& value)
506 template <typename U, typename V>
507 void multiply(const polynomial<U>& a, const polynomial<V>& b) {
513 std::vector<T> prod(a.size() + b.size() - 1, T(0));
514 for (unsigned i = 0; i < a.size(); ++i)
515 for (unsigned j = 0; j < b.size(); ++j)
516 prod[i+j] += a.m_data[i] * b.m_data[j];
521 polynomial& operator *=(const polynomial<U>& value)
523 this->multiply(*this, value);
527 template <typename U>
528 polynomial& operator /=(const polynomial<U>& value)
530 *this = quotient_remainder(*this, value).first;
534 template <typename U>
535 polynomial& operator %=(const polynomial<U>& value)
537 *this = quotient_remainder(*this, value).second;
541 template <typename U>
542 polynomial& operator >>=(U const &n)
544 BOOST_ASSERT(n <= m_data.size());
545 m_data.erase(m_data.begin(), m_data.begin() + n);
549 template <typename U>
550 polynomial& operator <<=(U const &n)
552 m_data.insert(m_data.begin(), n, static_cast<T>(0));
557 // Convenient and efficient query for zero.
560 return m_data.empty();
563 // Conversion to bool.
564 #ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
565 typedef bool (polynomial::*unmentionable_type)() const;
567 BOOST_FORCEINLINE operator unmentionable_type() const
569 return is_zero() ? false : &polynomial::is_zero;
572 BOOST_FORCEINLINE explicit operator bool() const
574 return !m_data.empty();
578 // Fast way to set a polynomial to zero.
584 /** Remove zero coefficients 'from the top', that is for which there are no
585 * non-zero coefficients of higher degree. */
588 #ifndef BOOST_NO_CXX11_LAMBDAS
589 m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end());
591 using namespace boost::lambda;
592 m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
597 template <class U, class R>
598 polynomial& addition(const U& value, R op)
600 if(m_data.size() == 0)
602 m_data[0] = op(m_data[0], value);
607 polynomial& addition(const U& value)
609 return addition(value, detail::plus());
613 polynomial& subtraction(const U& value)
615 return addition(value, detail::minus());
618 template <class U, class R>
619 polynomial& addition(const polynomial<U>& value, R op)
621 if (m_data.size() < value.size())
622 m_data.resize(value.size(), 0);
623 for(size_type i = 0; i < value.size(); ++i)
624 m_data[i] = op(m_data[i], value[i]);
629 polynomial& addition(const polynomial<U>& value)
631 return addition(value, detail::plus());
635 polynomial& subtraction(const polynomial<U>& value)
637 return addition(value, detail::minus());
641 polynomial& multiplication(const U& value)
643 #ifndef BOOST_NO_CXX11_LAMBDAS
644 std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x * value; });
646 using namespace boost::lambda;
647 std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 * value));
653 polynomial& division(const U& value)
655 #ifndef BOOST_NO_CXX11_LAMBDAS
656 std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x / value; });
658 using namespace boost::lambda;
659 std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 / value));
664 std::vector<T> m_data;
669 inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
671 polynomial<T> result(a);
675 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
677 inline polynomial<T> operator + (polynomial<T>&& a, const polynomial<T>& b)
683 inline polynomial<T> operator + (const polynomial<T>& a, polynomial<T>&& b)
689 inline polynomial<T> operator + (polynomial<T>&& a, polynomial<T>&& b)
697 inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
699 polynomial<T> result(a);
703 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
705 inline polynomial<T> operator - (polynomial<T>&& a, const polynomial<T>& b)
711 inline polynomial<T> operator - (const polynomial<T>& a, polynomial<T>&& b)
717 inline polynomial<T> operator - (polynomial<T>&& a, polynomial<T>&& b)
725 inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
727 polynomial<T> result;
728 result.multiply(a, b);
733 inline polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b)
735 return quotient_remainder(a, b).first;
739 inline polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b)
741 return quotient_remainder(a, b).second;
744 template <class T, class U>
745 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (polynomial<T> a, const U& b)
751 template <class T, class U>
752 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (polynomial<T> a, const U& b)
758 template <class T, class U>
759 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (polynomial<T> a, const U& b)
765 template <class T, class U>
766 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator / (polynomial<T> a, const U& b)
772 template <class T, class U>
773 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator % (const polynomial<T>&, const U&)
775 // Since we can always divide by a scalar, result is always an empty polynomial:
776 return polynomial<T>();
779 template <class U, class T>
780 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (const U& a, polynomial<T> b)
786 template <class U, class T>
787 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (const U& a, polynomial<T> b)
793 template <class U, class T>
794 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (const U& a, polynomial<T> b)
801 bool operator == (const polynomial<T> &a, const polynomial<T> &b)
803 return a.data() == b.data();
807 bool operator != (const polynomial<T> &a, const polynomial<T> &b)
809 return a.data() != b.data();
812 template <typename T, typename U>
813 polynomial<T> operator >> (polynomial<T> a, const U& b)
819 template <typename T, typename U>
820 polynomial<T> operator << (polynomial<T> a, const U& b)
826 // Unary minus (negate).
828 polynomial<T> operator - (polynomial<T> a)
830 std::transform(a.data().begin(), a.data().end(), a.data().begin(), detail::negate());
835 bool odd(polynomial<T> const &a)
837 return a.size() > 0 && a[0] != static_cast<T>(0);
841 bool even(polynomial<T> const &a)
847 polynomial<T> pow(polynomial<T> base, int exp)
850 return policies::raise_domain_error(
851 "boost::math::tools::pow<%1%>",
852 "Negative powers are not supported for polynomials.",
853 base, policies::policy<>());
854 // if the policy is ignore_error or errno_on_error, raise_domain_error
855 // will return std::numeric_limits<polynomial<T>>::quiet_NaN(), which
856 // defaults to polynomial<T>(), which is the zero polynomial
857 polynomial<T> result(T(1));
860 /* "Exponentiation by squaring" */
870 template <class charT, class traits, class T>
871 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
874 for(unsigned i = 0; i < poly.size(); ++i)
888 // Polynomial specific overload of gcd algorithm:
890 #include <boost/math/tools/polynomial_gcd.hpp>
892 #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP