...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/tools/polynomial.hpp>
namespace boost { namespace math { namespace tools { 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); template <class I> polynomial(I first, I last); template <class U> explicit polynomial(const U& point, typename std::enable_if<std::is_convertible<U, T> >::type* = nullptr); template <class Range> explicit polynomial(const Range& r, typename std::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = nullptr); // C++14 polynomial(std::initializer_list<T> l); // C++11 polynomial(std::vector<T>&& p); // access: size_type size()const; size_type degree()const; value_type& operator[](size_type i); const value_type& operator[](size_type i)const; std::vector<T> const& data() const; std::vector<T>& data(); explicit operator bool() const; polynomial<T> prime() const; polynomial<T> integrate() const; T operator()(T z) const; // modify: void set_zero(); // operators: template <class U> polynomial& operator +=(const U& value); template <class U> polynomial& operator -=(const U& value); template <class U> polynomial& operator *=(const U& value); template <class U> polynomial& operator /=(const U& value); template <class U> polynomial& operator %=(const U& value); template <class U> polynomial& operator +=(const polynomial<U>& value); template <class U> polynomial& operator -=(const polynomial<U>& value); template <class U> polynomial& operator *=(const polynomial<U>& value); template <class U> polynomial& operator /=(const polynomial<U>& value); template <class U> polynomial& operator %=(const polynomial<U>& value); }; template <class T> polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b); template <class T> polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b); template <class T> polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b); template <class T> polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b); template <class T> polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b); template <class T, class U> polynomial<T> operator + (const polynomial<T>& a, const U& b); template <class T, class U> polynomial<T> operator - (const polynomial<T>& a, const U& b); template <class T, class U> polynomial<T> operator * (const polynomial<T>& a, const U& b); template <class T, class U> polynomial<T> operator / (const polynomial<T>& a, const U& b); template <class T, class U> polynomial<T> operator % (const polynomial<T>& a, const U& b); template <class U, class T> polynomial<T> operator + (const U& a, const polynomial<T>& b); template <class U, class T> polynomial<T> operator - (const U& a, const polynomial<T>& b); template <class U, class T> polynomial<T> operator * (const U& a, const polynomial<T>& b); template <class T> polynomial<T> operator - (const polynomial<T>& a); template <class T> polynomial<T> operator >>= (const U& a); template <class T> polynomial<T> operator >> (polynomial<T> const &a, const U& b); template <class T> polynomial<T> operator <<= (const U& a); template <class T> polynomial<T> operator << (polynomial<T> const &a, const U& b); template <class T> bool operator == (const polynomial<T> &a, const polynomial<T> &b); template <class T> bool operator != (const polynomial<T> &a, const polynomial<T> &b); template <class T> polynomial<T> pow(polynomial<T> base, int exp); template <class charT, class traits, class T> std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly); template <typename T> std::pair< polynomial<T>, polynomial<T> > quotient_remainder(const polynomial<T>& a, const polynomial<T>& b); } // namespace tools }} // namespace boost { namespace math
This is a somewhat trivial class for polynomial manipulation.
See
Implementation is currently of the "naive" variety, with 𝑶(N2) multiplication, for example. This class should not be used in high-performance computing environments: it is intended for the simple manipulation of small polynomials, typically generated for special function approximation.
It does has division for polynomials over a field (here floating point, complex, etc) and over a unique factorization domain (integers). Division of polynomials over a field is compatible with Euclidean GCD.
Division of polynomials over a UFD is compatible with the subresultant algorithm for GCD (implemented as subresultant_gcd), but a serious word of warning is required: the intermediate value swell of that algorithm will cause single-precision integral types to overflow very easily. So although the algorithm will work on single-precision integral types, an overload of the gcd function is only provided for polynomials with multi-precision integral types, to prevent nasty surprises. This is done somewhat crudely by disabling the overload for non-POD integral types.
Advanced manipulations: the FFT, factorisation etc are not currently provided. Submissions for these are of course welcome :-)
First include the essential polynomial header (and others) to make the example:
#include <boost/math/tools/polynomial.hpp>
and some using statements are convenient:
using std::string; using std::exception; using std::cout; using std::abs; using std::pair; using namespace boost::math; using namespace boost::math::tools; // for polynomial using boost::lexical_cast;
Store the coefficients in a convenient way to access them, then create some polynomials using construction from an iterator range, and finally output in a 'pretty' formula format.
Tip | |
---|---|
Although we might conventionally write a polynomial from left to right in descending order of degree, Boost.Math stores in ascending order of degree. |
Read/write for humans: 3x^3 - 4x^2 - 6x + 10 Boost polynomial storage: [ 10, -6, -4, 3 ]
std::array<double, 4> const d3a = {{10, -6, -4, 3}}; polynomial<double> const a(d3a.begin(), d3a.end()); // With C++11 and later, you can also use initializer_list construction. polynomial<double> const b{{-2.0, 1.0}}; // formula_format() converts from Boost storage to human notation. cout << "a = " << formula_format(a) << "\nb = " << formula_format(b) << "\n\n";
for output:
a = 3x^3 - 4x^2 - 6x + 10 b = x - 2
// Now we can do arithmetic with the usual infix operators: + - * / and %. polynomial<double> s = a + b; cout << "a + b = " << formula_format(s) << "\n"; polynomial<double> d = a - b; cout << "a - b = " << formula_format(d) << "\n"; polynomial<double> p = a * b; cout << "a * b = " << formula_format(p) << "\n"; polynomial<double> q = a / b; cout << "a / b = " << formula_format(q) << "\n"; polynomial<double> r = a % b; cout << "a % b = " << formula_format(r) << "\n";
for output
a + b = 3x^3 - 4x^2 - 5x + 8 a - b = 3x^3 - 4x^2 - 7x + 12 a * b = 3x^4 - 10x^3 + 2x^2 + 22x - 20 a / b = 3x^2 + 2x - 2 a % b = 6
Division is a special case where you can calculate two for the price of one.
Actually, quotient and remainder are always calculated together due to the nature of the algorithm: the infix operators return one result and throw the other away.
If you are doing a lot of division and want both the quotient and remainder, then you don't want to do twice the work necessary.
In that case you can call the underlying function, quotient_remainder
,
to get both results together as a pair.
pair< polynomial<double>, polynomial<double> > result; result = quotient_remainder(a, b); // Reassure ourselves that the result is the same. BOOST_MATH_ASSERT(result.first == q); BOOST_MATH_ASSERT(result.second == r);
The source code is at polynomial_arithmetic.cpp